基于OMP算法的MIMO-OFDM信道估计

我们知道，对于OFDM系统，只要不发生载波间扰(ICI)，即能够保持子波之间的正交性，就能将每一个子载波看做独立的信道。

这种正交性使得接收信号的每个子载波分量可以被表示成发射信号与子载波的信道频率响应的乘积。因此，仅通过估计每个子载波的信道响应就可以恢复发射信号。

总的来说，可以使用发射机和接收机都已知的导频 (Pilot)符号进行信道估计，并且可以利用不同的插值技术来估计导频之间的子载波上的信道响应。

选择 OFDM 系统的信道估计技术时，必须考虑许多系统实现方面的问题，包括性能需求、计算复杂度和信道时变特性。

常用的信道估计方法有LS估计、LS-DFT估计、MMSE估计、OMP估计等。

在多输入多输出正交频分复用(MIMO-OFDM)系统中，相干检测和均衡需要接收端的信道状态信息(CSI)。然而，在真实的无线环境中，CSI是未知的。因此，信道估计在MIMO-OFDM系统中至关重要。

为简便起见，一般将导频辅助MIMO-OFDM信道估计分解为多个SISO-OFDM信道的同时估计。为了获得更好的信道估计性能，当第i个天线发送导频符号时，所有其他天线必须保持静默。

其中，

是导频位置的集合。

本文考虑了一个具有两个发射天线和两个接收天线的MIMO系统，并利用带有QR分解的OMP算法对MIMO-OFDM信道进行了估计。

OMP算法

设为信道系数的估计值。为非零信道系数的指示集合，为残差，并且，为是测量矩阵的列，是迭代次数。步骤如下。

算法1：OMP算法

Step1：初始化

Step2：通过选择与残差

Step3: 将新选择的索引与索引做并集

Step4:对做QR分解

Step5:计算新的残差

Step6:如果,回到Step2。

Step7: 用反代法求解，得到的k个非零系数。OMP以测量矩阵θ和测量向量y作为输入。在第一次迭代中，选择测量矩阵中与测量向量相关性最大的一列，对其进行QR分解。残差被更新并用于下一次迭代。从第二次迭代开始，选择与残差相关性最大的测量矩阵的列。最后经过K次迭代得到θ的正确列集。在θ上进行QR分解，得到最终的Q和R，然后使用Q, R和y作为反代入的输入得到。QR分解使用改进的Gram-Schmidt正交化进行。

算法2：QR分解

Step1:设

Step2:对，进行循环

Step3:对，进行循环

end for

end for

算法3：回代法

令

Step1:设,且

Step2:对，进行循环

MIMO-OFDM信道估计以2 × 2 MIMO-OFDM信道估计为例，该信道估计分解为4个SISO-OFDM信道,,,和，同时进行估计。

下列MIMO-OFDM系统的仿真参数为:

发射机数量(NT) = 2

接收器数量(NR) = 2

FFT点数量(N_fft) = 512;

每个发射天线导频数(N_P) = 128;

循环前缀长度(N_g) = 64;

信道抽头数(L) = 64;

非零信道系数数(S) = 16;

使用的调制方式:16-QAM

编写并执行下面的MATLAB程来估计2×2 MIMO-OFDM信道。2×2 MIMO-OFDM信道估计的Eb/N0与归一化均方误差(NMSE)如图所示。

从图中可以看出，2 × 2 MIMO-OFDM信道估计，OMP算法性能上优于LS估计。当然若从硬件实现角度看，则OMP算法复杂度更高，耗时更长。

注意，在实际通信系统中，需要根据实际应用场景的信道环境，设计所需的波形，选择合适的信道估计方法。

部分示例代码：

 

 

% Program to compute BER performance of OFDM in sparse Rayleigh fading channel


close all; clear ;clc;




N_fft=512; N_g=64;%N_fft/8;
N_ofdm=N_fft+N_g;
N_sym=100;N_ps=8;
N_p=N_fft/N_ps; N_d=N_fft-N_p; % Pilot spacing, Numbers of pilots and data per OFDM symbol
N_bps=4; M=2^N_bps; % Number of bits per (modulated) symbol
Es=1; A=sqrt(3/2/(M-1)*Es); % Signal energy& QAM normalization factor
EbN0s = [-520];


