mine16.m

16QAM仿真的程序

clear
global dt df N t f T
QAM=16
MASK=sqrt(QAM);
Ntotal=2^16;            %总采样点数
N=Ntotal/2 ;       %每支路采样点数
L=16;                   %每码元的采样点数
M=N/L;                  %每支路码元数            
Rb=2;                   %码速率是2Mb/s
Tb=1/Rb;                %码元间隔
Ts=2/Rb;                %每一支路码元间隔
dt=Ts/L;                %每一支路时域采样间隔
df=1/(N*dt);            %每一支路频域采样间隔
T=N*dt;                 %每一支路截短时间
Bs=N*df/2;              %每一支路系统带宽
Na=4;                   %示波器扫描宽度为4个码元
Again=20;

alpha=1;                %alpha变化参数，可方便修改
t=[-T/2+dt/2:dt:T/2];   %每一支路时域横坐标
f=[-Bs+df/2:df:Bs];     %每一支路频域横坐标
ttotal=[-T/2+dt/2:dt/2:T/2];
%-----------------------------------------------设置载波
fc=5*Rb;                %载波频率
carry1=sqrt(2)*cos(2*pi*fc*t);
carry2=-sqrt(2)*sin(2*pi*fc*t);
%-----------------------------------------------升余弦
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;
%------------------------------------------------对10个信噪比
for loop1=1:10
    Eb_N0(loop1)=loop1-1;
    eb_n0(loop1)=10^(Eb_N0(loop1)/10);
    Eb=1;	
    n0=Eb/eb_n0(loop1);         %信道的噪声谱密度   
    sita=n0*Bs;                 %信道中噪声功率
%------------------------------------------------产生模拟信号
    for loop2=1:5
        for i=1:M
            temp=rand;
            desource1(i)=1+floor(MASK*temp);
            temp=rand;
            desource2(i)=1+floor(MASK*temp);
        end
        d=1;
        analog=[-3*d,-d,d,3*d];
        for i=1:M
            AOut1(i)=analog(:,desource1(i));  %第一路模拟信号
            AOut2(i)=analog(:,desource2(i));  %第二路模拟信号
        end
%------------------------------------------------对信号进行抽样
        imp1=zeros(1,N);            
        imp2=zeros(1,N); 
        imp1(L/2:L:N)=AOut1/dt;     %生冲激序列1
        imp2(L/2:L:N)=AOut2/dt;     %生冲激序列2
        IMP1=t2f(imp1);
        IMP2=t2f(imp2);
%------------------------------------------------经过发送滤波器合并为一路信号进行信道传输
        message1=real(f2t(GR.*IMP1));
        message2=real(f2t(GR.*IMP2));
        message=[message1,message2];
        sendM1=message1.*carry1;  
        sendM2=message2.*carry2;
        sendM=sendM1+sendM2;
%------------得到接收端输出噪声,不同支路的噪声信号独立随机产生，但是他们的分布和方差是一致的
        noise1=sqrt(sita)*randn(1,N);               %产生的随机噪声
        noise2=sqrt(sita)*randn(1,N);
        noiseOut1=real(f2t(t2f(noise1).*GR));       %经过接收滤波器后的随机噪声
        noiseOut2=real(f2t(t2f(noise2).*GR));
%------------------------------------------------接收并解调
        receiveM1=sendM1+noise1;                    %接受噪声干扰后的信号
        receiveM2=sendM2+noise2;
        RECEIVEM1=t2f(receiveM1);
        RECEIVEM2=t2f(receiveM2);
        GETMessage1=t2f(receiveM1.*carry1).*GT;     %解调
        GETMessage2=t2f(receiveM2.*carry2).*GT;
        getMessage1=real(f2t(GETMessage1));
        getMessage2=real(f2t(GETMessage2));
        getMessage=[getMessage1,getMessage2];
     %y=receiveM;      %得到接收端输信号
%------------------------------------------------判决
        numoferr=0;%误比特的初始置0
        for i=1:M
            for j=1:MASK
                metrics1(j)=abs(getMessage1(-L/2+i*L)-analog(j));
                metrics2(j)=abs(getMessage2(-L/2+i*L)-analog(j));
            end
            [min_metric1 decis1(i)]=min(metrics1);
            [min_metric2 decis2(i)]=min(metrics2);
            if(decis1(i)~=desource1(i)),%出错误比特+1
                numoferr=numoferr+1;
            end
            if(decis2(i)~=desource2(i)),
                numoferr=numoferr+1;
            end
        end
    end
    Pe(loop1)=numoferr/(M*loop2);
