main1.m

自己变得OFDM迭代信道估计程序

% 本程序用以仿真块状导频时不同信噪比条件下的误码率和均方误差

echo off;clear all;
close all;
clc;
fprintf( 'OFDM仿真\n') ;
tic
% ---------------------------------------------%
%                   参数定义                    %
% --------------------------------------------- %
% Initialize the parameters
Num_frame=1000;  %仿真OFDM符号数
Bits_per_frame=60;  %每符号比特数
NumSubc=128;   %   载波数
Numcp=NumSubc/4;  %cp数约为载波数的1/4
mentor_times=20;  %蒙氏仿真  最大次数
%------------------------------------------------%
Bits_Tx = floor(rand(Num_frame,Bits_per_frame)*2); 
% Generate the random binary stream for transmit test
%------------------------------------------------%
% 信道编码   卷积码  对符号(即每行)进行编码
for i=1:Num_frame
    a=Bits_Tx(i,:);
    Bits_convenc(i,:)=cnv(a);   
end
%------------------------------------------------%
%  随机交织  对每符号进行交织   用系统函数生成
state=32;  %定义交织器初始参数
for i=1:Num_frame
    b=Bits_convenc(i,:);
Bits_Interleaved(i,:)=randintrlv(b,state); %随机交织
end
%-------------------------------------------------%
%  4QAM  星座图映射  对每帧进行映射
  for i=1:Num_frame 
      c=Bits_Interleaved(i,:);
      QAM_modulated(i,:)=QAM_modu(c); %调用映射函数 
  end
%----------------------------------------------------%
%  插入导频序列
[m,n]=size(QAM_modulated);
Pilot_seq=ones(1,n);  %导频序列 
Symbol_tx=[Pilot_seq;QAM_modulated];  %将其插入信息符号的前部， 成为发送符号 
%----------------------------------------------------%
%  串/并转换
 Symbol_paralleled=Symbol_tx.';
 %---------------------------------------------------%
 %  IFFT变换
 Symbol_ifft_temp=ifft( Symbol_paralleled,NumSubc);  %已变换
 %---------------------------------------------------%
 %  并/串转换
 Symbol_ifft= Symbol_ifft_temp.';
 %----------------------------------------------------%
 %  加cp前缀 
 [m,n]=size(Symbol_ifft);
 
 Symbol_cp=zeros(m,n+Numcp);
 Symbol_cp(:,(Numcp+1):(n+Numcp))=Symbol_ifft(:,1:n);%先把Symbol_ifft整体复制到Symbol_cp_temp的后半部分去
 Symbol_cp(:,1:Numcp)=Symbol_ifft(:,(n-Numcp+1):n);  %加cp  

 %-----------------------------------------------------%
%------------------------------------------------------%
 %  送入信道 
 h_time=[0.2901+0.6129j 0.4061+0.3881j 0.2103+0.1637j];  %假设的信道冲击响应
 L=length(h_time);  %信道长度
 H_act=fft([h_time zeros(1,NumSubc-L)],NumSubc);  %信道频率响应
 
 %-------------------------------------------------------%
 after_channel=zeros(m,n+Numcp);
 %  加高白噪声
 snrtable=zeros(11,2);
 MSE_table=zeros(11,1);
 for snr=0:2:20
     ber_ratio=0;  %初始化误码率
   snrtable(snr/2+1,1)=snr;
   estimation_error=zeros(1,NumSubc); %用以存储中间过程的估计误差 
   for mentor=1:mentor_times
 for i=1:m
     d=Symbol_cp(i,:);
 after_channel_temp=filter(h_time,[1],d);  %   过信道  用滤波器实现卷积
 after_channel(i,:)=awgn(after_channel_temp,snr,'measured');  %加噪声
 end
 %------------------------------------------------------% 
 %-----------------------------------------------------%
 % 去cp前缀
 Symbol_de_cp_temp1=after_channel.';  %串/并转换
 Symbol_de_cp=zeros(n,m);
 Symbol_de_cp(1:n,:)=Symbol_de_cp_temp1((Numcp+1):(n+Numcp),:);  %即要把Symbol_de_cp_temp1最上面的NumCP行去掉 变成n*m矩阵 已经作好串/并转换
 %-----------------------------------------------------%
 %  送入fft解调器
 Symbol_fft=fft(Symbol_de_cp,NumSubc);  %出来是并联的
 %-----------------------------------------------------%
 % 信道估计部分 
 Y1_Pilot_seq=Symbol_fft(:,1);  %接收到的导频序列
 X1_Pilot_seq=Pilot_seq(1,:).'; %这是发送端的导频序列
 H_p=zeros(NumSubc,1); 
 H_p(:)=Y1_Pilot_seq./X1_Pilot_seq; %由导频得到的初始信道估计
  H_temp1=ifft(H_p,NumSubc);  %反变换到时域
  H_temp2=H_temp1.';
  H_temp3=[H_temp2(:,1:L) zeros(1,NumSubc-L)];  %补零  
  H_est=fft(H_temp3,NumSubc).';  %变换到频域
  error=(abs(H_act-H_est.').^2);  %估计误差
  estimation_error(1,:)=estimation_error(1,:)+error;  %累计误差
%--------------------------------------------------------%
 for i=1:Num_frame     
     X(:,i)= Symbol_fft(:,i+1)./H_est;  %得到估计值
 end
 for i=1:Num_frame
     b=X(:,i).';
     transmit_bits_demapped(i,:)=QAM_demodu(b);  %解调
 end
 for i=1:Num_frame
     c=transmit_bits_demapped(i,:);
 transmit_bits_deinterleaved(i,:)=randdeintrlv(c,state);  %解交织
 end
  for i=1:Num_frame
      e=transmit_bits_deinterleaved(i,:);
      [transmit_bits_decoded,survivor_state,cumulated_metric]=viterbi(e);  %译码
      decode(i,:)=transmit_bits_decoded;
  end
  %-----------------------------------------------%
  %  误码率计算
  [number,ratio] = biterr(Bits_Tx,decode);
ber_ratio=ber_ratio+ratio;
   end
 snrtable(snr/2+1,2)=ber_ratio/mentor_times;
 estimation_error=estimation_error./mentor_times;
 MSE_table(snr/2+1,1)=sum(estimation_error(1,:))/NumSubc; %赋值
 
%---------------------------------
%各信噪比条件下估计值与真实值的比较
  figure(snr/2+1);
 plot(abs(H_act),'k-');
 grid on;
hold on;
 plot(abs(H_est),'k-.');
 legend('真实响应','估计的响应');
 title('拟合效果');
 hold off;
%-----------------------------------
 end

snrtable
 MSE_table
%%-----------------------------------------------%
figure;
semilogy(snrtable(:,1),snrtable(:,2),'k*-');
grid on;
xlabel('SNR/dB');
ylabel('BER');
title('信噪比/误码率曲线')

figure;
semilogy(snrtable(:,1),MSE_table(:,1),'k*-');
xlabel('SNR/dB');
ylabel('MSE');
title('信噪比/均方误差曲线');
grid on;
time_of_sim = toc
echo on;