dbbcp2.m

OFDM频偏估计优化算法代码

%function [ output_args ] = Untitled1( input_args )
%UNTITLED1 Summary of this function goes here
%  Detailed explanation goes here


%协方差矩阵的估计，采用接收数据分块进行估计
%
%


close all
clear all
%参数初始化
%
%
ofdm_syb_num = 100;         %要仿真的符号总数
N = 256;                       %每个符号中的点数
cp_ratio = 1/8;              %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 信道抽头滤波器的系数向量，列向量
SNR = [0:2:70];                     %信噪比向量(dB)
average_num = 100;          % 估计协方差矩阵的平均窗口长度
LEVEL = 0.99;              % 西电的多径信道数目估计算法中的门限值
simu_num_per_snr = 100;%每个信噪比下的仿真次数

estimit_xidian = zeros(2,length(SNR));%西电算法在不同信噪比下的估计参数平均值，第一行代表均值，第二行代表方差
estimit_MDL = zeros(2,length(SNR));%MDL算法在不同信噪比下的估计参数平均值，第一行代表均值，第二行代表方差

%不同信噪比的循环
%
%
for SNR_loop = 1:1:length(SNR)
    
    estimit_multi_channel_num = zeros(2,simu_num_per_snr);%生成多径数目估计矩阵，行表示两种多径数目测度方法，列表示每个协方差估计矩阵对应的估计数目
    
    %每个信噪比下仿真次数的循环,循环体内包含完整的算法及过程
    %
    %
    for simulation_loop = 1:1:simu_num_per_snr

        %生成初始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发送符号

        %信号通过多径高斯信道
        %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(1,SNR_loop));%加高斯白噪

        %信号的子空间分解算法（subspace_based）
        %output : s   协方差矩阵的奇异值列向量
        %
        lenght_H = length(H);%多径数目
        lenght_cpsyb = N+ cp_lenght;%带cp的总符号长度

        receive_signal_matrix = reshape(channal_awgn_signal,lenght_cpsyb,ofdm_syb_num);%将接收信号转换成矩阵形式          test point

        observe_receive_signal_matrix = receive_signal_matrix...
            (lenght_H:lenght_cpsyb , :);%取SB分解法的观察向量矩阵

        [M1,N1] = size(observe_receive_signal_matrix);%求接收信号观测矩阵的规模

        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))' ;%平均窗口长度的观察矩阵的和
        end

        estimit_covariance_matrix = (1/average_num)*estimit_covariance_matrix;%生成迭代算法时估计协方差矩阵的初始化矩阵

        [U,S,V] = svd(estimit_covariance_matrix);%对相关矩阵进行奇异值分解
        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,simulation_loop) = sum_loop;%估计的多径数目保存在估计多径数目矩阵中的相应位置
                break;%若大于门限则保存多径估计值，跳出循环
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%  MDL多径信道数目估计算法 %%%%%%%%%%%%%%%%%%%%%%%%%
        s1 = s(1:average_num);%
        lenght_s1 = length(s1);

        T_sph = zeros(1,(lenght_s1 -1));%生成中间参数向量规模
        F_MDL = zeros(1,(lenght_s1 -1));%生成测度向量规模
        A = log10(average_num);

        for MDL_loop = 1:1:(lenght_s1 -1);
            T_sph(1,MDL_loop) = ( sum( s1( (MDL_loop+1):lenght_s1 ) ) / (lenght_s1 - MDL_loop) ) /...
                (  prod( s1( (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) +1 ) +...
                ( MDL_loop*(2*lenght_s1 - MDL_loop)*A )/2;%测度
        end
        [minimum,index_min]=min(F_MDL);%找最小值及其索引
        estimit_multi_channel_num(2,simulation_loop) = index_min;%MDL方法估计出的多径数目

    end
    
    estimit_xidian(1,SNR_loop) = mean(estimit_multi_channel_num(1,:));%
    estimit_xidian(2,SNR_loop) = std(estimit_multi_channel_num(1,:));%
    
    estimit_MDL(1,SNR_loop) = mean(estimit_multi_channel_num(2,:));%
    estimit_MDL(2,SNR_loop) = std(estimit_multi_channel_num(2,:));%
    
    present_SNR = SNR(1,SNR_loop)
    
end


figure(1)
plot(SNR,estimit_xidian(1,:),'b*-');
xlabel('SNR');
ylabel('estimite numbers');
title('xidian mean');


figure(2)
plot(SNR,estimit_xidian(2,:),'b*-');
xlabel('SNR');
ylabel('standred');
title('xidian standred');


figure(3)
plot(SNR,estimit_MDL(1,:),'b*-');
xlabel('SNR');
ylabel('estimite numbers');
title('MDL mean');


figure(4)
plot(SNR,estimit_MDL(2,:),'b*-');
xlabel('SNR');
ylabel('standred');
title('MDL standred');