%------------------------------------------------画星座图   
    if(loop1==7) %只观察 某一分贝下星座图        
        figure(1)
        plot(getMessage1(L/2:L:N),getMessage2(L/2:L:N),'o')
        grid on
        hold on 
        title('16QAM星座图');
        for i=1:M
            if(decis1(i)~=desource1(i)|decis2(i)~=desource2(i)),
                plot(getMessage1(-L/2+L*i),getMessage2(-L/2+L*i),'rx');
                plot([getMessage1(-L/2+L*i), AOut1(i)],[getMessage2(-L/2+L*i), AOut2(i)],'--r');
            end
        end
        legend('正确判定点','错误判定点');
    end

%------------------------------------------------模拟发送基带信号
    SumP=0;
    if (loop1==10)
        SENDM1=t2f(sendM1);
        SumP=SumP+SENDM1.*conj(SENDM1)/T; 
        
        figure(2);
        set(2,'Position',[100,50,600,600])
                %设定窗口位置及大小
        subplot(6,1,1);
        plot(ttotal/Tb,message);
        axis([-50,50,-5,5]);
        legend('总发基带信号');
        subplot(6,1,2);
        plot(ttotal/Tb,getMessage,'g');
        axis([-50,50,-5,5]);
        legend('总收基带信号');
        subplot(6,1,3);  
        plot(t/Tb,message1,'b');
        axis([-50,50,-5,5]);
        legend('第1路发送基带信号');
        subplot(6,1,4)
        plot(t/Tb,getMessage1,'g'); 
        axis([-50,50,-5,5]);
        legend('第1路接收基带信号');
        subplot(6,1,5);  
        plot(t/Tb,message2,'b');
        axis([-50,50,-5,5]);
        legend('第2路发送基带信号');
        subplot(6,1,6); 
        plot(t/Tb,getMessage2,'g');
        axis([-50,50,-5,5]);
        legend('第2路接收基带信号');
%------------------------------------------------模拟发送调制信号    
        figure(3)
        subplot(4,1,1);  
        plot(t/Tb,sendM1);
        title('第1、2路发送\接收调制信号','fontsize',10);
        axis([-20,20,-5,5]);
        hold on
        %画包络
        plot(t/Tb,message1*sqrt(2),'--g','linewidth',2);    
        plot(t/Tb,-message1*sqrt(2),'--g','linewidth',2);  
        legend('第一路发送调制信号','第一路发送基带信号');
        subplot(4,1,2);  
        plot(t/Tb,receiveM1);
        axis([-20,20,-5,5]);
         hold on
        %画包络
        plot(t/Tb,getMessage1*sqrt(2),'--g','linewidth',2); 
        plot(t/Tb,-getMessage1*sqrt(2),'--g','linewidth',3); 
        legend('第一路接收调制信号','第一路接收基带信号');
        subplot(4,1,3);  
        plot(t/Tb,sendM2);
        axis([-20,20,-5,5]);
        hold on
        %画包络
        plot(t/Tb,message2*sqrt(2),'--g','linewidth',2);  
        plot(t/Tb,-message2*sqrt(2),'--g','linewidth',2);
        legend('第二路发送调制信号','第二路发送基带信号');
        subplot(4,1,4);  
        plot(t/Tb,receiveM2);
        axis([-20,20,-5,5]);
        hold on
        %画包络
        plot(t/Tb,getMessage2*sqrt(2),'--g','linewidth',2);  
        plot(t/Tb,-getMessage2*sqrt(2),'--g','linewidth',2);  
        legend('第二路接收调制信号','第二路接收基带信号');
%------------------------------------------------眼图  
        figure(4)
            tt=[0:dt:Na*L*dt];
            for ss=1:Na*L:N-Na*L
                plot(tt,getMessage1(ss:ss+Na*L));
                hold on             
            end
    end
end
figure(5) 
semilogy(Eb_N0,Pe);
hold on
semilogy(Eb_N0,(MASK-1)/MASK*0.5*erfc(sqrt(10.^(Eb_N0/10))),'g')
axis([0,9,1e-4,1])
title('16QAM误码率','fontsize',10);
xlabel('Eb/N0','fontsize',10);
ylabel('Pe','fontsize',10);
legend('16QAM实际误码率','16QAM理论误码率');
%------------------------------------------------功率普图
P=SumP/Again;
figure(6)
plot(f/Rb,30+10*log10(P));
title('16QAM功率普','fontsize',10);
xlabel('f/Rb','fontsize',10);
ylabel('P(f)','fontsize',10);