clear all
clc
format long
%本次仿真載頻為2GHz，帶寬1MHz，子載波數128個，cp為16
%子載波間隔為7.8125kHz
%一個ofdm符號長度為128us，cp長度為16us
%采用16QAM調制方式
%最大doppler頻率為132Hz
%多徑信道為5徑，功率延遲譜服從負指數分布~exp(-t/trms),trms=(1/4)*cp時長，各徑延遲取為delay=[0 2e-6 4e-6 8e-6 12e-6]
pilot_inter=5;%導頻符號間隔為10,可以調整，看不同導頻間隔下的BER情況，和理論公式比較
pilot_symbol_bit=[0 0 0 1];%導頻為常數，對應星座點1+3*j
cp_length=16;%cp長度為16
SNR_dB=[0 2 4 6 8 10 12 14 16];
ls_err_ber=zeros(1,length(SNR_dB));
lmmse_err_ber=zeros(1,length(SNR_dB));
lr_lmmse_err_ber=zeros(1,length(SNR_dB));
for i=1:length(SNR_dB)%每個SNR點上仿真若干次
    ls_error_bit=0;
    lmmse_error_bit=0;
    lr_lmmse_error_bit=0;
    total_bit_num=0;
loop_num=1; %共仿真10次
for l=1:loop_num
    ofdm_symbol_num=100;%每次仿真產生100個ofdm符號,則每次仿真共有100×128個星座映射符號；16QAM調制下，1個星座映射符號包含4個bit
    
    bit_source=input_b(128,ofdm_symbol_num);%為每次仿真產生100個ofdm符號的比特個數，128為每個ofdm符號的子載波個數
    [nbit,mbit]=size(bit_source);
    total_bit_num=total_bit_num+nbit*mbit;
    
    map_out=map_16qam(bit_source);%對一次仿真符號塊進行16QAM映射
    
    [insert_pilot_out,pilot_num,pilot_sequence]=insert_pilot(pilot_inter,pilot_symbol_bit,map_out);%按塊狀導頻結構，對映射后的結果插入導頻序列
    
    ofdm_modulation_out=ifft(insert_pilot_out,128);%作128點逆FFT運算，完成ofdm調制
    
    ofdm_cp_out=insert_cp(ofdm_modulation_out,cp_length);%插入循環前綴
        
    %********************** 以下過程為ofdm符號通過頻率選擇性多徑信道 *************************
    num=5;
    %假設功率延遲譜服從負指數分布~exp(-t/trms),trms=(1/4)*cp時長；
    %t在0~cp時長上均勻分布
    %若cp時長為16e-6s，可以取5徑延遲如下
    delay=[0 2e-6 4e-6 8e-6 12e-6];
    %a=[0.123 0.3 0.4 0.5 0.6];
    trms=4e-6;
    var_pow=10*log10(exp(-delay/trms));
    fd=132;%最大doppler頻率為132Hz
    t_interval=1e-6;%采樣間隔為1us
    counter=200000;%各徑信道的采樣點間隔，應該大于信道采樣點數。由以上條件現在信道采樣點數
    count_begin=(l-1)*(5*counter);%每次仿真信道采樣的開始位置
    trms_1=trms/t_interval;
    t_max=16e-6/t_interval;
    %信道采樣點數，每個調制符號采一個點
    passchan_ofdm_symbol=multipath_chann(ofdm_cp_out,num,var_pow,delay,fd,t_interval,counter,count_begin);
    
    %********************** 以上過程為ofdm符號通過頻率選擇性多徑信道 *************************
    
     %********************** 以下過程為ofdm符號加高斯白噪聲 *************************
    snr=10^(SNR_dB(i)/10);
    [nnl,mml]=size(passchan_ofdm_symbol);
    spow=0;
    for k=1:nnl
      for b=1:mml
        spow=spow+real(passchan_ofdm_symbol(k,b))^2+imag(passchan_ofdm_symbol(k,b))^2;
      end
    end
    spow1=spow/(nnl*mml);        
    sgma=sqrt(spow1/(2*snr));%sgma如何計算，與當前SNR和信號平均能量有關系
    receive_ofdm_symbol=add_noise(sgma,passchan_ofdm_symbol);%加入隨機高斯白噪聲，receive_ofdm_symbol為最終接收機收到的ofdm符號塊
    
    %********************** 以上過程為ofdm符號加高斯白噪聲 *************************
    cutcp_ofdm_symbol=cut_cp(receive_ofdm_symbol,cp_length);%去除循環前綴
    
    ofdm_demodulation_out=fft(cutcp_ofdm_symbol,128);%作128點FFT運算，完成ofdm解調
    
    %********************** 以下就是對接收ofdm信號進行信道估計和信號檢測的過程************************
    ls_zf_detect_sig=ls_estimation(ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num);%采用LS估計算法及迫零檢測得到的接收信號
    lmmse_zf_detect_sig=lmmse_estimation(ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,trms_1,t_max,snr);%采用LMMSE估計算法及迫零檢測得到的接收信號
    low_rank_lmmse_sig=lr_lmmse_estimation(ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,trms_1,t_max,snr,cp_length);%采用低秩LMMSE估計算法及迫零檢測得到的接收信號
    %********************** 以下就是對接收ofdm信號進行信道估計和信號檢測的過程************************
    
    ls_receive_bit_sig=de_map(ls_zf_detect_sig);%16QAM解映射
    lmmse_receive_bit_sig=de_map(lmmse_zf_detect_sig);
    lr_lmmse_receive_bit_sig=de_map(low_rank_lmmse_sig);
    
    %以下過程統計各種估計算法得到的接收信號中的錯誤比特數
    ls_err_num=error_count(bit_source,ls_receive_bit_sig);
    lmmse_err_num=error_count(bit_source,lmmse_receive_bit_sig);
    lr_lmmse_err_num=error_count(bit_source,lr_lmmse_receive_bit_sig);
    
    ls_error_bit=ls_error_bit+ls_err_num;
    lmmse_error_bit=lmmse_error_bit+lmmse_err_num;
    lr_lmmse_error_bit=lr_lmmse_error_bit+lr_lmmse_err_num;
end
%計算各種估計算法的誤比特率

ls_err_ber(i)=ls_error_bit/total_bit_num;
lmmse_err_ber(i)=lmmse_error_bit/total_bit_num;
lr_lmmse_err_ber(i)=lr_lmmse_error_bit/total_bit_num;

end
plot(SNR_dB,ls_err_ber,'b-*')
hold on

plot(SNR_dB,lmmse_err_ber,'r-o')

plot(SNR_dB,lr_lmmse_err_ber,'g-+')
hold off