for i=1:length(EbN0s)
    EbN0 = EbN0s(i); 
    % randn('seed',1)
    rng('default');
    NMSE_OMPi=0;
    NMSE_LSi=0;
    for nsym=1:N_sym
        X_p1= 2*(randn(1,N_p)>0)-1; % Pilot sequence generation
        msg_int=randi(1,N_fft-N_p,M); % bit generation
          Data = qammod(msg_int,M)*A; %QAM Modulated symbols
        ip1 = 0; pilot_loc1 = [];
            for k=1:N_fft
                if mod(k,N_ps)==1
                    X1(k) = X_p1(floor(k/N_ps)+1); 
                    pilot_loc1 = [pilot_loc1 k]; 


                    ip1 = ip1+1;
                else
                    X1(k) = Data(k-ip1);
                end
            end


            X_p2= 2*(randn(1,N_p)>0)-1; % Pilot sequence generation
            ip2 = 0;
            pilot_loc2 = [];
            for k=1:N_fft
                if mod(k,N_ps)==1
                    X2(k) = X_p2(floor(k/N_ps)+1);
                    pilot_loc2 = [pilot_loc2 k]; 
                    ip2 = ip2+1;
                else
                    X2(k) = Data(k-ip2);
                end
            end
        x1 = ifft(X1,N_fft); % IFFT
        xt1 = [x1(N_fft-N_g+1:N_fft) x1]; % Add CP
        x2 = ifft(X2,N_fft); % IFFT
        xt2 = [x2(N_fft-N_g+1:N_fft) x2]; % Add CP
        L = 64; %Total number of channel taps
        K =7; % non-zero channel taps
        T = randperm(L);T = T(1:K);
        h11 = zeros(L,1);
        h11(T) = randn(K,1) + 1i*randn(K,1);
        H11= fft(h11',N_fft); 
        channel_length = length(h11); 
        y_channel11 = conv(xt1, h11'); % Channel path (convolution) 
        h12 = zeros(L,1);
        T = randperm(L);T = T(1:K);
        h12(T) = randn(K,1) + 1i*randn(K,1);
        H12= fft(h12',N_fft); channel_length = length(h12); y_channel12= conv(xt1, h12'); % Channel path (convolution)
        h21 = zeros(L,1);T = randperm(L);T = T(1:K);h21(T) = randn(K,1) + 1i*randn(K,1);
        H21= fft(h21',N_fft); channel_length = length(h12); y_channel21 = conv(xt2, h21'); % Channel path (convolution)
        h22 = zeros(L,1);T = randperm(L);T = T(1:K);
        h22(T) = randn(K,1) + 1i*randn(K,1);
        H22= fft(h22',N_fft); channel_length = length(h22);


        y_channel22 = conv(xt2, h22'); % Channel path (convolution)
        yt11= awgn(y_channel11,EbN0,'measured'); %awgn noise added
        y11 = yt11(N_g+1:N_ofdm); % Remove CP
        Y11 = fft(y11); % FFT
        yt12= awgn(y_channel12,EbN0,'measured'); %awgn noise added 
        y12 = yt12(N_g+1:N_ofdm); % Remove CP
        Y12 = fft(y12); % FF
        yt21= awgn(y_channel21,EbN0,'measured'); %awgn noise added
        y21 = yt21(N_g+1:N_ofdm); % Remove CP
        Y21 = fft(y21); % FFT
        yt22= awgn(y_channel22,EbN0,'measured'); %awgn noise added
        y22 = yt22(N_g+1:N_ofdm); % Remove CP
        Y22 = fft(y22); % FFT


        k=1:N_p; Est_LSH11(k) = Y11(pilot_loc1(k))./X_p1(k);
        k=1:N_p; Est_LSH12(k) = Y12(pilot_loc1(k))./X_p1(k);
        k=1:N_p; Est_LSH21(k) = Y21(pilot_loc2(k))./X_p2(k);
        k=1:N_p; Est_LSH22(k) = Y22(pilot_loc2(k))./X_p2(k);
        Xp1=diag(X1(pilot_loc1));Xp2=diag(X2(pilot_loc2));
        yps11=Y11(pilot_loc1); yps12=Y12(pilot_loc1);
        yps21=Y21(pilot_loc2);
        yps22=Y22(pilot_loc2);


            for ii=1:N_p
                ypss11(ii,1)=yps11(ii); ypss12(ii,1)=yps12(ii); 
                ypss21(ii,1)=yps21(ii); ypss22(ii,1)=yps22(ii); 
            end
            xq=1:N_fft;
        % Est_HLS11 = interpolate(Est_LSH11,pilot_loc1,N_fft,'linear');
        Est_HLS11 = interp1(pilot_loc1,Est_LSH11,xq,'linear');
        h_estLS11 = ifft(Est_HLS11,L); 
        % Est_HLS12 = interpolate(Est_LSH12,pilot_loc1,N_fft,'linear');
        Est_HLS12 = interp1(pilot_loc1,Est_LSH12,xq,'linear');
        h_estLS12 = ifft(Est_HLS12,L); 
        % Est_HLS21 = interpolate(Est_LSH21,pilot_loc2,N_fft,'linear');
        Est_HLS21 = interp1(pilot_loc2,Est_LSH21,xq,'linear');
        h_estLS21 = ifft(Est_HLS21,L); 
        % Est_HLS22 = interpolate(Est_LSH22,pilot_loc2,N_fft,'linear');
        Est_HLS22 = interp1(pilot_loc2,Est_LSH22,xq,'linear');
        h_estLS22 = ifft(Est_HLS22,L); 


        W =exp(-1j*2*pi/N_fft*[0:N_fft-1]'*[0:N_fft-1]);
        WW=W(1N_fft,1:L); %sub dft matrix
        AA1=diag(X_p1)*WW; 
        AA2=diag(X_p2)*WW;


        homp11=ompgs(ypss11,AA1,K); homp12=ompgs(ypss12,AA1,K);
        homp21=ompgs(ypss21,AA2,K); homp22=ompgs(ypss22,AA2,K);
        H_org=[h11 h12 h21 h22];
        H_omp=[ homp11; homp12; homp21; homp22];
        H_LS=[h_estLS11;h_estLS12;h_estLS21;h_estLS22];
        NMSE_OMPi= NMSE_OMPi+ (norm(H_org-H_omp')/norm(H_org));
        NMSE_LSi= NMSE_LSi+ (norm(H_org-H_LS')/norm(H_org));
    end


    NMSE_OMP(i)=NMSE_OMPi;
    NMSE_LS(i)=NMSE_LSi;
end
NMSE_LS=NMSE_LS/(N_fft*N_sym);
NMSE_OMP=NMSE_OMP/(N_fft*N_sym);
figure, semilogy(EbN0s',NMSE_LS,'-x', EbN0s',NMSE_OMP,'-s')
xlabel ('Eb/N0 (dB)') ;ylabel ( 'NMSE');legend('LS','OMP');
title(' BER performance of OFDM ');
 编辑：黄飞