ofdm_sim.m

OFDM系统仿真

%%%%%%%%%%%  Developed for 16 QAM %%%%%%%%%%%%%%
clear all
close all
t0=clock;
%%%%%%%  输入参数 %%%%%%%%%%%
N=1024;%FFT点数
M=512;%OFDM子载波数目
L=10;%一帧的OFDM符号数目
MMMa=16;% QPSK 16QAM
MMM1a=log2(MMMa);%每个子载波可以携带的比特数目
Tg=0.1;%保护间隔
W=10;%%% 带宽(MHz)
DFF=W/M;%每个子载波的带宽
%%%%%% ts=, Effective Symbol Duration (Micro Sec) %%%%%%%
ts=M/W;%有效符号周期，即未加保护间隔的符号周期长度.(Micro Sec)
D=ts*Tg;%保护间隔长度
Ts=ts+D;%加保护间隔后的符号长度
st1=ts/N;%进行过采样后采样点之间的时域间隔.(Micro Sec)
Ng=round(D/(ts/N));%保护间隔内的采样点数
Nt=(N+Ng);%一个符号的（含保护间隔）的采样数
NNt=Nt*L;%一帧的采样点数
tbw=1/st1;%采样频率
%%%%%%% tsw=, Total Allocated Signal Bandwidth (MHz) %%%%%%%
tsw=(1/ts*M);%总的可分配的信号带宽
ddf=tbw/NNt; %频率分辩率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=0.1;%时延扩展
ds=st1;%空间时延
NP=10;  %多径数目

%%% 设置每条路径的指数衰减常数和延时 
amp(1)=0.7854;
amp(2)=0.4861;
amp(3)=0.3009;
amp(4)=0.1863;
amp(5)=0.1153;
amp(6)=0.0714;
amp(7)=0.0442;
amp(8)=0.0273;
amp(9)=0.0169;
amp(10)=0.0105;
dts1(1:NP)=0:9;%延时
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cn=30;%%% dB
cn1=10^(-cn/10);
 
LLL=10;%%% 试验次数
for KLL=1:LLL;
 if mod(KLL,1)==0,KLL
  end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  zxn1=zeros(1,NP); %NP多径数目，
  zyn1=zeros(1,NP);
  zn1=zeros(1,NP);
  zn2=zeros(1,NP);
  phn=zeros(1,NP);%相位
  amp1=zeros(1,NP);%噪声幅度
  qd3=zeros(1,NNt);%NNT是一帧的采样点数
  %产生瑞利衰退
  %增益
  zxn1(1:NP)=randn(1,NP);%噪声实部
  zyn1(1:NP)=randn(1,NP);%噪声虚部
  zn1(1:NP)=zxn1(1:NP)+j*zyn1(1:NP);%复数噪声
  zn2(1:NP)=abs(zn1(1:NP))/sqrt(2);%噪声幅度
  %%% Phase
  phn(1:NP)=atan2(zyn1(1:NP),zxn1(1:NP));%相位
  amp1(1:NP)=zn2(1:NP).*amp(1:NP).*exp(j*phn(1:NP)); %总幅度变化＝指数衰减×噪声幅度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ya为发送，x3为接收
ya=zeros(L,M);%未经调制数据符号
x3=zeros(L,N);%调制数据符号(已过采样)
L1=N/2+1-M/2;%过采样开始位置
L2=N/2+M/2;%过采样结束位置
%%%%%%%%%%%%%数据生成 %%%%%%%%%%%%%%
ya(1:L,1:M)=randint(L,M,MMMa);
%%%%%%%%%%%%% 调制数据,生成星座点 %%%%%%%%%
x3(1:L,L1:L2)=f_EnMQAM(ya(1:L,1:M),MMMa);%L1～L2之间为原始采样点，(0～L1)和(L2～Nff)之间补0
%%%%%%%%%%%频域信号经FFT变到时域信号%%%%%%%%%
 for k=1:L; 
 zz1(k,1:N)=ifft(fftshift(x3(k,1:N)))*sqrt(N);%fftshift使得过采样所补的0变到中间
 %增加保护间隔
 z3(k,1:Nt)=[zz1(k,N-Ng+1:N) zz1(k,1:N)];%把符号后面的Ng位移到前面，依次向后
 %z3为完成FFT过采样和加循环前缀后的信号
