%本程序中時間單位是微秒
%頻率單位為MHz
%碼元速率單位是Mb/s

global dt df t f N
close all

N=2^13;   %采樣點數
L=64;    %每碼元的采樣點數
M=N/L     %碼元數
Rb=2;
Ts=0.5		%碼元寬度是0.5us
dt=Ts/L;        
df=1/(N*dt)     %MHz
T=N*dt          %截短時間
Bs=N*df/2       %系統帶寬
n0=0.001;

t=[-T/2+dt/2:dt:T/2]; %時域橫坐標
f=[-Bs+df/2:df:Bs];   %頻域橫坐標

figure(1)
set(1,'Position',[10,50,400,400])
%設定窗口位置及大小 
figure(2)
set(2,'Position',[410,50,400,400])
%設定窗口位置及大小 
EP=zeros(size(f))+eps;
EPR=zeros(size(f))+eps;
for ii=1:30               
    a=round(rand(1,M));
    s=zeros(1,N);  
    for ii=1:0.5*L,
        s(ii+[0:M-1]*L)=a; ,
    end%產生單極性歸零碼
    
    sita=n0*Bs; %信道中噪聲功率
    n_ch=sqrt(sita)*randn(size(t)); %信道噪聲
    sr=s+n_ch;       %接收信號
    
    S=t2f(s);    %S是s的傅氏變換
    P=S.*conj(S)/T;    %功率譜
    SR=t2f(sr);    %接收信號傅氏變換
    PR=SR.*conj(SR)/T;    %功率譜
    EP=(EP*(ii-1)+P+eps)/ii;
    EPR=(EPR*(ii-1)+PR+eps)/ii;
end


figure(1)
subplot(2,1,1)
plot(t,s,'b')
grid
axis([-3,+3,-1.5,1.5'])
xlabel('t (ms)')
ylabel('s(t) (V)')

figure(1)
subplot(2,1,2)
aa=30+10*log10(EP+eps);    %加eps以避免除以零
plot(f,aa,'b');
grid
axis([-10,+10,-50,50])
xlabel('f  (MHz)')
ylabel('Ps(f)   (dBm/MHz)')


figure(2)
subplot(2,1,1)
plot(t,sr,'b')
grid
axis([-3,+3,-1.5,1.5'])
xlabel('t (ms)')
ylabel('s(t) (V)')

figure(2)
subplot(2,1,2)
aa=30+10*log10(EPR+eps);    %加eps以避免除以零
plot(f,aa,'b');
grid
axis([-10,+10,-50,50])
xlabel('f  (MHz)')
ylabel('Ps(f)   (dBm/MHz)')