fsk.m

包括RZ

%2FSK 调制与解调及功率谱
%信源为单极性不归零码（NRZ）FSK 
global dt df t f N T %定义全局变量 
close all    %关闭以前的应用窗口 
N=2^13;        %总的取样点数
L=32;	%L 为每个码元的取样点数
M=N/L;		%M 码元总数 
Rb=2;	%码元速率 
Rb=2   %Mb/s 
Ts=1/Rb;		%码元宽度 
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,500,400])  
figure(3)
set(3,'Position',[515,250,500,400]) 
figure(4)
  set(4,'Position',[500,10,500,400])  

A=1; 
f1=5/Ts; 
f2=2/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=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));    %信道噪声
  sr1=st.*m; %解调

  sr2=st.*(-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);
axis([-T/64,T/64,min(s)-0.5,max(s)+0.5])
grid
xlabel('t (ms)');
ylabel('s    (V/KHz)');
title('单极性不归零码信号波形');

subplot(2,2,2)
bb=30+10*log10(EP+eps);
plot(f,bb)
axis([-5,5,-20,max(bb)+0.5])
grid
xlabel('f(Mhz)');
ylabel('Ps(f) (dBm/MHz)');
title('信源功率谱密度图形');

subplot(2,2,3)
plot(t,e);
axis([-T/64,T/64,min(e)-0.5,max(e)+0.5])
grid
xlabel('t (ms)');
ylabel('e    (V/KHz)');
title('解调信号');

subplot(2,2,4)
bb=30+10*log10(EPo+eps);
plot(f,bb)
axis([-5,5,-20,max(bb)+0.5])
grid
xlabel('f(Mhz)');
ylabel('Po(f) (dBm/MHz)');
title('输出信号功率谱密度图形');

figure(2)
subplot(1,2,1)
plot(t,st);
axis([-1,1,min(st)-0.5,max(st)+0.5])
grid
xlabel('t (ms)');
ylabel('st    (V/KHz)');
title('FSK 信号波形');

subplot(1,2,2)
bb=30+10*log10(EPa+eps);
plot(f,bb)
axis([-15,15,-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])
grid
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])
grid
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 
title('1 路接收信号眼图');
grid

subplot(2,3,4)
plot(t,s2);
axis([-T/128,T/128,min(s2)-0.5,max(s2)+0.5])
grid
xlabel('t (ms)');
ylabel('s2    (V/KHz)');
title('2 路调制信号波形');

subplot(2,3,5)
plot(t,sr1);
axis([-T/128,T/128,min(sr2)-0.5,max(sr2)+0.5])
grid
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
title('2 路接收信号眼图');
grid
    
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')
title('误比特率曲线')