pam.m

包括RZ

% 8pam通信系统仿真及误码率分析
global dt df t f N
close all
N=3*2^13;
L=16;
M=N/L;
Rb=1;		%码速率是Mb/s
Ts=1/Rb;        %码元间隔
dt=Ts/L;        %时域采样间隔
df=1/(N*dt) ;    %频域采样间隔
T=N*dt ;         %截短时间
Bs=N*df/2;       %系统带宽
alpha=0.5; 	%滚降系数
t=linspace(-T/2,T/2,N); %时域横坐标
f=linspace(-Bs,Bs,N)+eps;   %频域横坐标

figure(1)
set(1,'Position',[10,50,500,400])
figure(2)
set(2,'Position',[515,50,500,400])

%升余弦
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;
EP=zeros(size(f))+eps;

for loop1=1:10
 Eb_N0(loop1)=(loop1-1) ;%Eb/N0 in dB
 eb_n0(loop1)=10^(Eb_N0(loop1)/10);
 Eb=1;	
 n0=Eb/eb_n0(loop1); %信道的噪声谱密度   
 sita=n0*Bs; %信道中噪声功率
 n_err=0;  
 for jj=1:20
     
  % 8pam信号的产生
  a=round(rand(1,M));
  copy=ones(1,L);
  b=a(copy,:);
  c=reshape(b,1,L*M);

  for ii=1:M/3;
   A=[a(3*ii-2),a(3*ii-1),a(3*ii)];
   B=1-A;
   pa=A(1)*2^2+A(2)*2+A(3);
   pb=B(1)*2^2+B(2)*2+B(3);
   pc=pa-pb;
   a(3*ii-2:3*ii)=[pc,pc,pc];
  end

  d=a(copy,:);
  pam8=reshape(d,1,L*M);

  figure(1)
  subplot(4,1,1)
  plot(t,c)
  grid on
  axis([-10,+10,-.5,+1.5])
  xlabel('输入0，1码')
  title('0对应－7,1对应－5,2对应-3,3对应-1,4对应1,5对应3,6对应5,7对应7');

  subplot(4,1,2)
  plot(t,pam8)
  grid on
  xlabel('8PAM信号波形')
  axis([-10,+10,-7.1,+7.1])    
  % 8pam信号功率谱
  Pam8=abs(t2f(pam8));
  PAM=real(Pam8.*conj(Pam8)/T);
  EP=(EP*(jj-1)+PAM+eps)/jj;

  figure(1)
  subplot(4,1,4)
  plot(f,log10(EP+eps))
  grid on
  xlabel('8PAM信号的功率谱')
  axis([-2,+2,-1.5,2])
  
  %产生冲激序列,幅度为8pam
  imp=zeros(1,N);    
  imp(L/2:L:N)=a/dt;
  
  %发送信号
  St=t2f(imp).*GT;
  st=f2t(St);
  
  %噪声
  n=sqrt(sita)*randn(size(t));
  
  %接收及输出
  r=st+n;
  Y1=t2f(r).*GR;
  y1=f2t(Y1);
  
  %判决
  y2=y1(L/2:L:N);
  for k=1:length(y2);
   if y2(k)>=6,y2(k)=7;
   elseif abs(y2(k)-5)<1,y2(k)=5;
   elseif abs(y2(k)-3)<1,y2(k)=3;
   elseif abs(y2(k)-1)<1,y2(k)=1;
   elseif abs(y2(k)+1)<1,y2(k)=-1;
   elseif abs(y2(k)+3)<1,y2(k)=-3;
   elseif abs(y2(k)+5)<1,y2(k)=-5;
   elseif y2(k)<=-6,y2(k)=-7;
   end
  end
  y3=y2(copy,:);
  y=reshape(y3,1,L*M);

  figure(1)
  subplot(4,1,3)
  plot(t,y)
  grid on
  xlabel('输出信号波形')
  axis([-10,+10,-7.1,+7.1])
  
  %误码率
  n_err=n_err+length(find(y2~=a));
 end
 Pe(loop1)=n_err/(M*20);
end

figure(2) 
semilogy(Eb_N0,Pe,'g');
hold on
semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)),'r');
axis([0,9,1e-4,1])
xlabel('Eb/N0')
ylabel('Pe (dBm/MHz)')
title('误比特率曲线')