yufengdx.m

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

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

global dt df t f N T %定义全局变量
close all  %关闭以前的应用窗口

k=input('取样点数＝2^k, k=[13]'); %采样点数自定义，默认为2^15 
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个码元
t=linspace(-T/2+Ts/4,T/2+Ts/4,N);%时域横坐标
f=linspace(-Bs,Bs,N)+eps; %频域横坐标
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;

figure(1)
 set(1,'Position',[10,250,500,400])
                %设定窗口位置及大小 
                
 figure(2)
 set(2,'Position',[10,10,400,300])
                %设定窗口位置及大小                
                
figure(3)
 set(3,'Position',[515,250,500,400])
                %设定窗口位置及大小 
            
figure(4)
 set(4,'Position',[500,10,400,300])
                %设定窗口位置及大小                 
A=1;
f1=5/Ts;
f2=3/Ts;
m=A*cos(2*pi*f1*t);%z载波
n=-A*cos(2*pi*f2*t);


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

%s1=real(f2t(t2f(s).*GT)); 
s1=s.*m; 
s2=ss.*n;
st=s1+s2;%经调制后的FSK信号
ST=t2f(st);
Pa=ST.*conj(ST)/T;
EPa=(EPa*(loop1-1)+Pa)/loop1; %经调制后的FSK信号功率谱

n_ch=sqrt(sita)*randn(size(t));  %信道噪声
%sr=real(f2t(t2f(h).*S1));
stt=st+n_ch;    %加噪声后的信号
sr1=stt.*m; %解调
sr2=stt.*(-n);
SR=t2f(sr1);
Pr=SR.*conj(SR)/T;
EPr=(EPr*(loop1-1)+Pr)/loop1;

y1=real(f2t(t2f(sr1).*GR)); %通过接收滤波器
y2=real(f2t(t2f(sr2).*GR)); %通过接收滤波器
Y1=t2f(y1);
Py=Y1.*conj(Y1)/T;
EPy=(EPy*(loop1-1)+Py)/loop1;%通过发送滤波器信号的功率谱

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
 Pe(loop1)=n_err/(M*loop2);
for ii=1:Ts/dt;  
e(ii+[0:M-1]*L)=y;   %将输出信号变成双极性不归零码
end
E=t2f(e);
Po=E.*conj(E)/T;
EPo=(EPo*(loop1-1)+Po)/loop1;%输出信号功率谱
end

figure(1)
 subplot(2,2,1)
 plot(t,s,'r');
 axis([-T/64,T/64,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,e,'b');
 axis([-T/64,T/64,min(e)-0.5,max(e)+0.5])
 xlabel('t (ms)');
 ylabel('e  (V/KHz)');
 title('解调信号');
 
 subplot(2,2,4)
  bb=30+10*log10(EPo+eps);
 plot(f,bb,'g')
 axis([-20,20,-20,max(bb)+0.5])
 grid
 xlabel('f(Mhz)');
 ylabel('Po(f) (dBm/MHz)');
 title('输出信号功率谱密度图形');
 
 figure(2)
 subplot(1,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('FSK信号波形');
 
subplot(1,2,2)
 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('FSK信号功率谱密度图形');
 
 
 
 figure(3)
 subplot(2,3,1)
 plot(t,s1,'b');
 axis([-T/128,T/128,min(s1)-0.5,max(s1)+0.5])
 xlabel('t (ms)');
 ylabel('s1  (V/KHz)');
 title('1路调制信号波形');
 
 subplot(2,3,2)
 plot(t,sr1,'b');
 axis([-T/128,T/128,min(sr1)-0.5,max(sr1)+0.5])
 xlabel('t (ms)');
 ylabel('sr1  (V/KHz)');
 title('1路解调信号波形');
 
 subplot(2,3,3)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,y1(jj:jj+Na*L)); %接收信号眼图
  end
  
 subplot(2,3,4)
 plot(t,s2,'b');
 axis([-T/128,T/128,min(s2)-0.5,max(s2)+0.5])
 xlabel('t (ms)');
 ylabel('s2  (V/KHz)');
 title('2路调制信号波形');
 
 subplot(2,3,5)
 plot(t,sr1,'b');
 axis([-T/128,T/128,min(sr2)-0.5,max(sr2)+0.5])
 xlabel('t (ms)');
 ylabel('sr2  (V/KHz)');
 title('2路解调信号波形');
 
  subplot(2,3,6)
  tt=[0:dt:Na*L*dt];
  hold on
  for jj=1:Na*L:N-Na*L
    plot(tt,y2(jj:jj+Na*L)); %接收信号眼图
  end
   


  figure(4)
  semilogy(Eb_N0,Pe,'g');
  eb_n0=10.^(Eb_N0/10);
  hold on
  plot(Eb_N0,0.5*erfc(sqrt(eb_n0/2)));
  axis([0,14,1e-4,1])
  xlabel('Eb/N0')
  ylabel('Pe')