dbbcp_fangzhen.m

OFDM频偏估计优化算法代码

%
%
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
%profile on

%参数初始化
%
%
ofdm_syb_num = 500;         %要仿真的符号总数
N = 128;                       %每个符号中的点数
cp_ratio = 1/4;              %cp所占符号的长度的百分比
FS = N;                         %fft/ifft点数
 H = [1,0.81,0.68,0.55,0.45,0.37,0.3,0.25,0.2,0.17,...
     0.34,0.18,0.1,0.06,0.04,0.02,0.01,0.002]' ;     %cost 259 信道抽头滤波器的系数向量，列向量,18 multi-channels

%H = [1,0.5,0.25];

SNR = 60;                     %信噪比向量(dB)
average_num = 480;          % 估计协方差矩阵的平均窗口长度
LEVEL = 0.99;              % 西电的多径信道数目估计算法中的门限值

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%生成初始QPSK信号矩阵
%output:    init_signal_matrix   N by ofdm_syb_num
%
init_signal_matrix = zeros(N,ofdm_syb_num);%初始QPSK信号矩阵规模

for signal_num_loop=1:ofdm_syb_num;
    %rand('state',signal_num_loop*SNR);%每个符号的种子数不同
    init_signal_matrix(:,signal_num_loop) = -1+2*round(rand(N,1))+...
        i*(-1+2*round(rand(N,1)));%将产生的QPSK信号添加到矩阵中，行数代表ofdm符号中的点数，列数代表ofdm符号个数
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %生成ifft后，加上了cp头的信号矩阵
 %output:   cp_signal_ifft   N+N*cp_ratio by ofdm_syb_num
 %
 signal_ifft =FS.* ifft(init_signal_matrix , FS);%
 cp_lenght = N*cp_ratio;%
 cp_block_matrix = signal_ifft((N-cp_lenght+1):N , :);%生成cp块矩阵
 cp_signal_ifft = [cp_block_matrix ; signal_ifft];%带cp头的ofdm发送符号
 
 %%%%%%%%%%%%%%%        completed OFDM symbel genetation (include CP)       %%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %信号通过多径高斯信道
 %output:   channal_awgn_signal   (N+N*cp_ratio)*ofdm_syb_num  by  1
 %
 cp_signal_ifft_serial = cp_signal_ifft(:);%信号矩阵变列向量    N*cp_ratio+N  by  1
 
 channal_signal = filter(H,1,cp_signal_ifft_serial);%过多径信道    N*cp_ratio+N   by  1
 
 channal_awgn_signal = awgn(channal_signal,SNR,'measured',SNR);%加高斯白噪
 %Y = AWGN(X,SNR,SIGPOWER) when SIGPOWER is numeric, it represents 
 %   the signal power in dBW. When SIGPOWER is 'measured', AWGN measures
 %   the signal power before adding noise.
 %Y = AWGN(X,SNR,SIGPOWER,STATE) resets the state of RANDN to STATE.
 
 %%%%===================observe  zone==========================%%%%
