作者：Duoqiang Liu，来源： FPGA算法工程师微信公众号
我们知道，对于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算法
算法2：QR分解
下列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 = [-5:5:20];
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(1:N_ps:N_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 ');