end
  q13=z3';
  q73=q13(:)';%变成行向量，完成并串变换
  %%%%%%%%%%在高斯白噪声的状态下进行信道估计
  q73a=fftshift(fft(q73)); %完成FFT，并把零频分量移到频谱中心
  pow1=q73a.*conj(q73a);%发送信号的功率
  pow2=10*log10(pow1/max(pow1));%衰落增益
  %%%%%%%%%%%
  po3=mean(q73.*conj(q73));%发送平均功率
  c3=po3/2*cn1; % 噪声功率
  %%%%%% 生成加性高斯白噪声  %%%%%
     zxc(1:L,1:Nt)=sqrt(c3).*randn(L,Nt);%噪声实部 
     zyc(1:L,1:Nt)=sqrt(c3).*randn(L,Nt);%噪声虚部
     znc(1:L,1:Nt)=zxc(1:L,1:Nt)+j*zyc(1:L,1:Nt);%复噪声
     znn1c=znc';
     znnc=znn1c(:)';% 变成行向量
     
     qd3(1:NNt)=amp1(1).*q73(1:NNt);%首先产生第一路多径
     %产生多径衰落信号
  for k=2:NP;
     qu23(k,1:NNt)=[q73(NNt+1-dts1(k):NNt) q73(1:NNt-dts1(k))];
     qu23(k,1:NNt)=amp1(k).*qu23(k,1:NNt);
     qd3(1:NNt)=qd3(1:NNt)+qu23(k,1:NNt); % NP条多径合成
  end
  %%%%%%%%%% 接收的多径信号功率谱 
  qd3a=fftshift(fft(qd3));
  pow3=qd3a.*conj(qd3a);
  pow4=10*log10(pow3/max(pow3));
  %%%% 加上高斯白噪声
  q73x=qd3+znnc;% 接收的经过多径信道的信号.
  q73y=q73+znnc; % 接收的经过白噪声信道的信号.
%%%%%%%%%%%%%%%%%%%%%%% 接收信号处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
imax=1;
%%% 去除保护间隔
for k=1:L;
     p93(k,1:Nt)=q73x(imax-1+(k-1)*Nt+(1:Nt)); % 串并变换
     p993(k,1:N)=[p93(k,Ng+1:N+Ng)]; % 去循环前缀
     %%%%%%%%%%  FFT %%%%%
     p813(k,1:N)=fftshift(fft(p993(k,1:N))/sqrt(N)); % 时域变到频域
     HH3(1:M)=p813(1,L1:L2)./x3(1,L1:L2);%？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？？
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     p8x3(k,1:M)=p813(k,L1:L2)./HH3(1:M);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     p93y(k,1:Nt)=q73y(imax-1+(k-1)*Nt+(1:Nt)); % 接收的多径未解调数据
     p993y(k,1:N)=[p93y(k,Ng+1:N+Ng)];
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
     p813y(k,1:N)=fftshift(fft(p993y(k,1:N))/sqrt(N));
     p8y3(k,1:M)=p813y(k,L1:L2);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
    
    %%% 解调
  pb1(1:L,1:M)=f_DeMQAM(p8x3(1:L,1:M),MMMa);%%%
  %%%理想误码率
  per1(KLL)=biterr(pb1(2:L,1:M),ya(2:L,1:M))/(M*(L-1)*MMM1a);
  %%% 多径衰落下的解调
  pb2(1:L,1:M)=f_DeMQAM(p8y3(1:L,1:M),MMMa);%%%
  per2(KLL)=biterr(pb2(2:L,1:M),ya(2:L,1:M))/(M*(L-1)*MMM1a);
end
%%%% Final BER
 ERR1=mean(per1(1:LLL));
 ERR2=mean(per2(1:LLL));
 
 %%%%%%%%%%%%%%%%%%%%%%%%
 figure(1)%%%%%% Information Real Data in Freq
 stem(real(x3(2,1:N)))
 grid on
 figure(2)%%%%%% Information Imag Data in Freq
 stem(imag(x3(2,1:N)))
 grid on
 figure(3)%%%%%% Time domain signal
 plot(abs(z3(1,1:Nt)))
 grid on
 figure(4)
 plot(abs(z3(1,1:Ng)))
 XX=abs(z3(1,1:Ng));
 grid on
 figure(5)
 plot(abs(z3(1,Nt-Ng:Nt)))
 YY=abs(z3(1,Nt-Ng:Nt));
 grid on
 figure(6)
 plot(pow2)
 grid on
 figure(7)
 plot(pow4)
 grid on
 figure(8)
 plot(abs(HH3))
 grid on
 figure(9)
 stem(real(p8x3(2,1:M)))
 grid on
 figure(10)
 stem(imag(p8x3(2,1:M)))
 grid on
 figure(11)
 plot(p8x3(2:L,1:M),'x')
 grid on
 figure(12)
 plot(p8y3(2:L,1:M),'x')
 grid on
 
 

etime(clock,t0)/60