work.m

基于MATLAB的8ASK数字调制系统设计的程序设计,别人编写的.

%本程序中时间单位是微秒
%频率单位为MHz
%码速率单位是Mb/s

global dt t f df N T
close all
clear Eb_N0 Pe

N1=3*2^14;                            %采样点数，为方便进行进制转换设为3的倍数
L=8;                                  %每码元的采样点数
M1=N1/L;                              %码元数

M=2*M1;                              %卷积编码后的码元数
N=L*M;                               %卷积编码后的总采样点数   

Rb=3;		                         %码速率是3Mb/s
Ts=1/Rb;                             %码元间隔
dt=Ts/L;                             %时域采样间隔
df=1/(N*dt);                         %频域采样间隔
T=N*dt;                              %截短时间
Bs=N*df/2;                           %系统带宽
T1=N1*dt;                            %截短时间

t=[-T/2+dt/2:dt:T/2];                %时域横坐标
f=[-Bs+df/2:df:Bs];                  %频域横坐标
t1=[-T1/2+dt/2:dt:T1/2];             %时域横坐标

figure(1)
 set(1,'Position',[10,500,400,300]); %设定窗口位置及大小 
figure(2)
 set(2,'Position',[600,500,400,300]);%设定窗口位置及大小 
figure(3)
 set(3,'Position',[10,50,400,300]);  %设定窗口位置及大小
figure(4)
 set(4,'Position',[600,50,400,300]); %设定窗口位置及大小
figure(5)
 set(5,'Position',[300,200,400,300]); %设定窗口位置及大小
 
FH=1;                                %截止频率
GR=LPF(FH);                          %低通滤波器 

%载波
A=1;
fc=3/Ts;

EP=zeros(size(f))+eps;
EPT=zeros(size(f))+eps;  
EPR=zeros(size(f))+eps; 
EPy=zeros(size(f))+eps; 

