global dt df t f N
close all
%初始化
N=2^16;%總取樣點數
L=64;%每碼元的采樣點數
M=N/L;%碼元數
Rb=2;%碼率
Rt=0.5;%占空比
Ts=0.5;%碼元寬度
dt=Ts/L;%系統時域采樣間隔
df=1/(N*dt);%系統頻域采樣間隔
T=N*dt;%截斷時間
Bs=N*df/2;

t=[-T/2+dt/2:dt:T/2];
f=[-Bs+df/2:df:Bs];

figure(1)
  set(1,'position',[10,50,300,200]);
figure(2)
  set(2,'position',[400,50,300,200]);
  
EP=zeros(size(f));
for ii=1:50
    signal_rz=round(rand(1,M));%產生信息碼元
    rz=zeros(1,N);%全零序列以備附值
    for jj=1:(Ts*Rt)/dt, rz(jj+[0:M-1]*L)=signal_rz;end%產生單極性RZ碼
    RZ=t2f(rz);%傅立葉變換
    P=RZ.*conj(RZ)/T;%求功率譜
    EP=(EP*(ii-1)+P)/ii;
    
    figure(1)
    aa=30+10*log10(EP+eps);
    plot(f,aa,'b');
    grid
    axis([-20,+20,-50,+50]);
    xlabel('f (MHz)');
    ylabel('Ps(f)  (dbm/MHz)');%作出功率譜密度隨頻率變化圖
    
    figure(2)
    plot(t,rz,'g')
    grid
    axis([-2,+2,-1.5,+1.5]);
    xlabel('t (us)');
    ylabel('s(t) (V)');%作出RZ信號時域波形圖
end