%   lenght_serial = length(cp_signal_ifft_serial)
%   lenght_channal_signal = length(channal_signal)
%  t = 1:1:lenght_channal_signal;
%  subplot(2,1,1)
%  h = plot(t,real(channal_signal),'r',t,real(channal_awgn_signal),'--b');
%  axis([1,100,-inf,inf])
%  xlabel('sample number')
%  ylabel('real part')
%  legend(h,'real(channal_signal)','real(channal_awgn_signal)')
%  
%  subplot(2,1,2)
%  h = plot(t,imag(channal_signal),'r',t,imag(channal_awgn_signal),'--b');
%  axis([1,100,-inf,inf])
%  xlabel('sample number')
%  ylabel('imaginary part')
%  legend(h,'imag(channal_signal)','imag(channal_awgn_signal)')
 
 
 %%%%================================================%%%%%%%
 
 
 %
 %
 %
 lenght_H = length(H);%多径数目
 lenght_cpsyb = N+N*cp_ratio;%带cp的总符号长度
 
 receive_signal_matrix = reshape(channal_awgn_signal,(lenght_cpsyb),ofdm_syb_num);%将接收信号转换成矩阵形式
 
 observe_receive_signal_matrix = receive_signal_matrix...
     (lenght_H:lenght_cpsyb , :);%取SB分解法的观察向量矩阵
 
  [M1,N1] = size(observe_receive_signal_matrix);%求接收信号观测矩阵的规模
 observe_receive_signal_matrix = [zeros(M1,1),observe_receive_signal_matrix];%扩充接收信号观测矩阵的规模
                                                                                                                %第一列补0，方便求估计协方差矩阵的迭代初始值
 
     
 estimit_covariance_matrix = zeros(M1,M1);%生成估计协方差矩阵的规模
 
 for loop = 1:average_num;
 estimit_covariance_matrix =  estimit_covariance_matrix + ...
     observe_receive_signal_matrix(:,loop)*( observe_receive_signal_matrix(:,loop))' ;%平均窗口长度减1的观察矩阵的和
 end
 
 estimit_covariance_matrix = (1/average_num)*estimit_covariance_matrix;%生成迭代算法时估计协方差矩阵的初始化矩阵
 %U_matrix = zeros(M1,M1,(ofdm_syb_num - average_num+1));%生成svd分解的特征值向量矩阵，径表示估计矩阵的个数，每个估计矩阵对应一个二维矩阵
 
 estimit_multi_channel_num = zeros(2,(ofdm_syb_num - average_num+1));%生成多径数目估计矩阵，行表示两种多径数目测度方法，列表示每个协方差估计矩阵对应的估计数目
 
 %ofdm symbels number is N,and average windows length is K,so the number of
 %estimition is N-K+1
 for covariance_loop = (average_num+1) : 1 : (N1+1);;%average_num==K     ,     N1==N
     
     estimit_covariance_matrix = estimit_covariance_matrix + 1/average_num*...
         (    observe_receive_signal_matrix(:,covariance_loop)*( observe_receive_signal_matrix(:,covariance_loop) )' -...
        observe_receive_signal_matrix(:,(covariance_loop-average_num))*( observe_receive_signal_matrix(:,(covariance_loop-average_num ) ) )'     );
                                                                                                                                                                  %迭代算法生成估计协方差矩阵
    
    [U,S,V] = svd(estimit_covariance_matrix);%对相关矩阵进行奇异值分解 
    %U_matrix(:,:,(covariance_loop - average_num)) = U;%
   s = diag(S);%将分解的奇异值对角阵转化为列向量

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%  西电的多径信道数目估计算法 %%%%%%%%%%%%%%%%%%%%%%%%
   
   ss = s.^2;%
   sum_ss = sum(ss);%列向量s的平方和
   lenght_s = length(s);%列向量s的长度
   sum_eigenvalue = 0;%特征值累和量初始值

   for sum_loop = 1:lenght_s;
       sum_eigenvalue = sum_eigenvalue + ss(sum_loop);%对s的平方向量逐值累加
         if( sqrt( sum_eigenvalue/sum_ss ) >= LEVEL )%归一化测度是否大于门限
             estimit_multi_channel_num(1,(covariance_loop - average_num)) = sum_loop;%估计的多径数目保存在估计多径数目矩阵中的相应位置
             break;%若大于门限则保存多径估计值，跳出循环
         end
   end
  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%  MDL多径信道数目估计算法 %%%%%%%%%%%%%%%%%%%%%%%%%
   s1 = s(1:60);
   lenght_s1 = length(s1);
   
   T_sph = zeros(1,(lenght_s1 -1));%生成中间参数向量规模
   F_MDL = zeros(1,(lenght_s1 -1));%生成测度向量规模
   
   for MDL_loop = 1:1:(lenght_s1 -1);
       T_sph(1,MDL_loop) = ( sum( s( (MDL_loop+1):lenght_s1 ) ) / (lenght_s1 - MDL_loop) ) /...
           (  prod( s( (MDL_loop+1):lenght_s1 ) )^( 1/(lenght_s1 - MDL_loop) )  );%中间参数
       
       F_MDL(1,MDL_loop) = average_num*(lenght_s1 - MDL_loop)*log10(T_sph(1,MDL_loop)) +...
           ( MDL_loop*(2*lenght_s1 - MDL_loop)*log10(average_num) )/2;%测度
   end
   [minimum,index_min]=min(F_MDL);%找最小值及其索引
    estimit_multi_channel_num(2,(covariance_loop - average_num)) = index_min;%MDL方法估计出的多径数目
   loop=covariance_loop
 end   
   
   xidian_mean = mean(estimit_multi_channel_num(1,:),2)
   MDL_mean = mean(estimit_multi_channel_num(2,:),2)
   xidian_standard = std(estimit_multi_channel_num(1,:),0,2)
   MDL_standard = std(estimit_multi_channel_num(1,:),0,2)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   %%%%===================observe  zone==========================%%%%
%    figure(1)
%    stem3(estimit_multi_channel_num);
%    title('estimite of mutil-channel numbers with xidian/MDL method');
%    xlabel('number of ofdm symbels');
%    ylabel('different method:xidian and MDL');
%    zlabel('estimite of mutil-channel numbers');
%    
   %figure(2)
   
   
   %%%%================================================%%%%%%%
   
   %profile off
%profile viewer
%
%save dbbfangzhen.mat