for loop1=1:20
    Eb_N0(loop1)=(loop1-1);             %Eb/N0 in dB
    eb_n0(loop1)=10^(Eb_N0(loop1)/10);
    Eb=1;	
    n0=Eb/eb_n0(loop1);                 %信道的噪声谱密度  
    sita=n0*Bs;                         %信道中噪声功率
    n_err=0;                            %误码计数
    for loop2=1:10  
        a_tem=round(rand(1,M1));             %产生随机码
        ss=zeros(1,N1);
        for loop3=1:Ts/dt;
            ss(loop3+[0:M1-1]*L)=a_tem;       %将信源变成单极性不归零码
        end;
        
        %卷积码编码     
        trel = poly2trellis(4,[11 13]);    % 定义g1=[1,0,1,1];g2=[1,1,1,1],约束长度为4
        a = convenc(a_tem,trel);           % 编码，编码后长度为2*N1即N.98304
 
        s=zeros(1,N);
        for loop3=1:Ts/dt;
            s(loop3+[0:M-1]*L)=a;       %将信号变成单极性不归零码
        end;
        
        % 产生八进制幅度序列  
        sa=zeros(1,N);      
        ii=zeros(1,M/3);   
        for l1=1:3:(M-1);
            if a(l1)==0&a(l1+1)==0&a(l1+2)==0,    %000=>+7
                ii((l1+2)/3)=+7;
            elseif a(l1)==0&a(l1+1)==0&a(l1+2)==1,%001=>+5
                ii((l1+2)/3)=+5;
            elseif a(l1)==0&a(l1+1)==1&a(l1+2)==0,%010=>+3
                ii((l1+2)/3)=+3;
            elseif a(l1)==0&a(l1+1)==1&a(l1+2)==1,%011=>+1
                ii((l1+2)/3)=+1;
            elseif a(l1)==1&a(l1+1)==0&a(l1+2)==0,%100=>-1
                ii((l1+2)/3)=-1;
            elseif a(l1)==1&a(l1+1)==0&a(l1+2)==1,%101=>-3
                ii((l1+2)/3)=-3;
            elseif a(l1)==1&a(l1+1)==1&a(l1+2)==0,%110=>-5
                ii((l1+2)/3)=-5;
            elseif a(l1)==1&a(l1+1)==1&a(l1+2)==1,%111=>-7
                ii((l1+2)/3)=-7;
            end;
        end;
        for l2=1:3*Ts/dt;
            sa(l2+[0:M/3-1]*3*L)=ii;    %八进制幅度序列
        end;
       
        S=t2f(s);
        P=S.*conj(S)/T;
        EP=(EP*(loop2-1)+P)/loop2;      %二进制信息序列的功率谱密度
        ST=t2f(sa);
        Pt=ST.*conj(ST)/T;
        EPT=(EPT*(loop2-1)+Pt)/loop2;   %八进制幅度序列的功率谱密度
      
        m=A*cos(2*pi*fc*t);
        s1=sa.*m;                       %经调制后的8ASK 信号        
        n_ch=sqrt(sita)*randn(size(t)); %信道噪声
        s2=s1+n_ch;                     %加噪声后的信号
        s3=s2.*m;                       %解调
        S3=t2f(s3);
        SR=S3.*GR;                      %接收信号的频谱
        Pr=SR.*conj(SR)/T;
        EPR=(EPR*(loop2-1)+Pr)/loop2;   %接收信号的功率谱密度  
        
        nr=real(f2t(t2f(n_ch.*m).*GR)); %输出噪声
        sr=real(f2t(SR));               %接收信号
         
        yi=sr(3*L/2:3*L:N);             %抽样判决
        yii=zeros(1,M);
        Se=zeros(1,M);
        
        for l3=1:3:M-1             %将抽样结果进行判决并译为对应的二进制序列
            k1=abs(yi((l3+2)/3)+4);
            k2=abs(yi((l3+2)/3)+3);
            k3=abs(yi((l3+2)/3)+1.8);
            k4=abs(yi((l3+2)/3)+0.6);
            k5=abs(yi((l3+2)/3)-0.6);
            k6=abs(yi((l3+2)/3)-1.8);
            k7=abs(yi((l3+2)/3)-3);
            k8=abs(yi((l3+2)/3)-4);
            k=[k1,k2,k3,k4,k5,k6,k7,k8];
            if k1==min(k)
                yii(l3)=1;yii(l3+1)=1;yii(l3+2)=1;   %-7=>111
            elseif k2==min(k)
                yii(l3)=1;yii(l3+1)=1;yii(l3+2)=0;   %-5=>110
            elseif k3==min(k)
                yii(l3)=1;yii(l3+1)=0;yii(l3+2)=1;   %-3=>101
            elseif k4==min(k)
                yii(l3)=1;yii(l3+1)=0;yii(l3+2)=0;   %-1=>100
            elseif k5==min(k)
                yii(l3)=0;yii(l3+1)=1;yii(l3+2)=1;   %+1=>011
            elseif k6==min(k)
                yii(l3)=0;yii(l3+1)=1;yii(l3+2)=0;   %+3=>010
            elseif k7==min(k)
                yii(l3)=0;yii(l3+1)=0;yii(l3+2)=1;   %+5=>001
            elseif k8==min(k)
                yii(l3)=0;yii(l3+1)=0;yii(l3+2)=0;   %+7=>000  
            end
        end 
                
        for l4=1:L;
            e(l4+[0:M-1]*L)=yii;   %判决后的二进制序列
        end
  
        Se=t2f(e);
        Py=Se.*conj(Se)/T;
        EPy=(EPy*(loop2-1)+Py)/loop2;      %判决后的信号功率谱
        n_err=n_err+length(find(yii~=a));                     
    
        
        % 维特比译码
        trel = poly2trellis(4,[11 13]); % 定义网格图结构
        tblen = 3;                      % Traceback length
        e_tem=vitdec(yii,trel,tblen,'cont','hard'); %%译码结果为e_tem，码长度为M1
        
        rec=zeros(1,N1);
        for loop4=1:Ts/dt;
            rec(loop4+[0:M1-1]*L)=e_tem;       %将信号变成单极性不归零码
        end;
      
    end   
     Pe(loop1)=n_err/(M*loop2);  %实际误码率    
