%Auther:楊登鋒
%Date:4/8/2008
%Describtion: 仿真數字基帶傳輸系統，包括輸入、輸出信號波形及其功率譜，
%                   眼圖（升余弦滾降系數a=0.5），Pe~Eb\No曲線，取樣時間無偏差

global dt df t f N T                            %定義全局變量
close all                                             %關閉以前的應用窗口
Rt=0.5;

%輸出提示信息
disp('請選擇輸入信號模式：');              
disp('1.單極性不歸零碼，');
disp('2.雙極性不歸零碼，');
disp('3.單極性歸零碼，');
disp('4.雙極性歸零碼.');
mode=input('mode=[1]:');                  %選擇模式
if isempty(mode), mode=1;end

if (mode==3) ||(mode==4)
    Rt=input('歸零碼的占空比,Rt=[0.5]:'); %Rt為占空比
            if isempty(Rt),Rt=0.5;end
end

k=input('取樣點數＝2^k, k=[14]');
if isempty(k), k=14; end
N=2^k;                                              %采樣點數
L=64;                                                 %每碼元的采樣點數
M=N/L;                                             %碼元數
Ts=0.5;                                              %碼元寬度是0.5us
dt=Ts/L;                                            %時域取樣間隔
df=1/(N*dt);                                      %頻域取樣間隔，單位MHz
T=N*dt;                                             %截短時間
Bs=N*df/2 ;                                       %系統帶寬
Na=2;                                               %示波器掃描寬度為2個碼元
alpha=0.5;                                         %升余弦滾降系數

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,250,500,400])        %設定窗口位置及大小 
figure(2)
 set(2,'Position',[515,250,500,400])      %設定窗口位置及大小                 
figure(3)
 set(3,'Position',[360,10,300,200])        %設定窗口位置及大小 
                
%功率譜矩陣，初值為零
EP=zeros(size(f));                               %信源功率譜
EPa=zeros(size(f));                             %信道前信號功率譜
EPr=zeros(size(f));                              %接收信號功率譜
EPo=zeros(size(f));                             %輸出信號功率譜

for loop1=1:20
    Eb_N0(loop1)=3*(loop1-1) ;          %Eb/N0 in dB，橫坐標
    eb_n0(loop1)=10^(Eb_N0(loop1)/10);%求其功率
    Eb=0.1;
    n0=Eb/eb_n0(loop1);                     %信道的噪聲譜密度
    sita=n0*Bs;                                    %信道中噪聲功率
    n_err=0;                                        %誤碼計數
    
    for ii=1:10

        switch mode
            case 1                                    %單極性不歸零碼
                rcode=round(rand(1,M));    %產生隨機碼
                bc=ones(1,L);                     %定義復制的次數
                c=rcode(bc,:);                     %將a的第一行復制L次，生成L*M矩陣
                e=reshape(c,1,L*M);           %將c重排成1*L*M數組
                s=e;
            case 2                                    %雙極性不歸零碼
                rcode=sign(randn(1,M));    %產生隨機碼
                bc=ones(1,L);                     %定義復制的次數
                c=rcode(bc,:);                     %將a的第一行復制L次，生成L*M矩陣
                e=reshape(c,1,L*M);           %將c重排成1*L*M數組
                s=e;
            case 3                                    %單極性歸零碼
                rcode=round(rand(1,M));    %產生隨機碼
                s=zeros(1,N);
                for jj=1:M
                    for tt=1:Rt*Ts/dt;
                        s(tt+(jj-1)*L)=rcode(jj);
                    end
                end
            case 4                                    %雙極性歸零碼
                rcode=sign(randn(1,M));    %產生隨機碼
                s=zeros(1,N);
                 for jj=1:M
                    for tt=1:Rt*Ts/dt;
                        s(tt+(jj-1)*L)=rcode(jj);
                    end
                end
               otherwise
                    disp('輸入錯誤！');
        end

        S=t2f(s);                                     %信源信號的傅式變換
        P=S.*conj(S)/T;                           %信源信號的功率譜
        EP=(EP*(loop1-1)+P)/loop1;       %信源功率譜

        s1=real(f2t(S.*GT));                     %經發送濾波后的信源信號
        S1=t2f(s1);
        Pa=S1.*conj(S1)/T;
        EPa=(EPa*(loop1-1)+Pa)/loop1; %信道前信號功率譜

        n_ch=sqrt(sita)*randn(size(t));   %信道噪聲
        nr=real(f2t(t2f(n_ch).*GR));        %輸出噪聲

        dtt=Ts/2;
        sr=real(f2t(S.*HR))+nr;             %接收信號
        SR=t2f(sr);
        Pr=SR.*conj(SR)/T;
        EPr=(EPr*(loop1-1)+Pr)/ loop1; %接收信號功率譜
        y=sr(L/2:L:N);                            %采樣
        
        %判決輸出信號
        if((mode==1)||(mode==3))
            aa=sign(y-0.25*Rt);                    
            aa=(aa+1)/2;                        
        else
            aa=sign(y-0.25*Rt);
        end
        
        n_err=n_err+length(find(aa~=rcode)); %誤碼率計數
    end
    
    Pe(loop1)=n_err/(M*ii);                 %誤碼率
    e=zeros(1,N);
    
        switch mode
            case 1                                   %單極性不歸零碼
                bc=ones(1,L);                    %定義復制的次數
                c=aa(bc,:);                         %將a的第一行復制L次，生成L*M矩陣
                e=reshape(c,1,L*M);          %將c重排成1*L*M數組
            case 2                                    %雙極性不歸零碼
                bc=ones(1,L);                     %定義復制的次數
                c=aa(bc,:);                         %將a的第一行復制L次，生成L*M矩陣
                e=reshape(c,1,L*M);           %將c重排成1*L*M數組
            case 3                                    %單極性歸零碼
                for jj=1:M
                    for tt=1:Rt*Ts/dt;
                        e(tt+(jj-1)*L)=aa(jj);
                    end
                end
            case 4                                    %雙極性歸零碼
                for jj=1:M
                    for tt=1:Rt*Ts/dt;
                        e(tt+(jj-1)*L)=aa(jj);
                    end
                end
               otherwise
                    disp('輸入錯誤！');
        end
    
    E=t2f(e);
    Po=E.*conj(E)/T;
    EPo=(EPo*(loop1-1)+Po)/loop1;     %輸出信號功率譜
