%function ofdmchannelest4()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%對每幀OFDM增加循環(huán)前綴，確定性多徑信道系數(shù)
%幀結構：前面1－4個訓練符號，接著是（symbols_per_carrier＋1）個數(shù)據(jù)符號。
%該主程序OFDM信道估計,是對參考文獻中程序修改后得到的。
%參考文獻呂愛琴1 ,田玉敏1 ,朱明華05OFDM信道估計帶MATLAB程序
%ofdmchannelest4.m 修改自ofdmchannelest3.m
%現(xiàn)版本作者Qingmin Meng tested on May4,2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc
fprintf('確定性多徑信道系數(shù)下,OFDM信道估計仿真,如果仿真幀數(shù)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; %選擇調制階數(shù)，如2 為4PSK或DQPSK
% symbols_per_carrier=1;
symbols_per_carrier=50;
TotalLength_training=1; %訓練序列所形成的符號個數(shù)
% TotalLength_training=4;
%  fprintf('TotalLength_training必須為1 或者4,symbols_per_carrier必須為1 或者50.\n');
Lenght_CP=IFFT_bin_length/16;
Rayleigh_channel_model=1;
%------多徑信道系數(shù)定義
% 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多徑數(shù)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，存放共軛符號
    % 信號發(fā)射
    baseband_out=round(rand(1,baseband_out_length)) ;%產(chǎn)生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)變?yōu)橐粋€復數(shù)
    % 加訓練序列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);%然后級聯(lián)形成更多行的訓練序列矩陣
        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);%訓練序列矩陣與發(fā)射信號矩陣級聯(lián)，訓練序列加在前行
    IFFT_modulation =zeros(TotalLength_training+symbols_per_carrier+1,IFFT_bin_length); %在頻域緩沖區(qū)，置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);%--完成準靜態(tài)信道系數(shù)的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; %時域信道信道沖擊響應或系數(shù)
            for p=1:Path_number
                ss(p,:)=[zeros(1,delay(p)) s zeros(1,delay(Path_number)-delay(p))];                            
            end ; 
            S=CH_Data*ss; %信道沖擊響應與發(fā)射信號的卷積。注意多徑信號相加，實際是采用移位相乘來實現(xiàn)卷積 .           
            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，得到頻域信道系數(shù)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;%分子除以分母，得到頻域信道系數(shù)H的倒數(shù)，即1/Hk
    %--存儲估計
                H_fbuffer_estimate(frame,:)=1./Rx_training_symbols0;%存儲估計的頻域信道系數(shù),Hk             
                H_fbuffer_theory(frame,:)=H_f(carriers);           %存儲理想的頻域信道系數(shù),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，需要發(fā)現(xiàn)接收信號處于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進制變?yōu)?進制
            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);%解調后的二進制序列與發(fā)射端的二進制序列進行比較
        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;%估計的頻域信道系數(shù),Hk 和理想的頻域信道系數(shù),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 !確定性多徑信道系數(shù)，仿真結束----------------- ');