ofdmchannelest4.m

块状导频的OFDM信道估计ofdmchannelest4为2008年本科毕业设计

%function ofdmchannelest4()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%对每帧OFDM增加循环前缀，确定性多径信道系数
%帧结构：前面1－4个训练符号，接着是（symbols_per_carrier＋1）个数据符号。
%该主程序OFDM信道估计,是对参考文献中程序修改后得到的。
%参考文献吕爱琴1 ,田玉敏1 ,朱明华05OFDM信道估计带MATLAB程序
%ofdmchannelest4.m 修改自ofdmchannelest3.m
%现版本作者Qingmin Meng tested on May4,2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc
fprintf('确定性多径信道系数下,OFDM信道估计仿真,如果仿真帧数Num_repeatsimulation_frame=1000，等10分钟!\n\n') ;
fpname = strcat('ofdmchannelest4result1.txt');
fp=fopen(fpname,'a+');fprintf(fp,'\n');
IFFT_bin_length=1024;
carrier_count=200;
bits_per_symbol=2; %选择调制阶数，如2 为4PSK或DQPSK
% symbols_per_carrier=1;
symbols_per_carrier=50;
TotalLength_training=1; %训练序列所形成的符号个数
% TotalLength_training=4;
%  fprintf('TotalLength_training必须为1 或者4,symbols_per_carrier必须为1 或者50.\n');
Lenght_CP=IFFT_bin_length/16;
Rayleigh_channel_model=1;
%------多径信道系数定义
% d1 = 4; a1=0.2; d2 =5; a2 =0.3; d3=6; a3 =0.4;d4 =7;a4 =0.5;
Path_number=2                                                            % Number of paths多径数1到3
delay=[0 2 4 6 9 13];                                                   % Delays of paths 多径时延
Path_Gain=[0.8084 0.462 0.253 0.259 0.0447 0.01];                         % Profile of channel model多径幅度
% Fc = 2.4e9;                                                              % Carrier frequency  载频
% %V =200;                                                                   % moving speed in km/h
% Tc = 1/5e6;	                                                          % Chip width  符号时间
% Time_Begin = 0;                                                           % Initializing the time
% Phase = 2*pi*rand(1,2);                                                   % Initializing the phase
channel_power=sum(abs(Path_Gain(1:Path_number)).^2);
Path_Gain_norm=Path_Gain(1:Path_number)/sqrt(channel_power);%功 率 归一 化
%------------------------------------主程序
Num_repeatsimulation_frame=1000;
Num_simulationSNR =5;%5 
for nEN=1:Num_simulationSNR
    SNR(nEN)=(nEN-1)*5;
    bit_error_count(nEN)=0;
  for frame=1:Num_repeatsimulation_frame    
    baseband_out_length=carrier_count*symbols_per_carrier*bits_per_symbol;
    carriers=(1:carrier_count) +(floor(IFFT_bin_length/4)-floor(carrier_count/2)); %非零子载波的序号1
    conjugate_carriers=IFFT_bin_length-carriers+2;%非零子载波的序号2，存放共轭符号
    % 信号发射
    baseband_out=round(rand(1,baseband_out_length)) ;%产生0，1序列
    convert_matrix=reshape(baseband_out,bits_per_symbol,length(baseband_out)/bits_per_symbol);
    for k=1:(length(baseband_out)/bits_per_symbol)
        integerbase(k)=0;
        for i=1:bits_per_symbol
            integerbase(k)=integerbase(k)+convert_matrix(i,k)*2^(bits_per_symbol-i);
        end;
    end;
    carrier_matrix=reshape(integerbase,carrier_count,symbols_per_carrier)';
    % QDPSK调制
    carrier_matrix =[zeros(1,carrier_count);carrier_matrix];%第一行为已知参考量，即全为零
    for i = 2:(symbols_per_carrier+1) %从第二行开始，第i行与第(i-1)行进行模4加（在QPSK时）
        carrier_matrix(i,:)=rem(carrier_matrix(i,:)+carrier_matrix(i-1,:),2^bits_per_symbol);
    end
    carrier_matrix = carrier_matrix*((2*pi)/(2^bits_per_symbol)) ;%0,1,2,3被映射为QPSK符号
    [X, Y]=pol2cart(carrier_matrix,ones( size(carrier_matrix, 1),size(carrier_matrix,2)));%把级坐标转化为直角坐标（x,y)
    complex_carrier_matrix=complex(X,Y) ;%使得（x,y)变为一个复数
    % 加训练序列1*200,
    training_symbols = [ 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j ...
    -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 ...
    -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j ...
    1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j ...
    -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j ...
    -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 ...
    -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 ...
    j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j ...
    -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 ...
    -j -j -1 1 j j 1 -1 -j -j -1 ];
    if TotalLength_training==1
        training_symbols = training_symbols;%
    elseif TotalLength_training==4
        training_symbols = cat(1, training_symbols,training_symbols);%然后级联形成更多行的训练序列矩阵
        training_symbols = cat(1, training_symbols,training_symbols);
    else
        fprintf('TotalLength_training必须为1 或者4.\n');
    end;
    complex_carrier_matrix=cat(1,training_symbols,complex_carrier_matrix);%训练序列矩阵与发射信号矩阵级联，训练序列加在前行
    IFFT_modulation =zeros(TotalLength_training+symbols_per_carrier+1,IFFT_bin_length); %在频域缓冲区，置0

    IFFT_modulation(:,carriers) =complex_carrier_matrix;%(TotalLength_training+symbols_per_carrier+1)*200矩阵， 当symbols_per_carrier=50时
    IFFT_modulation(:,conjugate_carriers)=conj(complex_carrier_matrix);%子载波共轭对称
    time_wave_matrix=ifft(IFFT_modulation.');%对每列(OFDM符号)进行ifft (长度IFFT_bin_length＝1024)

    for i = 1:TotalLength_training+symbols_per_carrier+1
        windowed_time_wave_matrix(:,i)=real(time_wave_matrix(:,i));%取实部,
    end;
    ofdm_symbol_cp=[windowed_time_wave_matrix(end-Lenght_CP+1:end,:);windowed_time_wave_matrix];%加CP
    % [num_row1,num_col1]=size(ofdm_symbol_cp);
    ofdm_modulation=reshape(ofdm_symbol_cp,...
        1,(IFFT_bin_length+Lenght_CP)*(TotalLength_training+symbols_per_carrier+1));
    % ofdm_modulation=reshape(windowed_time_wave_matrix,1,IFFT_bin_length*(TotalLength_training+symbols_per_carrier+1));

    Tx_data = ofdm_modulation;%1行，(TotalLength_training+symbols_per_carrier+1)*1024列,if TotalLength_training=1 and symbols_per_carrier=1
    % 信道
    length_signal0=(IFFT_bin_length+Lenght_CP)*(TotalLength_training+symbols_per_carrier+1);
    Signal_length=length_signal0+delay(Path_number);%
%     if Rayleigh_channel_model==1%多径信道用瑞利衰落信道
            s=Tx_data;
            if nEN==1 %Creat both the time domain and frequency domain channel coefficients
                H_tap_time=zeros(1,IFFT_bin_length);
                for i=1:Path_number
                    H_tap_time(delay(i)+1)=Path_Gain_norm(i);%
                end;%
                H_buffer(frame,:)=H_tap_time;
                H_f=fft(H_tap_time);%--完成准静态信道系数的fft变换
                H_fbuffer(frame,:)=H_f;  %频域信道信道冲击响应      
            else %recover the coefficients
                H_tap=squeeze(H_buffer(frame,:)); H_f=squeeze(H_fbuffer(frame,:)); 
            end ;% 
            CH_Data =Path_Gain_norm; %时域信道信道冲击响应或系数
            for p=1:Path_number
                ss(p,:)=[zeros(1,delay(p)) s zeros(1,delay(Path_number)-delay(p))];                            
            end ; 
            S=CH_Data*ss; %信道冲击响应与发射信号的卷积。注意多径信号相加，实际是采用移位相乘来实现卷积 .           
            Tx_data_multipath=S(1:length_signal0);
%     end;
    %-----------------------
    Tx_signal_power=var(Tx_data_multipath); %估计信号功率
    linear_SNR =10^(SNR(nEN)/10);
    noise_sigma=Tx_signal_power/linear_SNR;%求噪声功率
    noise_scale_factor = sqrt(noise_sigma);
    noise = randn(1,length(Tx_data_multipath))*noise_scale_factor;
    Rx_Data =Tx_data_multipath+noise; %加入AWGN
    % 信号接收
    Rx_Data_matrix_CP=reshape(Rx_Data,(IFFT_bin_length+Lenght_CP),TotalLength_training+symbols_per_carrier+1);
    Rx_Data_matrix=Rx_Data_matrix_CP(Lenght_CP+1:end,:);%除去每一列中的CP
    % Rx_Data_matrix=reshape(Rx_Data,IFFT_bin_length,TotalLength_training+symbols_per_carrier+1);

    Rx_spectrum = fft(Rx_Data_matrix);%每一列进行FFT
    Rx_carriers = Rx_spectrum(carriers,:).'; %--------提取对称子载波一部分，仅采用非零子载波的序号1

    Rx_training_symbols=Rx_carriers((1:TotalLength_training),:) ;%提取前四行训练序列Y

    Rx_carriers=Rx_carriers((TotalLength_training+1: end),:) ;%%--------提取5行以后的信号

    % 信道估计
    Rx_training_symbols = Rx_training_symbols./training_symbols;%Y除以已知的训练序列X，得到频域信道系数H，H为[h1 h2]的变换
     Rx_training_symbols_deno = abs(Rx_training_symbols).^2;    
%     Rx_training_symbols_deno = Rx_training_symbols.^2; %原文章中错误，缺少取模运算
    Rx_training_H_square=0;
    Rx_training_H=0;
    for i_train=1:TotalLength_training
        Rx_training_H_square=Rx_training_H_square+Rx_training_symbols_deno(i_train,:);
        Rx_training_H=Rx_training_H+Rx_training_symbols(i_train,:);
    end;
    Rx_training_symbols_deno=Rx_training_H_square;%前1-4行平均得到分母部分
    Rx_training_symbols_nume=conj(Rx_training_H);%前1-4行平均和共轭得到分子部分

    Rx_training_symbols0 = Rx_training_symbols_nume./Rx_training_symbols_deno;%分子除以分母，得到频域信道系数H的倒数，即1/Hk
    %--存储估计
                H_fbuffer_estimate(frame,:)=1./Rx_training_symbols0;%存储估计的频域信道系数,Hk             
                H_fbuffer_theory(frame,:)=H_f(carriers);           %存储理想的频域信道系数,Hk_theory  
    %--                
    Rx_training_symbols = repmat(Rx_training_symbols0,symbols_per_carrier+1,1); %如果一个训练序列所形成的符号，拷贝（symbols_per_carrier+1)次
    Rx_carriers = Rx_training_symbols.*Rx_carriers;%迫零运算，这里Rx_training_symbols含有1/Hk,Rx_carriers含有Hk*zk
    %-------------------对于QPSK，需要发现接收信号处于4个象限的哪一个象限
    Rx_phase = angle(Rx_carriers)*(180/pi) ;
    phase_negative = find(Rx_phase<0) ;
    Rx_phase(phase_negative) =rem(Rx_phase(phase_negative)+360,360);
    Rx_decoded_phase = diff(Rx_phase) ;%差分DIFF(X), for a vector X, is [X(2)-X(1)  X(3)-X(2) ... X(n)-X(n-1)].
    phase_negative = find(Rx_decoded_phase<0) ;
    Rx_decoded_phase(phase_negative)= rem(Rx_decoded_phase(phase_negative)+360,360);
    % QDPSK解调
    base_phase =360/2^bits_per_symbol;
    delta_phase =base_phase/2;%对于QPSK,角度为90度
    Rx_decoded_symbols =zeros(size(Rx_decoded_phase,1),size(Rx_decoded_phase,2));
    %
    for i = 1:(2^bits_per_symbol-1)
        center_phase = base_phase*i;
        plus_delta =center_phase+delta_phase;
        minus_delta =center_phase-delta_phase;
        decoded =find((Rx_decoded_phase<=plus_delta) &(Rx_decoded_phase>minus_delta));
        Rx_decoded_symbols(decoded)=i;
    end;
    Rx_serial_symbols = reshape(Rx_decoded_symbols',1,size(Rx_decoded_symbols,1)*size(Rx_decoded_symbols,2));
    for i = bits_per_symbol:-1:1
        if i~=1
            Rx_binary_matrix(i,:)=rem(Rx_serial_symbols,2);%4进制变为2进制
            Rx_serial_symbols= floor(Rx_serial_symbols/2);
        else
            Rx_binary_matrix(i,:)=Rx_serial_symbols;
        end;
    end;
        baseband_in = reshape(Rx_binary_matrix,1,size(Rx_binary_matrix,1)*size(Rx_binary_matrix,2));
        %误码率计算
        bit_errors = find(baseband_in~=baseband_out);%解调后的二进制序列与发射端的二进制序列进行比较
        bit_error_count(nEN)=bit_error_count(nEN)+size(bit_errors,2);
        total_bits = size(baseband_out,2);
        frame_simulation=frame; 
  end;%for frame=1;Num_repeatsimulation_frame
  SNR_simulation=SNR(nEN) 
  %------------------计算当前信噪比下的信道估计的均方误差
  MSE_Hk(nEN)=0;%估计的频域信道系数,Hk 和理想的频域信道系数,Hk_theory
  for i_f=1:Num_repeatsimulation_frame
        MSE_Hk(nEN)=MSE_Hk(nEN)+sum(abs(H_fbuffer_estimate(i_f,:)-H_fbuffer_theory(i_f,:)).^2);
  end;
  MSE_Hk1(nEN)=1/Num_repeatsimulation_frame/carrier_count*MSE_Hk(nEN);   
    %------------------      
end;% for nEN=1;Num_simulationSNR    
    bit_error_rate = bit_error_count/total_bits/Num_repeatsimulation_frame;
    fprintf(fp,'IFFT_bin_length=%5d,TotalLength_training=%3d,symbols_per_carrier=%5d,Rayleigh_channel_model=%1d,Num_repeatsimulation_frame=%4d.\n ',...
        IFFT_bin_length,TotalLength_training,symbols_per_carrier,Rayleigh_channel_model,Num_repeatsimulation_frame);
    fprintf('IFFT_bin_length=%5d,TotalLength_training=%3d,symbols_per_carrier=%5d,Rayleigh_channel_model=%1d,Num_repeatsimulation_frame=%4d.\n ',...
        IFFT_bin_length,TotalLength_training,symbols_per_carrier,Rayleigh_channel_model,Num_repeatsimulation_frame);
    fprintf(fp,'Lenght_CP=%3d,TotalLength_training=%4d,Path_number=%3d,,symbols_per_carrier=%4d,\n ',...
        Lenght_CP,TotalLength_training,Path_number,symbols_per_carrier);
    fprintf('Lenght_CP=%3d,TotalLength_training=%4d,Path_number=%3d,,symbols_per_carrier=%4d,\n ',...
        Lenght_CP,TotalLength_training,Path_number,symbols_per_carrier);
    fprintf('SNR=%3.1f dB.\n ',SNR);
    fprintf(fp,'BER=%f8.7\n',bit_error_rate);fprintf('BER=%f8.7\n',bit_error_rate);
    fclose(fp);
%------------------
semilogy(SNR,bit_error_rate);
grid;%plot 
xlabel('SNR   (dB) ');ylabel(' BER'); 
figure;
semilogy(SNR,MSE_Hk1);
grid;%plot 
xlabel('SNR   (dB) ');ylabel(' MSE '); 
%------------------        
fprintf('-----------------ofdmchannelest4 routine is over !确定性多径信道系数，仿真结束----------------- ');