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;
alpha=input('滾降系數（默認[0.5]）');
if alpha==[],alpha=0.5;end
Na=4;

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

hr1=sin(pi*t/Ts)./(pi*t/Ts);
hr2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
hr=hr1.*hr2;  %升余弦濾波器沖擊響應
HR=abs(t2f(hr));%升余弦濾波器傳遞函數
GT=sqrt(HR);
GR=GT;

figure(1)
  set(1,'position',[10,50,300,200]);
figure(2)
  set(2,'position',[400,50,300,200]);

 
EP2=zeros(size(f))+eps;
for ii=1:100
    a=round(rand(1,M));%產生雙極性信息碼
    aa=2*a;
    imp=zeros(1,N);
    imp(L/2:L:N)=aa/dt;%產生沖擊序列
    IMP=t2f(imp);
     
    s=f2t(IMP.*HR);%沖擊序列經過升余弦滾降信道
    S=t2f(s);%傅立葉變換
    P2=S.*conj(S)/T;%產生功率譜密度
    EP2=(EP2*(ii-1)+P2+eps)/ii;
   end    
    figure(1)
    plot(f,30+10*log10(EP2),'r');    
    grid
    axis([-3,3,-100,50]);
    xlabel('f (MHz)');
    ylabel('Ps(f) (dBm/MHz)');
    title('發送信號升余弦滾降的功率譜');
   
    figure(2)
    tt=[0:dt:Na*L*dt];
    hold on
    for jjj=1:Na*L:N-Na*L
        plot(tt,s(jjj:jjj+Na*L));
    end
    title('升余弦滾降的眼圖');