ofdm.m

本代码是基于OFDM系统的一种转换域估计算法

%   ----------------------------------
%   A:  ++++     Definitions    ++++
%   ----------------------------------

% 1. set ofdm system parameters

% Band_system = 20MHZ;
% T_max_multipath = 16us ; mutilpath time delay of outdoor mobile communication is ~us level 
% Tg >= T_max_multipath 这里取 Tg = 16us;
% T_OFDM = 8Tg = 144us; Tb = 128us; Ts = Tb/N = 128/128 = 1us;
% Ts is sample time(抽样率)
% Tg=(1/4;1/8;1/16;1/32)Tb ; here Tg is 1/8
% Ng = Tg/Ts = 16; N_OFDM = N + Ng = 128 + 16 = 144; %采样点数
% fd = v*fc/vc ; fd is doppler frequency formula, vc = 3e8, fc is carrier frequency 
% in the 802.16e，there is 192 data subcarriers,8 pilot subcarriers，55 safeguard subcarriers(27 in high frequency，
% 28 in low frequency)，1 DC subcarrier。
%在用matlab研究OFDM时，往往会用到ifft（）和fft（）这两个函数，但是由于matlab软件本身的算法原因，在使用IFFT和FFT进行
%OFDM调制与解调时，一定要注意能量的归一化。具体来说，在使用IFFT后，要乘以sqrt(N)，而在FFT后则要除以sqrt(N)。

number_ofdm_signal=50;       %Number of ofdm symbol for one loop 250
SNR_dB=[0 4 8 12 16 20]; % SNR per bit, i.e. Eb/N0, in dB

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));
 DFT_err_ber=zeros(1,length(SNR_dB));
DCT_err_ber=zeros(1,length(SNR_dB));

lmmse_mse=zeros(1,length(SNR_dB));
ls_mse=zeros(1,length(SNR_dB));
lr_lmmse_mse=zeros(1,length(SNR_dB));
DFT_mse=zeros(1,length(SNR_dB));
DCT_mse=zeros(1,length(SNR_dB));

% 2. MODULATION MODE

BPSK = 1;   %modulation type is bpsk
QPSK = 2;   %modulation type is qpsk
QAM16 = 3;  %modulation type is 16qam
mod_type = QAM16; 

% 4. PILOT MODE
BLOCK=1;
COMB=2;
pilot_type=BLOCK;
pilot_insert_flag=0;
%pilot_type=COMB;

nloop=1;  %Number of simulation loops

