信道估计matlab代码_信道估计代码

主函数：

clear;clc;close all;
Nfft = 256;          % 子载波数量
Ng = Nfft/8;         % 循环前缀长度
Nofdm = Nfft + Ng;   % OFDM符号长度
Nsym = 100;
Nps = 4;             % 导频间隔
Np = Nfft/Nps;       % 导频数
Nbps = 4;            % 每符号=4bit
M = 2^Nbps;
A = sqrt(3/2/(M-1)*1);  % QAM归一化因子
SNR = 30;
sq2 = sqrt(2);
MSE = zeros(1,6);
nose = 0;
for nsym = 1:Nsym
    Xp = 2*(randn(1,Np)>0)-1;
    msgint = randi([0,15],Nfft-Np,1);
    Data = A*qammod(msgint,M,'gray');
    ip=0;
    pilot_loc = zeros(8,1);
    X = zeros(Nfft,1);
    kk=1;
    for k=1:Nfft
        if mod(k,Nps)==1
            X(k,:)=Xp(floor(k/Nps)+1);
            pilot_loc(kk)=k;
            kk = kk+1;
            ip=ip+1;
        else
            X(k,:) = Data(k-ip,:);
        end
    end
    x=ifft(X,Nfft);                   % IFFT
    xt = [x(Nfft-Ng+1:Nfft);x];       % 加CP
    h = [(randn+1j*randn),(randn+1j*randn)/2];
    H = fft(h,Nfft);
    channel_length = length(h);
    H_power_dB = 10*log10(abs(H.*conj(H)));
    y_channel = conv(xt,h);
    yt = awgn(y_channel,SNR,'measured');
    y = yt(Ng+1:Nofdm);               % 去CP
    Y = fft(y);
    for m = 1:3
        if m==1
            H_est = LS_CE(Y,Xp,pilot_loc,Nfft,Nps,'linear');
            method = 'LS-linear';
        elseif m==2
            H_est = LS_CE(Y,Xp,pilot_loc,Nfft,Nps,'spline');
            method = 'LS-spline';
        else
            H_est = MMSE_CE(Y,Xp,pilot_loc,Nfft,Nps,h,SNR);
            method = 'MMSE';
        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);
        H_DFT_power_dB = 10*log10(abs(H_DFT.*conj(H_DFT)));
        if nsym==1
            subplot(319+2*m)
            plot(H_power_dB,'b','linewidth',1);
            hold on;
            plot(H_est_power_dB,'r:+','Markersize',4);
            axis([0 33 -10 10]);
            xlabel('Subcarrier Index');
            ylabel('Power(dB)');
            legend('True Channel',method);
            set(gca,'fontsize',9);
            subplot(320+2*m)
            plot(H_power_dB,'b','linewidth',1);
            hold on;
            plot(H_DFT_power_dB,'r:+','Markersize',4);
            axis([0 33 -10 10]);
            xlabel('Subcarrier Index');
            ylabel('Power(dB)');
            legend('True Channel','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;
    if nsym>=Nsym-10
        figure
        subplot(211)
        plot(Y,'.','Markersize',5);
        axis([-1.5 1.5 -1.5 1.5]);
        set(gca,'fontsize',9);
        hold on;
        subplot(212)
        plot(Y_eq,'.','Markersize',5);
        axis([-1.5 1.5 -1.5 1.5]);
        set(gca,'fontsize',9);
    end
    ip=0;
    Data_extracted = zeros(Nfft-Np,1);
    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,'gray');
    nose = nose + sum(msg_detected~=msgint);   
    MSEs = MSE/(Nfft*Nsym);
end

LS算法

% LS信道估计
function [H_LS]=LS_CE(Y,Xp,pilot_loc,Nfft,Nps,int_opt)
% 输入
% Y           频域接收信号
% Nfft        载波数量
% Nps         导频间隔
% Xp          导频信号
% pilot_loc   导频位置
% int_opt     插值方式
 
% 输出
% H_LS        LS信道估计          
 
Np = Nfft/Nps;  % 导频个数
k = 1:Np;
LS_est(k) = Y(pilot_loc(k))'./Xp(k);
H_LS = interpolate(LS_est,pilot_loc',Nfft,int_opt);
 
end

MMSE算法

% LS信道估计
function [H_MMSE]=MMSE_CE(Y,Xp,pilot_loc,Nfft,Nps,h,SNR)
% 输入
% Y           频域接收信号
% Nfft        载波数量
% Nps         导频间隔
% Xp          导频信号
% pilot_loc   导频位置
% h           信道冲击相应
% SNR         信噪比
% 输出
% H_MMSE      LS信道估计          
snr = 10^(SNR*0.1);
Np=Nfft/Nps; 
k = 1:Np;
H_tilde = Y(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 = 1j*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*((Rpp)\H_tilde'));  
end

插值算法

function [H_interpolate] = interpolate(H_est,pilot_loc,Nfft,method)
 
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
 
H_interpolate = interp1(pilot_loc,H_est,(1:Nfft),method);