数字基带传输系统.m

数字基带传输系统的matlab仿真

%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('误码率曲线')