%   ---------------------------------------------
%   B:  %   +++++   TRANSMITTER     +++++
%   ---------------------------------------------
for i=1:length(SNR_dB) %每个SNR点上仿真若干次
    ls_error_bit=0;
    lmmse_error_bit=0;
  lr_lmmse_error_bit=0;  
  DFT_error_bit=0;
 DCT_error_bit=0;
    total_bit_num=0;
    
  for iii=1:nloop
    
    % 1. chose the modulation type
     [seridata,para,fftlen]=bitcreate(mod_type,number_ofdm_signal);     
    
     [dSource,pilot_symbol_bit,symbol_bit] = baseMapping(seridata,mod_type);  %按照不同的方式进行调制
     
    % 2. convert serial data to parall data    
    
     paradata=reshape(dSource,fftlen,number_ofdm_signal);  

    % 3. insert pilot 
          
     [insert_pilot_out,pilot_num,pilot_sequence,pilot_space]=insert_pilot(pilot_symbol_bit,paradata,pilot_type,mod_type);
  
    % 4. IFFT transform
    
     ifft_out=ifft(insert_pilot_out);
        % Nomalize the power
    %ifft_out=ifft_out.*sqrt(fftlen); 

    % 5. Guard intervals were inserted to eliminate ISI caused by multipath fading    
     
     guard_internal=ceil(fftlen*1/8);   %Length of guard internal
     ofdm_guard_interval_out=guard_interval_insert(ifft_out,fftlen,guard_internal,number_ofdm_signal);

    % 6. convert  paradata into seridata  
   
     ofdm_guard_interval_out_serial=reshape(ofdm_guard_interval_out,1,(fftlen+guard_internal)*(number_ofdm_signal+pilot_num));

    %   ---------------------------------------------
    %   C:  %   +++++   CHANNEL     +++++
    %   ---------------------------------------------  

    % 1. multipath Rayleigh fading and doppler spread channel

        %假设功率延迟谱服从负指数分布~exp(-t/trms),trms=(1/4)*cp时长；
        %t在0~cp时长上均匀分布
        %若cp时长为16e-6s，可以取5径延迟如下
        multipath_number=5;
        %delay=[0 2e-6 4e-6 8e-6 12e-6];
        delay=[0 1e-6 2e-6 3e-6 4e-6];
        %trms=4e-6;   %trms为多径平均延迟
        trms=4e-6; %trms为多径平均延迟
        var_pow=10*log10(exp(-delay/trms)); %var_pow是各条多径的功率系数
        doppler_frequency=5;%最大doppler频率为132Hz       
        t_interval=1e-6;%采样间隔为1us
        counter=300000;%各径信道的采样点间隔，应该大于信道采样点数。由以上条件现在信道采样点数
        count_begin=(iii-1)*(5*counter);%每次仿真信道采样的开始位置
        trms_1=trms/t_interval;
        t_max=4e-6/t_interval;
        %信道采样点数，每个调制符号采一个点
        passchan_ofdm_symbol=multipath_chann(ofdm_guard_interval_out_serial,multipath_number,var_pow,delay,doppler_frequency,t_interval,counter,count_begin);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Contaminating the transmitted data with AWGN
    % spow=sum(ich3.^qch3)./number_ofdm_signal./para;
    % attn=0.5*spow*sr/br*10.^(-EbN0/10);
    % attn=sqrt(attn) + eps;
    % [ich4,qch4]=comb(ich3,qch3,attn);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % 2. add AWGN to multipath-signal
     
     snr=10^(SNR_dB(i)/10);   
     code_power=norm(passchan_ofdm_symbol)^2/(length(passchan_ofdm_symbol)); %信号的符号功率 =var(passchan_ofdm_symbol)
     
    
  % [receive_ofdm_symbol,awgn]=add_noise_1(code_power,passchan_ofdm_symbol,symbol_bit,snr); % use  wgn
    sgma=sqrt(code_power/(2*snr));%sgma如何计算，与当前SNR和信号平均能量有关系
  % [receive_ofdm_symbol,awgn]=add_noise(sgma,passchan_ofdm_symbol); % use  normrnd
    [receive_ofdm_symbol,awgn]=addnoise(passchan_ofdm_symbol,sgma); % use randn
    %   ---------------------------------------------
    %   D:  %   +++++   RECEIVER     +++++
    %   ---------------------------------------------  

    % 1. The guard interval was removed from received signals ich4 and qch4
    guard_interval_remove_out=guard_interval_remove(receive_ofdm_symbol,(fftlen+guard_internal),guard_internal,(number_ofdm_signal+pilot_num));
    
    % 2. FFT transform
    
    fft_out=fft(guard_interval_remove_out);
    % nomalize the power
    % fft_out=fft_out./sqrt(fftlen);

    % 3. channel estimation and signal detection

    % LS and zero-force
    ls_zf_detect_sig=ls_estimation(fft_out,pilot_space,pilot_sequence,pilot_num);
    % LMMSE and zero-force
    lmmse_zf_detect_sig=lmmse_estimation(fft_out,pilot_space,pilot_sequence,pilot_num,trms_1,t_max,snr);
    %low_rank_lmmse and zero-force
    [low_rank_lmmse_sig,dlamdaa,dlamda_sor]=lr_lmmse_estimation(fft_out,pilot_space,pilot_sequence,pilot_num,trms_1,t_max,snr,guard_internal);
    % transform domain: DFT
  [DFT_sig,L]=DFT_estimation(fft_out,pilot_space,pilot_sequence,pilot_num);
    % transform domain: DCT
