sim_8psk.m

通原实验课的作业

function [p]=sim_8psk(snr_in_dB)
E=1;
SNR=exp(snr_in_dB*log(10)/10); %信噪比
sgma=E/sqrt(2*SNR);	 %噪声的标准偏差
M=10000;  %码元数 
L=32;  %每个码元采样点数
N=M*L;  %采样点数
Rb=2;   %码速率
Ts=1/Rb;  %码元间隔
dt=Ts/L;  %时域采样间隔
df=1/(N*dt);  %频域采样间隔
d=16;  %半个码元采样点数
T=N*dt;  %截短时宽
Bs=N*df/2;  %频带宽度
t=[-T/2+dt/2:dt:T/2];   %时域坐标
t1=[-T/2+dt/2:32*dt:T/2];   %时域坐标
f=[-Bs+df/2:df:Bs];   %频域坐标

alpha=0.5;	%滚降系数
g1=sin(pi*t/Ts)./(pi*t/Ts);
g2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
g=g1.*g2;
G=sqrt(t2f(g,dt));  %发送和接收滤波器频谱
%Cf = 1;
Cf=1-0.5*exp(-2*i*pi*f*d);  %信道传输函数
Gf=G.^2.*Cf;  %系统总的频谱

s000=[1 0];
s001=[cos(pi/4) sin(pi/4)];
s011=[0 1];
s010=[cos(3*pi/4) sin(3*pi/4)];
s110=[-1 0];
s111=[cos(5*pi/4) sin(5*pi/4)];
s101=[0 -1];
s100=[cos(7*pi/4) sin(7*pi/4)];
%%  产生序列
c3 = (sign(randn(1,M))+1)/2;
c2 = (sign(randn(1,M))+1)/2;
c1 = (sign(randn(1,M))+1)/2;
for ii = 1:M
if((c1(ii)==0)&(c2(ii)==0)&(c3(ii)==0)),
      a(ii)=s000(1);
      b(ii)=s000(2);
   elseif((c1(ii)==0)&(c2(ii)==0)&(c3(ii)==1)),
      a(ii)=s001(1);
      b(ii)=s001(2);
   elseif((c1(ii)==0)&(c2(ii)==1)&(c3(ii)==0)),
      a(ii)=s010(1);
      b(ii)=s010(2);
   elseif((c1(ii)==0)&(c2(ii)==1)&(c3(ii)==1)),
      a(ii)=s011(1);
      b(ii)=s011(2);
   elseif((c1(ii)==1)&(c2(ii)==0)&(c3(ii)==0)),
      a(ii)=s100(1);
      b(ii)=s100(2);
   elseif((c1(ii)==1)&(c2(ii)==0)&(c3(ii)==1)),
      a(ii)=s101(1);
      b(ii)=s101(2);
   elseif((c1(ii)==1)&(c2(ii)==1)&(c3(ii)==0)),
      a(ii)=s110(1);
      b(ii)=s110(2);
   else
      a(ii)=s111(1);
      b(ii)=s111(2);
   end;
end
%%  a路经过发送端的成型滤波器和信道
imp=zeros(1,N);
imp(L/2:L:N)=a/dt;
S=t2f(imp,dt).*Gf;
s=f2t(S,dt);

s1=zeros(1,N);       %此段程序将经傅立叶变换得到的数字序列
s1(N/2+1:N)=s(1:N/2);%在中轴反转,如111-1-11,经变换变成-1-11111
s1(1:N/2)=s(N/2+1:N);%因此应把它反过来,这主要是由
s1=real(s1);         %f2t和t2f函数引起的.
gr=s1+sgma*randn(1,N);  %信号加噪声，实现SNR
gr1=gr(16:32:N);  %取每个码元中间的一个样本
%%  b路经过发送端的成型滤波器和信道
imp=zeros(1,N);
imp(L/2:L:N)=b/dt;
S=t2f(imp,dt).*Gf;
s=f2t(S,dt);

s1=zeros(1,N);       %此段程序将经傅立叶变换得到的数字序列
s1(N/2+1:N)=s(1:N/2);%在中轴反转,如111-1-11,经变换变成-1-11111
s1(1:N/2)=s(N/2+1:N);%因此应把它反过来,这主要是由
s1=real(s1);         %f2t和t2f函数引起的.
gr=s1+sgma*randn(1,N);  %信号加噪声，实现SNR
gr2=gr(16:32:N);  %取每个码元中间的一个样本

%%  判决和误码概率的计算
numofbiterror = 0;
for ii = 1:M
    r = [gr1(ii) gr2(ii)];
    
    c000=dot(r,s000);
   c001=dot(r,s001);
   c010=dot(r,s010);
   c011=dot(r,s011);
   c100=dot(r,s100);
   c101=dot(r,s101);
   c110=dot(r,s110);
   c111=dot(r,s111);
   
   c_max=max([c000 c001 c010 c011 c100 c101 c110 c111]);
   
   if(c000==c_max),
      decis1=0;decis2=0;decis3=0;
   elseif(c001==c_max),
      decis1=0;decis2=0;decis3=1;
   elseif(c010==c_max),
      decis1=0;decis2=1;decis3=0;
   elseif(c011==c_max),
      decis1=0;decis2=1;decis3=1;
   elseif(c100==c_max),
      decis1=1;decis2=0;decis3=0;
   elseif(c101==c_max),
      decis1=1;decis2=0;decis3=1;
   elseif(c110==c_max),
      decis1=1;decis2=1;decis3=0;
   else
      decis1=1;decis2=1;decis3=1;
   end;
   
   if(decis1~=c1(ii)),
      numofbiterror=numofbiterror+1;
   end;
   
   if(decis2~=c2(ii)),
      numofbiterror=numofbiterror+1;
   end;
   
   if(decis3~=c3(ii)),
      numofbiterror=numofbiterror+1;
   end;
end
   p=numofbiterror/(3*M);