基于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

算法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 ');
?編輯：黃飛