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

我們知道，對于OFDM系統(tǒng)，只要不發(fā)生載波間擾(ICI)，即能夠保持子波之間的正交性，就能將每一個子載波看做獨立的信道。

這種正交性使得接收信號的每個子載波分量可以被表示成發(fā)射信號與子載波的信道頻率響應的乘積。因此，僅通過估計每個子載波的信道響應就可以恢復發(fā)射信號。

總的來說，可以使用發(fā)射機和接收機都已知的導頻 (Pilot)符號進行信道估計，并且可以利用不同的插值技術(shù)來估計導頻之間的子載波上的信道響應。

選擇 OFDM 系統(tǒng)的信道估計技術(shù)時，必須考慮許多系統(tǒng)實現(xiàn)方面的問題，包括性能需求、計算復雜度和信道時變特性。

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

在多輸入多輸出正交頻分復用(MIMO-OFDM)系統(tǒng)中，相干檢測和均衡需要接收端的信道狀態(tài)信息(CSI)。然而，在真實的無線環(huán)境中，CSI是未知的。因此，信道估計在MIMO-OFDM系統(tǒng)中至關重要。

為簡便起見，一般將導頻輔助MIMO-OFDM信道估計分解為多個SISO-OFDM信道的同時估計。為了獲得更好的信道估計性能，當?shù)趇個天線發(fā)送導頻符號時，所有其他天線必須保持靜默。

其中，

是導頻位置的集合。

本文考慮了一個具有兩個發(fā)射天線和兩個接收天線的MIMO系統(tǒng)，并利用帶有QR分解的OMP算法對MIMO-OFDM信道進行了估計。

OMP算法

設為信道系數(shù)的估計值。為非零信道系數(shù)的指示集合，為殘差，并且，為是測量矩陣的列，是迭代次數(shù)。步驟如下。

算法1：OMP算法

Step1：初始化

Step2：通過選擇與殘差

Step3: 將新選擇的索引與索引做并集

Step4:對做QR分解

Step5:計算新的殘差

Step6:如果,回到Step2。

Step7: 用反代法求解，得到的k個非零系數(shù)。OMP以測量矩陣θ和測量向量y作為輸入。在第一次迭代中，選擇測量矩陣中與測量向量相關性最大的一列，對其進行QR分解。殘差被更新并用于下一次迭代。從第二次迭代開始，選擇與殘差相關性最大的測量矩陣的列。最后經(jīng)過K次迭代得到θ的正確列集。在θ上進行QR分解，得到最終的Q和R，然后使用Q, R和y作為反代入的輸入得到。QR分解使用改進的Gram-Schmidt正交化進行。

算法2：QR分解

Step1:設

Step2:對，進行循環(huán)

Step3:對，進行循環(huán)

end for

end for

算法3：回代法

令

Step1:設,且

Step2:對，進行循環(huán)

MIMO-OFDM信道估計以2 × 2 MIMO-OFDM信道估計為例，該信道估計分解為4個SISO-OFDM信道,,,和，同時進行估計。

下列MIMO-OFDM系統(tǒng)的仿真參數(shù)為:

發(fā)射機數(shù)量(NT) = 2

接收器數(shù)量(NR) = 2

FFT點數(shù)量(N_fft) = 512;

每個發(fā)射天線導頻數(shù)(N_P) = 128;

循環(huán)前綴長度(N_g) = 64;

信道抽頭數(shù)(L) = 64;

非零信道系數(shù)數(shù)(S) = 16;

使用的調(diào)制方式:16-QAM

編寫并執(zhí)行下面的MATLAB程來估計2×2 MIMO-OFDM信道。2×2 MIMO-OFDM信道估計的Eb/N0與歸一化均方誤差(NMSE)如圖所示。

從圖中可以看出，2 × 2 MIMO-OFDM信道估計，OMP算法性能上優(yōu)于LS估計。當然若從硬件實現(xiàn)角度看，則OMP算法復雜度更高，耗時更長。

注意，在實際通信系統(tǒng)中，需要根據(jù)實際應用場景的信道環(huán)境，設計所需的波形，選擇合適的信道估計方法。

部分示例代碼：

?

?

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

?

?

?