[DCT_sig,LL]=DCT_estimation(fft_out,pilot_space,pilot_sequence,pilot_num); 
    
    % 4. choose the demodulation type

    ls_receive_bit_sig=deMapping(ls_zf_detect_sig,mod_type,fftlen*number_ofdm_signal);
    lmmse_receive_bit_sig=deMapping(lmmse_zf_detect_sig,mod_type,fftlen*number_ofdm_signal);
  lr_lmmse_receive_bit_sig=deMapping(low_rank_lmmse_sig,mod_type,fftlen*number_ofdm_signal);
 DFT_receive_sig=deMapping(DFT_sig,mod_type,fftlen*number_ofdm_signal);
 DCT_receive_sig=deMapping(DCT_sig,mod_type,fftlen*number_ofdm_signal);
    % 5. count error bit number and ber
   
    ls_err_num=error_count(seridata,ls_receive_bit_sig);
    lmmse_err_num=error_count(seridata,lmmse_receive_bit_sig);
    lr_lmmse_err_num=error_count(seridata,lr_lmmse_receive_bit_sig);
 DFT_num=error_count(seridata,DFT_receive_sig);
DCT_num=error_count(seridata,DCT_receive_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;
 DFT_error_bit=DFT_error_bit+DFT_num;
 DCT_error_bit=DCT_error_bit+DCT_num;
    
    % 6 count packet error number and per
 end
    
    % 6. get BER

    total_bit_num=total_bit_num+para*number_ofdm_signal;
    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;
 DFT_err_ber(i)=DFT_error_bit/total_bit_num;
DCT_err_ber(i)=DCT_error_bit/total_bit_num;
    
   % 7. get mse
    ls_mse(i)=mse_ls(snr);
    lmmse_mse(i)=mse_lmmse( dlamdaa,snr );
 lr_lmmse_mse(i)=mse_lr_lmmse_mse(dlamda_sor,snr);
  DFT_mse(i)=mse_DFT(fftlen,L,snr);
DCT_mse(i)=mse_DCT(fftlen,LL,snr);
    
end

figure(1);
semilogy(SNR_dB,ls_mse,'y-*');
hold on
semilogy(SNR_dB,lmmse_mse,'r-o');
semilogy(SNR_dB,lr_lmmse_mse,'g-+');
semilogy(SNR_dB,DFT_mse,'b-s');
semilogy(SNR_dB,DCT_mse,'k-d');
grid on
hold off


figure(2);

semilogy(SNR_dB,ls_err_ber,'y-*');
hold on
semilogy(SNR_dB,lmmse_err_ber,'r-o');
semilogy(SNR_dB,lr_lmmse_err_ber,'g-+');
semilogy(SNR_dB,DFT_err_ber,'b-s');
semilogy(SNR_dB,DCT_err_ber,'k-d');
hold off

% %instantaneous number of errors and data bits
% noe2=sum(abs(seridata-demodata1));
% nod2=length(seridata);
% 
% %cumulative number of errors and data bits in noe and nod
% noe=noe+noe2;
% nod=nod+nod2;
% 
% %calculating PER
% if noe2~=0
%      eop=eop+1;
% else
%      eop=eop;

% eop;
% nop=nop+1;
% 
% fprintf('%d\t%e\t%d\n',iii,noe2/nod2,eop);
% %fprintf :built in function
% end % end of nloop
% %obtain the BER and PER by using the following operation
% ber=noe/nod;
% per=eop/nop;
% 
% fprintf('%f\t%e\t%e\t%d\t\n',ebn0,ber,per,nloop);
% fid=fopen('BERofdm.dat','a');
% fprintf(fid,'%f\t%e\t%e\t%d\t\n',ebn0,ber,per,nloop);
% fclose(fid);
% %********************* end of file ****************************************