yufengcx.m

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

%信源为双极性不归零码（NRZ），取样点数自己定义 。。。。。QPSK


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

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个码元
alpha=0;   %升余弦滚降系数

t=linspace(-T/2,T/2,N); %频域横坐标
f=linspace(-Bs,Bs,N)+eps;%时域横坐标

%升余弦
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,250,500,400])
                %设定窗口位置及大小 
figure(2)
 set(2,'Position',[515,250,500,400])
                %设定窗口位置及大小 
                
figure(3)
 set(3,'Position',[200,10,500,400])
                %设定窗口位置及大小 
%。。。。。。。。。。。。。。两个载波
A=1;
fc=5/Ts;
m=A*cos(2*pi*fc*t);
n=-A*sin(2*pi*fc*t);  

EP=zeros(size(f));      %发送功率谱
EPa=zeros(size(f)); 
EPt=zeros(size(f));  
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)=q(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);           %信源信号的傅式变换
P=S.*conj(S)/T;     
EP=(EP*(loop1-1)+P)/loop1; %发送功率谱

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

st=Ia+Qa;  %并串转换
ST=t2f(st);           %信源信号的傅式变换
PT=ST.*conj(ST)/T;     
EPt=(EPt*(loop1-1)+PT)/loop1; %发送功率谱
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
Pe(loop1)=n_err/(M*loop2);
Y=t2f(y);           %信源信号的傅式变换
Pa=Y.*conj(Y)/T;      %发送功率谱
EPa=(EPa*(loop1-1)+Pa)/loop1;
end

figure(1)
 subplot(2,2,1)
 plot(t,s,'r');
 axis([-T/16,T/16,min(s)-0.5,max(s)+0.5])
 xlabel('t (ms)');
 ylabel('s  (V/KHz)');
 title('双极性不归零码信号波形');
 
 subplot(2,2,2)
 bb=30+10*log10(EP+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Ps(f) (dBm/MHz)');
 title('信源功率谱密度图形');
 
 subplot(2,2,3)
 plot(t,y,'b');
 axis([-T/16,T/16,min(y)-0.5,max(y)+0.5])
 xlabel('t (ms)');
 ylabel('y  (V/KHz)');
 title('输出信号波形');
 
 subplot(2,2,4)
 bb=30+10*log10(EPa+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Pa(f) (dBm/MHz)');
 title('输出功率谱密度图形');
 
 figure(2)
 subplot(2,3,1)
 plot(t,Ia,'b');
 axis([-T/256,T/256,min(Ia)-0.5,max(Ia)+0.5])
 xlabel('t (ms)');
 ylabel('Ia  (V/KHz)');
 title('并联调制信号波形(I)路');
 
 subplot(2,3,2)
 plot(t,Ir,'b');
 axis([-T/16,T/16,min(Ir)-0.5,max(Ir)+0.5])
 xlabel('t (ms)');
 ylabel('Ir  (V/KHz)');
 title('解调信号波形(I)路');
 
 subplot(2,3,3)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,Ir(jj:jj+Na*L)); 
  end
  title('解调信号波形眼图(I)路');

 subplot(2,3,4)
 plot(t,Qa,'b');
 axis([-T/256,T/256,min(Qa)-0.5,max(Qa)+0.5])
 xlabel('t (ms)');
 ylabel('Qa  (V/KHz)');
 title('并联调制信号波形(Q)路');
 
 subplot(2,3,5)
 plot(t,Qr,'b');
 axis([-T/16,T/16,min(Qr)-0.5,max(Qr)+0.5])
 xlabel('t (ms)');
 ylabel('Qr  (V/KHz)');
 title('解调信号波形(Q)路');
 
 subplot(2,3,6)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,Qr(jj:jj+Na*L)); 
  end
  title('解调信号波形眼图(Q)路');
  
  figure(3)
 subplot(2,2,1)
 plot(t,st,'b');
 axis([-T/256,T/256,min(st)-0.5,max(st)+0.5])
 xlabel('t (ms)');
 ylabel('st  (V/KHz)');
 title('调制信号波形');

 subplot(2,2,2)
 plot(t,ei,'b');
 axis([-T/16,T/16,min(ei)-0.5,max(ei)+0.5])
 xlabel('t (ms)');
 ylabel('ei  (V/KHz)');
 title('判决输出信号波形');
 
 subplot(2,2,3)
 bb=30+10*log10(EPt+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Pt(f) (dBm/MHz)');
 title('调制信号功率谱密度图形');
  
  subplot(2,2,4)
  semilogy(Eb_N0,Pe,'g');
  eb_n0=10.^(Eb_N0/10);
  hold on
  plot(Eb_N0,0.5*erfc(sqrt(eb_n0)));
  axis([0,14,1e-4,1])
  xlabel('Eb/N0')
  ylabel('Pe')
  title('误码率曲线')