信道估计算法matlab代码汇总

最新推荐文章于 2024-02-28 21:35:35 发布

成电卷王

最新推荐文章于 2024-02-28 21:35:35 发布

阅读量3.5k

点赞数 5

文章标签：算法 matlab

本文链接：https://blog.csdn.net/weixin_45199852/article/details/125381531

版权

OFDM信道估计主函数：

% 16QAM modulation, Block-type pilot channel estimation, 
% Doppler frequency shift=40Hz
% BER VS SNR
clc; 
clear;
close all;
pilot_inter = 5;
cp_length = 16;
rankp = [5,16,20];        % rank SVD-LMMSE 
method='QAM16';
method2='Block';
SNR_dB=0:5:40;
fd=40;
BER1 = zeros(1,length(SNR_dB));
BER2 = zeros(1,length(SNR_dB));
BER3 = zeros(1,length(SNR_dB));
BER32 = zeros(1,length(SNR_dB));
BER33 = zeros(1,length(SNR_dB));
for i=1:length(SNR_dB)
    [BER1(i),BER2(i),BER3(i)]=OFDM(pilot_inter,cp_length,SNR_dB(i),fd,rankp(1),method,method2);
end

for i=1:length(SNR_dB)
    [BER1(i),BER2(i),BER32(i)]=OFDM(pilot_inter,cp_length,SNR_dB(i),fd,rankp(2),method,method2);
end

for i=1:length(SNR_dB)
    [BER1(i),BER2(i),BER33(i)]=OFDM(pilot_inter,cp_length,SNR_dB(i),fd,rankp(3),method,method2);
end

figure
% h=axes;
% set(h,'fontsize',25,'fontweight','b')
semilogy(SNR_dB,BER1,'b-o','markerfacecolor','b')
hold on
semilogy(SNR_dB,BER2,'r-s','markerfacecolor','r')
semilogy(SNR_dB,BER3,'g-d','markerfacecolor','g')
semilogy(SNR_dB,BER32,'g->','markerfacecolor','g')
semilogy(SNR_dB,BER33,'g-<','markerfacecolor','g')
grid on
legend('Block-LS','Block-LMMSE',['Block-LMMSE(SVD) rank p=',num2str(rankp(1))],...
    ['Block-LMMSE(SVD) rank p=',num2str(rankp(2))],['Block-LMMSE(SVD) rank p=',num2str(rankp(3))]);
title([method,' ',method2])
xlabel('SNR(dB)')
ylabel('BER(dB)')

%%
%16QAM modulation, Comb-type pilot channel estimation 
%Doppler frequency shift=80Hz
%BER VS SNR
clear;clc;
pilot_inter=5;
cp_length=16;
rankp=[];
method='QAM16';
method2='Comb';
SNR_dB=0:5:40;
fd=80;

for i=1:length(SNR_dB)
[BER1(i),BER2(i),BER3(i)]=OFDM(pilot_inter,cp_length,SNR_dB(i),fd,rankp,method,method2);
end
figure
% h=axes;
% set(h,'fontsize',25,'fontweight','b')
semilogy(SNR_dB,BER1,'r-v','markerfacecolor','r')
hold on
semilogy(SNR_dB,BER2,'b-^','markerfacecolor','b')
semilogy(SNR_dB,BER3,'g-<','markerfacecolor','g')
grid on
h=legend('Comb-LsLinear ','Comb-LsSecondOrder','Comb-LsSpline');
title([method,' ',method2])
xlabel('SNR(dB)')
ylabel('BER(dB)')

%%
%16QAM modulation, block-type pilot channel estimation 
%SNR=40dB
%BER VS Doppler frequency shift
clc;clear
pilot_inter=5;
cp_length=16;
rankp=5; % rank SVD LMMSE 
method='QAM16';
method2='Block';
SNR_dB=40;
fd=20:10:140;

for i=1:length(fd)
    [BER1(i),BER2(i),BER3(i)]=OFDM(pilot_inter,cp_length,SNR_dB,fd(i),rankp,method,method2);
end
% for i=1:length(fd)
%     [BER1(i),BER2(i),BER32(i)]=OFDM(pilot_inter,cp_length,SNR_dB,fd(i),rankp(2),method,method2);
% end
% for i=1:length(fd)
%     [BER1(i),BER2(i),BER33(i)]=OFDM(pilot_inter,cp_length,SNR_dB,fd(i),rankp(3),method,method2);
% end