end

%Figure 1
figure(1)
 subplot(2,2,1)
 plot(t,s,'r');
 axis([-T/16,T/16,min(s)-0.5,max(s)+0.5])
 xlabel('t (ms)');
 ylabel('s  (V/KHz)');
 title('信源信號波形');
 
 subplot(2,2,2)
 bb=30+10*log10(EP+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Ps(f) (dBm/MHz)');
 title('信源功率譜密度圖形');
 
 subplot(2,2,3)
 plot(t,e,'r');
 axis([-T/16,T/16,min(e)-0.5,max(e)+0.5])
 xlabel('t (ms)');
 ylabel('e  (V/KHz)');
 title('輸出信號波形');
 
 subplot(2,2,4)
 bb=30+10*log10(EPo+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Po(f) (dBm/MHz)');
 title('輸出信號功率譜密度');

%Figure 2
figure(2)
 subplot(2,3,1)
 plot(t,s1,'b');
 axis([-T/16,T/16,min(s1)-0.5,max(s1)+0.5])
 xlabel('t (ms)');
 ylabel('s1  (V/KHz)');
 title('信道前信號波形');
 
 subplot(2,3,2)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,s1(jj:jj+Na*L));                    %信道前信號眼圖
  end
  title('信道前升余弦滾降眼圖')
  
 subplot(2,3,3)
 bb=30+10*log10(EPa+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Pa(f) (dBm/MHz)');
 title('信道前信號功率譜密度圖形');
 
subplot(2,3,4)
 plot(t,sr,'b');
 axis([-T/16,T/16,min(sr)-0.5,max(sr)+0.5])
 xlabel('t (ms)');
 ylabel('sr  (V/KHz)');
 title('接收信號波形');
 
  subplot(2,3,5)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,sr(jj:jj+Na*L));                    %接收信號眼圖
  end
  title('信道后升余弦滾降眼圖')
  
 subplot(2,3,6)
 bb=30+10*log10(EPr+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Pr(f) (dBm/MHz)');
 title('接收信號功率譜密度圖形');

%Figure 3
figure(3) 
  semilogy(Eb_N0,Pe,'r:x');
  eb_n0=10.^(Eb_N0/10);
  hold on
  semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0/2)));
  axis([0,19,1e-4,1])
  xlabel('Eb/N0')
  ylabel('Pe')
  title('誤碼率曲線')