yufengx.m

yufengax-----基带 yufengbx-----PSK yufengcx-----QPSK yufengdx-----FSK yufengex-----16QAM yufengx-

%信源为单极性归零码，取样点数自己定义 
global dt df t f N T %定义全局变量
close all  %关闭以前的应用窗口
clear Eb_N0 Pe

k=input('取样点数＝2^k, k=[13]'); %采样点数自定义，默认为2^13
if k==[], k=13; end

%Rt=input('占空比[0：1](默认0.5,为1则是NRZ码)');
%if Rt==[],Rt=0.5;
%end;
N=2^k;    %总的取样点数
L=32;		%L为每个码元的取样点数
M=N/L;		%M码元总数	
Rb=2;		%码元速率Rb=2Mb/s
Ts=1/Rb;		%码元宽度Ts
dt=Ts/L;	%时域取样间隔
df=1/(N*dt);%频域取样间隔
T=N*dt;     %时间截断长度
Bs=N*df/2;  %带宽
Na=2;      %示波器扫描宽度为2个码元
t=linspace(-T/2,T/2,N); %频域横坐标
f=linspace(-Bs,Bs,N)+eps;%时域横坐标
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    BPSK

%升余弦
alpha=0;    %理想低通
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;

%。。。。。。。。。。。。。载波
A=1;
fc=1/Ts;
m=A*cos(2*pi*fc*t);
for loop1=1:15  
 Eb_N0(loop1)=2*(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;     %误码计数
 for loop2=1:5  
a=sign(randn(1,M));   %产生随机码
s=zeros(1,N);
for ii=1:Ts/dt;  
s(ii+[0:M-1]*L)=a;   %将信源变成双极性不归零码
end
S=t2f(s);           %信源信号的傅式变换
s1=s.*m; %经调制后的PSK信号
S1=t2f(s1);
n_ch=sqrt(sita)*randn(size(t));  %信道噪声

s2=s1+n_ch;    %加噪声后的信号
sr=s2.*m; %解调
SR=t2f(sr);
y0=real(f2t(t2f(sr).*GR)); %通过接收滤波器
Y0=t2f(y0);
y=y0(L/2:L:N);
y=sign(y);%抽样判决

n_err=n_err+length(find(y~=a));  
end
 Pe1(loop1)=n_err/(M*loop2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% QPSK

A=1;
fc=1/Ts;
m=A*cos(2*pi*fc*t);
n=-A*sin(2*pi*fc*t);  
for loop1=1:15    %误码率曲线
 Eb_N0(loop1)=2*(loop1-1);  %Eb/N0 in dB
 eb_n0(loop1)=10^(Eb_N0(loop1)/10);
 Eb=0.25;	
 n0=Eb/(eb_n0(loop1)); %信道的噪声谱密度   
 sita=n0*Bs; %信道中噪声功率
 n_err=0;     %误码计数

 for loop2=1:5  
 a=sign(randn(1,M));   %产生随机码
 i=zeros(1,M);         %串并转换
 q=zeros(1,M);         
y0=zeros(1,M);
s=zeros(1,N);
I=zeros(1,N);
Q=zeros(1,N);
y=zeros(1,N);
for nn=1:M/2          
    i(2*nn-1)=a(2*nn);
    q(2*nn-1)=a(2*nn-1);
    i(2*nn)=i(2*nn-1);
    q(2*nn)=i(2*nn-1);
end;   
for ii=1:Ts/dt;   
s(ii+[0:M-1]*L)=a;   %将信源变成不归零码
I(ii+[0:M-1]*L)=i;
Q(ii+[0:M-1]*L)=q;
end

S=t2f(s);           %信源信号的傅式变换
Ia=I.*m;    %两路分别调制
Qa=Q.*n;

st=Ia+Qa;  %并串转换

n_ch=sqrt(sita)*randn(size(t));  %信道噪声

st=st+n_ch; %加噪声
Iaa=st.*m;  %对两路信号分别解调
Qaa=st.*n;

Ir=real(f2t(t2f(Iaa).*GR)); %分别通过滤波器
Qr=real(f2t(t2f(Qaa).*GR)); 

%sr=Ir+Qr;
yi=Ir(L/2:L:N);  %抽样判决
yq=Qr(L/2:L:N);
yi=sign(yi);
yq=sign(yq);
bc=ones(1,L);        %定义复制的次数
yii=yi(bc,:);           %将a的第一行复制L次，生成L*M矩阵
ei=reshape(yii,1,L*M);  %将c重排成1*L*M数组
yqq=yq(bc,:);           %将a的第一行复制L次，生成L*M矩阵
eq=reshape(yqq,1,L*M);  %将c重排成1*L*M数组
for nn=1:M/2      %并串转换
    y0(2*nn)=yi(2*nn);
    y0(2*nn-1)=yq(2*nn-1);
end; 
for ii=1:L;  
                   
y(ii+[0:M-1]*L)=y0;  %变成不归零码

end
n_err=n_err+length(find(y0~=a));  
end
Pe2(loop1)=n_err/(M*loop2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   FSK

A=1;
f1=1/Ts;
f2=2/Ts;
m=A*cos(2*pi*f1*t);%z载波
n=-A*cos(2*pi*f2*t);
for loop1=1:15  
 Eb_N0(loop1)=2*(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;     %误码计数
 for loop2=1:5  
a=round(rand(1,M));   %产生随机码  
b=a-1;
s=zeros(1,N);
ss=zeros(1,N);
for ii=1:Ts/dt;  
s(ii+[0:M-1]*L)=a;   %将信源变成单极性不归零码
ss(ii+[0:M-1]*L)=b;
end
S=t2f(s);           %信源信号的傅式变换
s1=s.*m; 
s2=ss.*n;
st=s1+s2;%经调制后的FSK信号
ST=t2f(st);
n_ch=sqrt(sita)*randn(size(t));  %信道噪声
stt=st+n_ch;    %加噪声后的信号
sr1=stt.*m; %解调
sr2=stt.*(-n);
SR=t2f(sr1);
y1=real(f2t(t2f(sr1).*GR)); %通过接收滤波器
y2=real(f2t(t2f(sr2).*GR)); %通过接收滤波器
Y1=t2f(y1);
y11=y1(L/2:L:N);
y22=y2(L/2:L:N);
y=y11-y22;
y=(sign(y)+1)/2;;%抽样判决
n_err=n_err+length(find(y~=a));  
end
 Pe3(loop1)=n_err/(M*loop2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   16QAM
A=1;
fc=1/Ts;
m=A*cos(2*pi*fc*t);
n=-A*sin(2*pi*fc*t);      % lianggezaibo
for loop1=1:15    %误码率曲线
 Eb_N0(loop1)=2*(loop1-1);  %Eb/N0 in dB
 eb_n0(loop1)=10^(Eb_N0(loop1)/10);
 Eb=0.25;	
 n0=Eb/(eb_n0(loop1)); %信道的噪声谱密度   
 sita=n0*Bs; %信道中噪声功率
 n_err=0;     %误码计数

 for loop2=1:5  
 a=round(rand(1,M));   %产生随机码
 
for l3=1:L;  
s(l3+[0:M-1]*L)=a;   %将输出信号变成双极性不归零码
end
S=t2f(s);
 i=zeros(1,M/2);         %串并转换
 q=zeros(1,M/2);  
 I=zeros(1,N);         %串并转换
 Q=zeros(1,N); 
 
 i=a(1:2:M-1);
 q=a(2:2:M);
 ii=zeros(1,M/4);         %串并转换
 qq=zeros(1,M/4); 
 yii=zeros(1,M/2);         %串并转换
 yqq=zeros(1,M/2); 
 y=zeros(1,M); 
 for l1=1:2:(M/2-1);
        if i(l1)-i(l1+1)==1,
             ii((l1+1)/2)=1;
        elseif  i(l1)-i(l1+1)==-1
                ii((l1+1)/2)=-1;
        elseif  i(l1)==0
               ii((l1+1)/2)=-3;
           else  ii((l1+1)/2)=3;
        end
     
       if q(l1)-q(l1+1)==1,
             qq((l1+1)/2)=1;
        elseif  q(l1)-q(l1+1)==-1
                qq((l1+1)/2)=-1;
        elseif  q(l1)==0
               qq((l1+1)/2)=-3;
           else  qq((l1+1)/2)=3;
        end
 end

for l2=1:4*Ts/dt;   
   I(l2+[0:M/4-1]*4*L)=ii;
   Q(l2+[0:M/4-1]*4*L)=qq;
end;

Ia=I.*m;    %两路分别调制
Qa=Q.*n;

st=Ia+Qa;  %并串转换--矩形星座
ST=t2f(st);           %信源信号的傅式变换

n_ch=sqrt(sita)*randn(size(t));  %信道噪声

st=st+n_ch; %加噪声

Iaa=st.*m;  %对两路信号分别解调
Qaa=st.*n;

Ir=real(f2t(t2f(Iaa).*GR)); %分别通过滤波器
Qr=real(f2t(t2f(Qaa).*GR)); 

yi=Ir(2*L:4*L:N);  %抽样判决
yq=Qr(2*L:4*L:N);

for l2=1:2:M/2-1
    k1=abs(yi((l2+1)/2)+1);
    k2=abs(yi((l2+1)/2)+0.5);
    k3=abs(yi((l2+1)/2)-0.5);
    k4=abs(yi((l2+1)/2)-1);
    k=[k1,k2,k3,k4];
    
    p1=abs(yq((l2+1)/2)+1);
    p2=abs(yq((l2+1)/2)+0.5);
    p3=abs(yq((l2+1)/2)-0.5);
    p4=abs(yq((l2+1)/2)-1);
    p=[p1,p2,p3,p4];
    
    if k1==min(k)
        yii(l2)=0;yii(l2+1)=0;
    elseif k2==min(k)
        yii(l2)=0;yii(l2+1)=1;
    elseif k3==min(k)
        yii(l2)=1;yii(l2+1)=0;
    elseif k4==min(k)
        yii(l2)=1;yii(l2+1)=1;
    end
    
    
     if p1==min(p)
        yqq(l2)=0;yqq(l2+1)=0;
     elseif p2==min(p)
        yqq(l2)=0;yqq(l2+1)=1;
     elseif p3==min(p)
        yqq(l2)=1;yqq(l2+1)=0;
     elseif p4==min(p)
         yqq(l2)=1;yqq(l2+1)=1;
     end
    
 end
 y(1:2:M-1)=yii;
 y(2:2:M)=yqq;
n_err=n_err+length(find(y~=a));  
end
Pe4(loop1)=n_err/(M*loop2);
end
figure(1) 
  semilogy(Eb_N0,Pe1,'r');
  hold on
  semilogy(Eb_N0,Pe2,'y');
  hold on
  semilogy(Eb_N0,Pe3,'b');
  hold on
  semilogy(Eb_N0,Pe4,'g');
  axis([0,14,1e-4,1])
  xlabel('Eb/N0')
  ylabel('Pe')
  title('误码率曲线')