figure
% h=axes;
% set(h,'fontsize',25,'fontweight','b')
semilogy(fd,BER1,'r-o','markerfacecolor','r')
hold on
semilogy(fd,BER2,'y-d','markerfacecolor','y')
semilogy(fd,BER3,'g-v','markerfacecolor','g')
% semilogy(fd,BER32,'g->','markerfacecolor','g')
% semilogy(fd,BER33,'g-<','markerfacecolor','g')
grid on
legend('Block-LS','Block-LMMSE',['Block-LMMSE(SVD) rank p=',num2str(rankp(1))])
% set(h,'fontsize',20,'fontweight','b')
title([method,' ',method2])
xlabel('Doppler frequency shift(Hz)')
ylabel('BER(dB)')
%%
%QAM16 modulation, comb-type pilot channel estimation 
%SNR=40dB
%BER VS Doppler frequency shift

clear;clc;
pilot_inter=5;
cp_length=16;
rankp=[];
method='QAM16';
method2='Comb';
SNR_dB=40;
fd=20:10:140;

for i=1:length(fd)
[BER1(i),BER2(i),BER3(i)]=OFDM(pilot_inter,cp_length,SNR_dB,fd(i),rankp,method,method2);
end

figure
% h=axes;
% set(h,'fontsize',25,'fontweight','b')
semilogy(fd,BER1,'b->','markerfacecolor','b')
hold on
semilogy(fd,BER2,'r-<','markerfacecolor','r')
semilogy(fd,BER3,'g-v','markerfacecolor','g')
grid on
h=legend( 'Comb-LsLinear(16QAM) ','Comb-LsSecondOrder(16QAM)','Comb-LsSpline(16QAM)');
% set(h,'fontsize',16,'fontweight','b')
title([method, ' ',method2])
xlabel('Doppler frequency shift(Hz)')
ylabel('BER(dB)')

OFDM子函数

function [BER1,BER2,BER3]=OFDM(pilot_inter,cp_length,SNR_dB,fd,rankp,method,method2)

% simulation parameters:
% channel bandwith: 1MHZ
% frequency of carrier fc=1.9GHZ
% number of sub-carrier: 128
% number of OFDM symbols, symbols/carrier: 100
% interval of sub-carrier 7.8125KHZ
% length of cyclic prefix: cp_length
% pilot symbol: pilot_symbol
% number of pilot interval: pilot_inter
% number of multipath: 5
% maximum Doppler frequency shift: fd
% sub-channel time delay: delay=[0 2e-6 4e-6 8e-6 12e-6]
% method: modulation model QAM16 or QPSK
% method2: Block-type pilot or Comb-type pilot

ofdm_symbol_num=100;               % symbol/carrier, also is number of OFDM symbols

switch method
    case 'QAM16'
        pilot_symbol_bit=[0 0 0 1];   % pilot symbol
        bit_per_carrier=4;            % bit/symbol, for 16QAM is 4bits, while for QPSK is 2bits
    case 'QPSK'
        pilot_symbol_bit=[0 1];       % pilot symbol
        bit_per_carrier=2;            % bit/symbol, for 16QAM is 4bits, while for QPSK is 2bits
end

