DFT信道估计步骤及实例

clc;clear all; close all; figure(1), clf, figure(2), clfNfft=1024;Ng=Nfft/8;Nofdm=Nfft+Ng;Nsym=100;Nps=4; Np=Nfft/Nps; Nd=Nfft-Np; % 导频间隔、每个 OFDM 符号的导频数和数据Nbps=4; M=2^Nbps; % 每个（调制）符号的位数Es=5; A=sqrt(3

ddatalent

6681人浏览 · 2022-02-08 14:45:57

ddatalent · 2022-02-08 14:45:57 发布

在OFDM系统中，基于DFT的信道估计算法也经常被使用。DFT信道估计主要有以下几个步骤:
(1)在实载波导频处，利用已知导频符号，使用LS信道估计算法计算得到实载波导频处信道频域响应：（原理详见此处）
$\hat{H}_{P\_AC}(k)=\frac{Y_{AC}(mk)}{X_{AC}(mk)}，k=0,1,...,N_{P\_AC}-1$
其中， $Y_{AC}(mk)$ 为实载波导频处接收频域信号， $X_{AC}(mk)$ 为实载波导频处发送频域信号， $A C$ 是实载波的意思， $P\_AC$ 是实导频载波。
(2) 通过 IDFT 变换获得实际导频处的信道冲击响应：
$\hat{h}_{P\_AC}(n)=IDFT{\{\hat{H}_{P\_AC}(k)\}}，n=0,1,...,N_{P\_AC}-1$
(3) 对 CP 长度后的冲击响应进行置零，对信道冲击响应进行降噪：
$h_{P\_AC}(n)= \begin{cases} \hat{h}_{P\_AC}(n) & n=0,1,...,N_{CP}-1 \\ 0 & n=N_{CP},N_{CP}+1,...,N_{P\_AC}-1 \end{cases}$
(4) 在 $h_{P\_AC}(n)$ 中进行插零得到 $h_{AC}(n)$ ，插值后长度为实载波数量 $N_{AC}$ ，再通过DFT变换，得到实载波处的频域响应 $H_{AC}(k)$
$h_{P\_AC}(n)\overset{zero\_padding}{\implies}h_{AC}(n)，n=0,1,...,N_{AC}-1 \\ H_{AC}(k)=DFT\{h_{AC}(n)\}，k=0,1,...,N_{AC}-1$
在这里插入图片描述
经过上述步骤，即可得到DFT优化结果。
以下为实际操作例子

在这里插入图片描述

clc;clear all; close all; figure(1), clf, figure(2), clf
Nfft=1024;  Ng=Nfft/8;  Nofdm=Nfft+Ng;  Nsym=100;
Nps=4; Np=Nfft/Nps; Nd=Nfft-Np; % 导频间隔、每个 OFDM 符号的导频数和数据
Nbps=4; M=2^Nbps; % 每个（调制）符号的位数
Es=5; A=sqrt(3/2/(M-1)*Es); % 信号能量和 QAM 归一化因子
SNRs = 30;  
for i=1:length(SNRs)
   SNR = SNRs(i); 
   rand('seed',1); randn('seed',1);
   MSE = zeros(1,6); nose = 0;
   for nsym=1:Nsym
      Xp = 2*(randn(1,Np)>0)-1;    % Pilot sequence generation
      msgint=randi([0,M-1],1,Nd);    % bit generation
      Data = qammod(msgint,M)*A;
      %Data = modulate(mod_object, msgint); Data = modnorm(Data,'avpow',1)*Data;   % normalization
      ip = 0;    pilot_loc = [];
      for k=1:Nfft
         if mod(k,Nps)==1
             X(k) = Xp(floor(k/Nps)+1); pilot_loc = [pilot_loc k]; ip = ip+1;
         else       
             X(k) = Data(k-ip);
         end
      end
      x = ifft(X,Nfft);                            % IFFT
      xt = [x(Nfft-Ng+1:Nfft) x];                  % Add CP
      h = [(randn+j*randn)*3 (randn+j*randn)*0.7 (randn+j*randn)*0.2 (randn+j*randn)*0.01];     % generates a  channel
      H = fft(h,Nfft); channel_length = length(h); % True channel and its time-domain length
      H_power_dB = 10*log10(abs(H.*conj(H)));      % True channel power in dB
      y_channel = conv(xt, h);                     % Channel path (convolution)

      yt = awgn(y_channel,SNR,'measured');  
      y = yt(Ng+1:Nofdm);                   % Remove CP
      Y = fft(y);                           % FFT
      for m=1:1
         if m==1
             H_est = LS_CE(Y,Xp,pilot_loc,Nfft,Nps,'linear'); method='LS信道估计'; % LS estimation with linear interpolation
         elseif m==2
             H_est = LS_CE(Y,Xp,pilot_loc,Nfft,Nps,'spline'); method='LS-spline'; % LS estimation with spline interpolation
         else  
             H_est = MMSE_CE(Y,Xp,pilot_loc,Nfft,Nps,h,SNR); method='MMSE'; % MMSE estimation
         end
         H_est_power_dB = 10*log10(abs(H_est.*conj(H_est)));
         h_est = ifft(H_est); h_DFT = h_est(1:channel_length); 
         H_DFT = fft(h_DFT,Nfft); % DFT-based channel estimation
         H_DFT_power_dB = 10*log10(abs(H_DFT.*conj(H_DFT)));
         if nsym==1
%            figure(1), subplot(319+2*m), plot(H_power_dB,'b','linewidth',1); grid on; hold on;
%            plot(H_est_power_dB,'r+','Markersize',4,'linewidth',1); 
%            title(method); xlabel('子载波编号'); ylabel('功率[dB]');
%            legend('True Channel',method);  set(gca,'fontsize',9)
%            subplot(320+2*m), plot(H_power_dB,'b','linewidth',1); grid on; hold on;
%            plot(H_DFT_power_dB,'r+','Markersize',4,'linewidth',1); 
%            title([method ' with DFT']); xlabel('子载波编号'); ylabel('Power[dB]');
%            legend('True Channel',[method ' with DFT']);  set(gca,'fontsize',9)
        
            figure(1),subplot(1,2,1)
             plot(H_power_dB,'b','linewidth',1); grid on; hold on;
             plot(H_est_power_dB,'r+','Markersize',4,'linewidth',1); 
             title(method); xlabel('子载波编号'); ylabel('功率[dB]');
             legend('True Channel',method);  set(gca,'fontsize',9)
             subplot(1,2,2)
             plot(H_power_dB,'b','linewidth',1); grid on; hold on;
             plot(H_DFT_power_dB,'r+','Markersize',4,'linewidth',1); 
             title([method ' +  DFT']); xlabel('子载波编号'); ylabel('Power[dB]');
             legend('True Channel',[method ' with DFT']);  set(gca,'fontsize',9)
         end
         MSE(m) = MSE(m) + (H-H_est)*(H-H_est)';
         MSE(m+3) = MSE(m+3) + (H-H_DFT)*(H-H_DFT)';
      end
      Y_eq = Y./H_est;
      Y_eqDFT = Y./H_DFT;
      if nsym>=Nsym-10
        figure(2)
        subplot(221), plot(X,'.','Markersize',5),title('原星座图');axis('equal'), set(gca,'fontsize',9), axis([-5 5 -5 5]),hold on,
        subplot(222), plot(Y,'.','Markersize',5),title('信道补偿前星座图');axis('equal'), set(gca,'fontsize',9), axis([-20 20 -20 20]),hold on,       
        subplot(223), plot(Y_eq,'.','Markersize',5),title('信道补偿后星座图');axis('equal'), set(gca,'fontsize',9),axis([-5 5 -5 5]), hold on,       
        subplot(224), plot(Y_eqDFT,'.','Markersize',5),title('信道补偿后星座图（DFT）');axis('equal'), set(gca,'fontsize',9), axis([-5 5 -5 5]),hold on,       
      end
      ip = 0;
      for k=1:Nfft
         if mod(k,Nps)==1, ip=ip+1;  else  Data_extracted(k-ip)=Y_eq(k);  end
      end
      msg_detected = qamdemod(Data_extracted/A,M);
      nose = nose + sum(msg_detected~=msgint);
   end   
   MSEs(i,:) = MSE/(Nfft*Nsym)
end
Number_of_symbol_errors=nose

% figure(3), clf, semilogy(SNRs',MSEs(:,1),'-x', SNRs',MSEs(:,4),'-o')
% legend('LS-linear','LS_DFT')
% xlabel('SNR(dB)')
% ylable('MSE')
% legend('LS-linear','MMSE')
% fprintf('MSE of LS-linear/LS-spline/MMSE Channel Estimation = %6.4e/%6.4e/%6.4e\n',MSEs(end,1:3));
% fprintf('MSE of LS-linear/LS-spline/MMSE Channel Estimation with DFT = %6.4e/%6.4e/%6.4e\n',MSEs(end,4:6));

function H_LS = LS_CE(Y,Xp,pilot_loc,Nfft,Nps,int_opt)
% LS channel estimation function
% Inputs:
%       Y         = Frequency-domain received signal
%       Xp        = Pilot signal
%       pilot_loc = Pilot location
%       N         = FFT size
%       Nps       = Pilot spacing
%       int_opt   = 'linear' or 'spline'
% output:
%       H_LS      = LS channel etimate

%MIMO-OFDM Wireless Communications with MATLAB㈢   Yong Soo Cho, Jaekwon Kim, Won Young Yang and Chung G. Kang
%?2010 John Wiley & Sons (Asia) Pte Ltd

Np=Nfft/Nps; k=1:Np; LS_est(k) = Y(pilot_loc(k))./Xp(k);  % LS channel estimation
if  lower(int_opt(1))=='l', method='linear'; else method='spline';  end
H_LS = interpolate(LS_est,pilot_loc,Nfft,method); % Linear/Spline interpolation
end

function [H_MMSE]= MMSE_CE(Y,Xp,pilot_loc,Nfft,Nps,h,SNR)% MMSE信道估计
%输入:
% Y 频域接收信号
% Xp 导频信号
% pilot_loc:导频位置
% Nfft  FFT大小
% Nps 导频间隔
% h 信道脉冲响应
% SNR:信噪比[dB]
%输出:
% H_MMSE :MMSE信道估计
snr = 10^(SNR*0.1);
Np=Nfft/Nps; k=1:Np; H_tilde = Y(1,pilot_loc(k))./Xp(k);  % LS estimate
k=0:length(h)-1; %k_ts = k*ts; 
hh = h*h'; tmp = h.*conj(h).*k; %tmp = h.*conj(h).*k_ts;
r = sum(tmp)/hh;    r2 = tmp*k.'/hh; %r2 = tmp*k_ts.'/hh;
tau_rms = sqrt(r2-r^2);     % rms delay
df = 1/Nfft;  %1/(ts*Nfft);
j2pi_tau_df = j*2*pi*tau_rms*df;
K1 = repmat([0:Nfft-1].',1,Np); K2 = repmat([0:Np-1],Nfft,1);
rf = 1./(1+j2pi_tau_df*(K1-K2*Nps));
K3 = repmat([0:Np-1].',1,Np); K4 = repmat([0:Np-1],Np,1);
rf2 = 1./(1+j2pi_tau_df*Nps*(K3-K4));
Rhp = rf;
Rpp = rf2 + eye(length(H_tilde),length(H_tilde))/snr;
H_MMSE = transpose(Rhp*inv(Rpp)*H_tilde.');  % MMSE channel estimate

end

function H_interpolated = interpolate(H_est,pilot_loc,Nfft,method)
% Input:        H_est    = Channel estimate using pilot sequence
%           pilot_loc    = location of pilot sequence
%                Nfft    = FFT size
%              method    = 'linear'/'spline'
% Output: H_interpolated = interpolated channel

%MIMO-OFDM Wireless Communications with MATLAB㈢   Yong Soo Cho, Jaekwon Kim, Won Young Yang and Chung G. Kang
%?2010 John Wiley & Sons (Asia) Pte Ltd

if pilot_loc(1)>1
  slope = (H_est(2)-H_est(1))/(pilot_loc(2)-pilot_loc(1));
  H_est = [H_est(1)-slope*(pilot_loc(1)-1)  H_est]; pilot_loc = [1 pilot_loc];
end
if pilot_loc(end)<Nfft
  slope = (H_est(end)-H_est(end-1))/(pilot_loc(end)-pilot_loc(end-1));  
  H_est = [H_est  H_est(end)+slope*(Nfft-pilot_loc(end))]; pilot_loc = [pilot_loc Nfft];
end
if lower(method(1))=='l', H_interpolated = interp1(pilot_loc,H_est,[1:Nfft]);   
 else      H_interpolated = interp1(pilot_loc,H_est,[1:Nfft],'spline');
end  
end