end
        
figure(1)
   subplot(2,2,1)
        plot(t,s,'LineWidth',2.5);
        grid 
        axis([-10,+10,-2,2]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('卷积编码后的二进制信息序列');
   subplot(2,2,2)
        P1=30+10*log10(EP+eps);%加eps以避免除以零
        plot(f,P1);
        grid
        axis([-12,12,-20,max(P1)+0.5]);
        xlabel('f  (MHz)');
        ylabel('Ps(f)   (dBm/MHz)');
        title('二进制信息序列的功率谱');     
    subplot(2,2,3)
        plot(t,sa,'LineWidth',2.5);
        grid 
        axis([-10,+10,-8,8]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('八进制幅度序列');
    subplot(2,2,4)
        P2=30+10*log10(EPT+eps);%加eps以避免除以零
        plot(f,P2);
        grid
        axis([-12,12,-20,max(P2)+0.5]);
        xlabel('f  (MHz)');
        ylabel('Ps(f)   (dBm/MHz)');
        title('八进制幅度序列的功率谱');
    
figure(2)    
    subplot(2,2,1)  
        plot(t,s1,'LineWidth',2.5);
        grid 
        axis([-10,+10,-8,8]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('八进制幅度调制序列');
    subplot(2,2,2)
        plot(t,s3,'LineWidth',2.5);
        grid 
        axis([-10,+10,-8,8]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('八进制幅度解调序列') 
    subplot(2,2,3)  
        plot(t,sr,'LineWidth',2.5);
        grid 
        axis([-10,+10,-5,5]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('经低通后的接收信号');
    subplot(2,2,4) 
        P3=30+10*log10(EPR+eps);%加eps以避免除以零
        plot(f,P3);
        grid
        axis([-8,8,-20,max(P3)+0.5]);
        xlabel('f  (MHz)');
        ylabel('Ps(f)   (dBm/MHz)');
        title('经低通后的接收信号功率');
 
figure(3) 
    subplot(2,2,1)
        plot(t,s,'LineWidth',2.5);
        grid 
        axis([-10,+10,-2,2]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('卷积编码后的二进制信息序列');
    subplot(2,2,2)
        P4=30+10*log10(EP+eps);%加eps以避免除以零
        plot(f,P4);
        grid
        axis([-12,12,-20,30]);
        xlabel('f  (MHz)');
        ylabel('Ps(f)   (dBm/MHz)');
        title('二进制信息序列的功率谱');     
    subplot(2,2,3)  
        plot(t,e,'LineWidth',2.5);
        grid 
        axis([-10,+10,-2,2]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('信道传输判决后的二进制序列');
    subplot(2,2,4) 
        P5=30+10*log10(EPy+eps);%加eps以避免除以零
        plot(f,P5);
        grid
        axis([-12,12,-20,30]);
        xlabel('f  (MHz)');
        ylabel('Ps(f)   (dBm/MHz)');
        title('信道传输判决后的信号功率');       
        
figure(4) 
    subplot(2,1,1)
        plot(t1,ss,'LineWidth',2.5);
        grid 
        axis([-10,+10,-2,2]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('二进制输入信息序列');    
    subplot(2,1,2)  
        plot(t1,rec,'LineWidth',2.5);
        grid 
        axis([-10,+10,-2,2]);
        xlabel('t (us)');
        ylabel('s(t) (V)');
        title('维特比译码判决后的二进制序列');
            
 figure(5) 
    semilogy(Eb_N0,Pe,'g','LineWidth',2.5); %实际误码率曲线 
    eb_n0=10.^(Eb_N0/10);
    hold on
    semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0/2)),'LineWidth',2.5);%理论误码率曲线
    axis([0,12,1e-4,1]);
    xlabel('Eb/N0');
    ylabel('Pe');
    legend('实际值','理论值');
    title('系统误码率曲线');