% for i=1:length(SNR_dB)
    
    bit_source = input_b(128,ofdm_symbol_num,bit_per_carrier); % obtain the original binary signal waiting for transmit
    [nbit,mbit]=size(bit_source);
    total_bit_num=nbit*mbit;   % the total bit transmitted 
    switch method
        case 'QAM16'
            map_out = map_16qam(bit_source);         % QAM16 modulation
            switch method2
                case 'Block'
                    [insert_pilot_out,pilot_num,pilot_sequence]=insert_pilot(pilot_inter,pilot_symbol_bit,map_out,'QAM16');%Block-type pilot
                case 'Comb'
                    [insert_pilot_out,pilot_num,pilot_sequence]=insert_pilot(pilot_inter,pilot_symbol_bit,map_out','QAM16');%Comb-type pilot
                    insert_pilot_out=insert_pilot_out';
            end
        case 'QPSK'
            % QPSK modulation：(1 1)->pi/4;(0 1)->3*pi/4;(0 0)->-3*pi/4;(1,0)->-pi/4;
            X=reshape(bit_source,1,[]);
            s=(X.*2-1)/sqrt(2);
            sreal=s(1:2:length(X));
            simage=s(2:2:length(X));
            X1=sreal+1j.*simage;
            map_out=reshape(X1,128,[]); % QPSK modulation
            clear X sreal simage X1
            switch method2
                case 'Block'
                    [insert_pilot_out,pilot_num,pilot_sequence]=insert_pilot(pilot_inter,pilot_symbol_bit,map_out,'QPSK');%Block-type pilot
                case 'Comb'
                    [insert_pilot_out,pilot_num,pilot_sequence]=insert_pilot(pilot_inter,pilot_symbol_bit,map_out','QPSK');%Comb-type pilot
                    insert_pilot_out=insert_pilot_out';
            end
    end
    
    ofdm_modulation_out=ifft(insert_pilot_out);            % IFFT
    ofdm_cp_out=insert_cp(ofdm_modulation_out,cp_length);  % insert CP
    
    %======================fading rayleigh channel===============
    %suppose the power spectrum follows a negative exponential
    %distribution: exp(-t/trms), where trms=(1/4)*cp
    %t is uniformly distributed in [0,cp]
    
    num=5;                   % number of sub-channels
    delay=[0 2e-6 4e-6 8e-6 12e-6];
    trms = 0.25e-6*cp_length;
    var_pow = 10*log10(exp(-delay/trms));
    t_interval=1e-6;         % sampling interval 1e-6
    counter=200000;          % number of samples
    %     count_begin=(l-1)*(5*counter);%location for begaining the sampling
    count_begin=0;
    trms_1=trms/t_interval;
    t_max=12e-6/t_interval;
    % pass the channel
    passchan_ofdm_symbol = multipath_chann(ofdm_cp_out,num,var_pow,delay,fd,t_interval,counter,count_begin);
    % add noise
    snr=10^(SNR_dB/10);
    [nnl,mml]=size(passchan_ofdm_symbol);
    spow=0;
    for k=1:nnl
        for b=1:mml
            spow=spow+real(passchan_ofdm_symbol(k,b))^2+imag(passchan_ofdm_symbol(k,b))^2;
        end
    end
    spow1 = spow/(nnl*mml);
    sgma = sqrt(spow1/(2*snr));
    receive_ofdm_symbol = add_noise(sgma,passchan_ofdm_symbol);           % recieved OFDM signal
    
    cutcp_ofdm_symbol=cut_cp(receive_ofdm_symbol,cp_length);  % remove CP
    
    ofdm_demodulation_out=fft(cutcp_ofdm_symbol);             % FFT
    
    switch method2
        
        case 'Comb'
            
            switch method
                case 'QPSK'
                    bit_source=reshape(s,bit_per_carrier*128,[]);
            end
            
            ls_zf_detect_sig_Linear=ls_estimation(ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,method2,'Linear');%LS with linear interpolation
            ls_zf_detect_sig_SecondOrder=ls_estimation(ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,method2,'SecondOrder');%LS with second-Order interpolation
            ls_zf_detect_sig_Spline=ls_estimation(ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,method2,'Spline');%LS with Cubic spline  interpolation

            ls_receive_bit_sig_Linear=de_map(ls_zf_detect_sig_Linear,method);
            ls_receive_bit_sig_SecondOrder=de_map(ls_zf_detect_sig_SecondOrder,method);
            ls_receive_bit_sig_Spline=de_map(ls_zf_detect_sig_Spline,method);
                        
            ls_err_Linear=error_count(bit_source,ls_receive_bit_sig_Linear)/total_bit_num;
            ls_err_SecondOrder=error_count(bit_source,ls_receive_bit_sig_SecondOrder)/total_bit_num;
            ls_err_Spline=error_count(bit_source,ls_receive_bit_sig_Spline)/total_bit_num;
           
        case 'Block'
            
            switch method
                case 'QPSK'
                    bit_source=reshape(s,bit_per_carrier*128,[]);
            end
            
            ls_zf_detect_sig_None=ls_estimation(ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,method2,'None');%LS
            lmmse_zf_detect_sig=lmmse_estimation(map_out,ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,trms_1,t_max,snr);%LMMS
            low_rank_lmmse_sig=lr_lmmse_estimation(map_out,ofdm_demodulation_out,pilot_inter,pilot_sequence,pilot_num,trms_1,t_max,snr,rankp);%采用低秩LMMSE估计算法及迫零检测得到的接收信号
            
            ls_receive_bit_sig_None=de_map(ls_zf_detect_sig_None,method);
            lmmse_receive_bit_sig=de_map(lmmse_zf_detect_sig,method);
            lr_lmmse_receive_bit_sig=de_map(low_rank_lmmse_sig,method);
            
            ls_err_None=error_count(bit_source,ls_receive_bit_sig_None)/total_bit_num;
            lmmse_err=error_count(bit_source,lmmse_receive_bit_sig)/total_bit_num;
            lr_lmmse_err=error_count(bit_source,lr_lmmse_receive_bit_sig)/total_bit_num;
            
    end
    
% end

switch method2
    
    case 'Comb'
        
        BER1=ls_err_Linear;
        BER2=ls_err_SecondOrder;
        BER3=ls_err_Spline;
%         semilogy(SNR_dB,ls_err_Linear,'r-o')
%         hold on
%         semilogy(SNR_dB,ls_err_SecondOrder,'g-d')
%         semilogy(SNR_dB,ls_err_Spline,'b-d')
%         
%         grid on
%         legend('LsLinear ','LsSecondOrder','LsSpline')
%         title([method,' ',method2])
%         xlabel('SNR(dB)')
%         ylabel('BER(dB)')
%         
    case 'Block'
                
        BER1=ls_err_None;
        BER2=lmmse_err;
        BER3=lr_lmmse_err;

%         
%         semilogy(SNR_dB,ls_err_None,'b-*')
%         hold on
%         semilogy(SNR_dB,lmmse_err,'r-o')
%         semilogy(SNR_dB,lr_lmmse_err,'g-p')
%         grid on
%         legend('LS','LMMSE','LMMSE(SVD)')
%         title([method,' ',method2])
%         xlabel('SNR(dB)')
%         ylabel('BER(dB)')
end

LMMSE函数

function output=lmmse_estimation(bit_source,input,pilot_inter,pilot_sequence,pilot_num,trms,t_max,snr)

temp=reshape(bit_source,1,[]);
beta=mean(abs(temp).^2)*mean((1./abs(temp)).^2);
% beta=17/9;
[N,NL]=size(input);
Rhh=zeros(N,N);
for k=1:N
    for l=1:N
        Rhh(k,l)=(1-exp((-1)*t_max*((1/trms)+1j*2*pi*(k-l)/N)))./(trms*(1-exp((-1)*t_max/trms))*((1/trms)+1j*2*pi*(k-l)/N));
    end
end
output=zeros(N,NL-pilot_num);
i=1;
count=0;
while i<=NL
    Hi=input(:,i)./pilot_sequence;
    Hlmmse=Rhh*(Rhh+(beta/snr)\eye(N))*Hi;
    count=count+1;
    if  count*pilot_inter<=(NL-pilot_num)
        for p=((count-1)*pilot_inter+1):count*pilot_inter
            output(:,p)=input(:,(i+p-(count-1)*pilot_inter))./Hlmmse;
        end
    else
        for p=((count-1)*pilot_inter+1):(NL-pilot_num)
            output(:,p)=input(:,(i+p-(count-1)*pilot_inter))./Hlmmse;
        end
    end
    i=i+pilot_inter+1;
end
end

LS代码

function output=ls_estimation(input,pilot_inter,pilot_sequence,pilot_num,method2,interpolation)

switch method2
    
    case 'Block'
        
        switch interpolation
            case 'None'
                [N,NL]=size(input);
                output=zeros(N,NL-pilot_num);
                i=1;
                count=0;
                while i<=NL
                    Hi=input(:,i)./pilot_sequence;
                    count=count+1;
                    if  count*pilot_inter<=(NL-pilot_num)
                        for j=((count-1)*pilot_inter+1):count*pilot_inter
                            output(:,j)=input(:,(i+j-(count-1)*pilot_inter))./Hi;
                        end
                    else
                        for j=((count-1)*pilot_inter+1):(NL-pilot_num)
                            output(:,j)=input(:,(i+j-(count-1)*pilot_inter))./Hi;
                        end
                    end
                    i=i+pilot_inter+1;
                end
        end
        
    case 'Comb'
        [N,NL]=size(input);
        output=zeros(size(input));
        pilot_local=1:pilot_inter+1:(pilot_num-1)*(pilot_inter+1)+1;  %location of pilots
        H=zeros(size(input)); 
        H(pilot_local,:)=input(pilot_local,:)./repmat(pilot_sequence',pilot_num,1); %response in the location of pilots
        
        switch interpolation
            
            case 'Linear'
                
                adjacent=zeros(size(input,1),2);
                idx=removerows((1:N)',pilot_local);
                for i=1:length(idx)
                    [~,temp1]=sort(abs(idx(i)-pilot_local));
                    adjacent(idx(i),:)=pilot_local(temp1(1:2));
                end
                for i=1:length(idx)
                    H(idx(i),:)=H(adjacent(idx(i),1),:)+(H(adjacent(idx(i),2),:)-H(adjacent(idx(i),1),:))*(idx(i)-adjacent(idx(i)))/pilot_inter;
                end
                
                output=input./H;
                output=removerows(output, pilot_local);
                
            case 'SecondOrder'
                
                alfa=1/N;
                c0=alfa*(alfa-1)/2;
                c1=-(alfa-1)*(alfa+1);
                c2=alfa*(alfa+1)/2;
                
                adjacent=zeros(size(input,1),3);
                idx=removerows((1:N)',pilot_local);
                for i=1:length(idx)
                    [~,temp1]=sort(abs(idx(i)-pilot_local));
                    adjacent(idx(i),:)=pilot_local(temp1(1:3));
                end
                for i=1:length(idx)
                    H(idx(i),:)=c1*H(adjacent(idx(i),1),:)+c0*H(adjacent(idx(i),2),:)+c2*H(adjacent(idx(i),3),:);
                end
                
                output=input./H;
                output=removerows(output, pilot_local);
                
            case 'Spline'
                
                idx=removerows((1:N)',pilot_local);
                for j=1:NL
                    H(idx,j)=spline(pilot_local,H(pilot_local,j),idx);
                end
                output=input./H;
                output=removerows(output, pilot_local);
                
        end
end

SVD-LMMSE

function output=lr_lmmse_estimation(bit_source,input,pilot_inter,pilot_sequence,pilot_num,trms,t_max,snr,cp)

temp=reshape(bit_source,1,[]);
beta=mean(abs(temp).^2)*mean((1./abs(temp)).^2);

% beta=17/9;
[N,NL]=size(input);
Rhh=zeros(N,N);
for k=1:N
    for l=1:N
        Rhh(k,l)=(1-exp((-1)*t_max*((1/trms)+1j*2*pi*(k-l)/N)))./(trms*(1-exp((-1)*t_max/trms))*((1/trms)+1j*2*pi*(k-l)/N));
    end
end

[U,D]=eig(Rhh);
dlamda=diag(D);
[dlamda_sort,IN]=sort(dlamda);
lamda_sort = zeros(1,N);
IN_new = zeros(1,N);
for k=1:N
    lamda_sort(k)=dlamda_sort(N-k+1);
    IN_new(k)=IN(N-k+1);
end
U_new=zeros(N,N);
for k=1:N
    U_new(:,k)=U(:,IN_new(k));
end

delta=zeros(N,1);
for k=1:N
    if k<=cp
       delta(k)=lamda_sort(k)/(lamda_sort(k)+beta/snr);
    else
       delta(k)=0;
    end
end
D_new=diag(delta);
output=zeros(N,NL-pilot_num);
i=1;
count=0;
while i<=NL
    Hi=input(:,i)./pilot_sequence;
    Hlr=U_new*D_new*(U_new')*Hi;
    count=count+1;
    if  count*pilot_inter<=(NL-pilot_num)
        for p=((count-1)*pilot_inter+1):count*pilot_inter
            output(:,p)=input(:,(i+p-(count-1)*pilot_inter))./Hlr;
        end
    else
        for p=((count-1)*pilot_inter+1):(NL-pilot_num)
            output(:,p)=input(:,(i+p-(count-1)*pilot_inter))./Hlr;
        end
    end
    i=i+pilot_inter+1;
end
end