qpsk.m

加性高斯白噪声信道条件下的QPSK信号的调制与解调

%本程序中时间单位是微秒
%频率单位为MHz
%码速率单位是Mb/s
global dt df t f N 
close all 
%初始化 
N=2^14;%总取样点数 
L=64;%每码元的采样点数 
M=N/L;%码元数 
Rb=2;%码率 
fc=4;%载频
Rt=1;%占空比 
Ts=0.5;%码元宽度 
dt=Ts/L;%系统时域采样间隔 
df=1/(N*dt);%系统频域采样间隔 
T=N*dt;%截断时间 
Bs=N*df/2; 
alpha=0; %理想低通

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,200,500,400]); 
figure(2) 
  set(2,'position',[520,250,300,200]); 
figure(3) 
  set(3,'position',[850,250,300,200]); 

    
    for loop1=1:10
         Eb_N0(loop1)=(loop1-1)  %Eb/N0 in dB
         eb_n0(loop1)=10^(Eb_N0(loop1)/10);
         Eb=1;	
         n0=Eb/(4*eb_n0(loop1)); %信道的噪声谱密度   
         sita=n0*Bs; %信道中噪声功率
         n_err=0;     %误码计数
           EP=zeros(size(f)); 
  
    for ii=1:10 
        a=sign(randn(1,M));
        ai=a(2:2:M)
        aq=a(1:2:M)
        i=zeros(1,N);    %产生冲激序列
        q=zeros(1,N);
        ss=zeros(1,N);
        for jj=1:Ts/dt
            ss(jj+[0:M-1]*L)=a;
        end
        for jj=1:2*Ts/dt
            i(jj+[0:(M-1)/2]*2*L)=ai;
            q(jj+[0:(M-1)/2]*2*L)=aq;
        end
     
      s1=cos(2*pi*fc*t)/sqrt(2);
      s2=sin(2*pi*fc*t)/sqrt(2);
      I=i.*s1;
      Q=q.*s2;
      qpsk=I+Q;
      QPSK=t2f(qpsk);
      P=QPSK.*conj(QPSK)/T;
       EP=(EP*(ii-1)+P)/ii; 
       n_ch=sqrt(sita)*randn(size(t)); %信道噪声
       r=qpsk+n_ch;
     %解调
     si=r.*s1;
     sq=r.*s2;
     yi=real(f2t(t2f(si).*GR));
     yq=real(f2t(t2f(sq).*GR));
     out1i=yi(L/2:L:N);		  %取样
      outi=sign(out1i);        %判决 
     out1q=yq(L/2:L:N);		  %取样
     outq=sign(out1q);        %判决 
 
    bi=outi(2:2:M);
    bq=outq(1:2:M);
    b=zeros(1,M)
    for jj=1:M/2
        b(2*jj)=bi(jj);
        b(2*jj-1)=bq(jj);end

     n_err=n_err+length(find(b~=a)) ;

 end
          Pe(loop1)=n_err/(M*ii);
           for jj=1:Ts/dt
              sout(jj+[0:M-1]*L)=b;
           end
end
    figure(1)
    subplot(2,2,1)
    plot(t,ss,'r') 
    grid 
    axis([-8,+8,-1.5,+1.5]); 
    xlabel('t (us)'); 
    ylabel('ss(t) (V)');
    title('输入双极性不归零码信号波形');
    subplot(2,2,2)
     plot(t,sout,'r') 
    axis([-8,+8,-1.5,+1.5]); 
    xlabel('t (us)'); 
    ylabel('sout(t) (V)');%作出QPSK信号时域波形图
    title('输出双极性不归零码信号波形');  
    subplot(2,2,3)
    plot(t,i,'r') 
    axis([-8,+8,-1.5,+1.5]); 
    xlabel('t (us)'); 
    ylabel('i(t) (V)');%作出I路信号时域波形图
    title('串并变换后第I路信号时域波形');  
    subplot(2,2,4)
    plot(t,q,'r') 
    grid 
    axis([-8,+8,-1.5,+1.5]); 
    xlabel('t (us)'); 
    ylabel('q(t) (V)');%作出Q路信号时域波形图
   title('串并变换后第Q路信号时域波形');  
        
   figure(2) 
    aa=30+10*log10(EP+eps); 
    plot(f,aa,'b'); 
    grid 
    axis([-8,+8,-50,50]); 
    xlabel('f (MHz)'); 
    ylabel('Ps(f)  (dbm/MHz)');%作出功率谱密度随频率变化图
    title('QPSK功率谱密度随频率变化图');
     figure(3) 
     semilogy(Eb_N0,Pe,'r');
    eb_n0=10.^(Eb_N0/10);
     hold on
    semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)));
     axis([0,10,1e-4,100])
     xlabel('Eb/N0')
     ylabel('Pe')
     title('误码率曲线图');