基于压缩感知的ofdm信道估计

clc;clear; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 生成导频位置 %%%%%%%%%%%%%%%%%%%%%%%%%%%% %随机导频 256个子载波插100个导频 P_num1 = 100; %导频个数 N_fft = 256; %256个子载波 M = 4; %调制阶数 frame = 1000; %重复运行次数 snr_begin = 0; %信噪比起始位置 snr_end = 30; signal1=zeros(256,1); q=randperm(256); q1 = find((q > 0)&(q < 256)); q4 = [q1]; qq1 = q(q4); q1 = qq1; signal1(q1(1:P_num1)) = sign(randn(P_num1,1)); a1=abs(signal1);%产生稀疏信号 b1=diag(a1); c1=[1:1:256];%导频坐标 pilot_index1=c1*b1; pilot_index1(find(pilot_index1==0))=[]; pilot_pos = [pilot_index1];%导频位置 n = 1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 生成导频 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i = snr_begin:1:snr_end error_num_bin = 0; mse0 = 0; for j = 1:frame data = randi([0 M-1],(N_fft - length(pilot_pos)),1);%生产随机信号 symbol_bin = qammod(data,M); data_freq = zeros(256,1); data_freq(pilot_pos,:) = 1;%插入导频 data_pos = find(data_freq == 0);%获取数据位，等0就是要插的数据位置 data_freq(data_pos,:) = symbol_bin; data_freq = ifft(data_freq,256,1)*sqrt(256);%傅里叶变换 data_freq_ad_cp = [data_freq(256-30+1:end,:);data_freq];%加cp data_freq_ad_cp = reshape(data_freq_ad_cp,[],1);%串并 %%%%%%%%%%%%%%%%%%%%%%%%%%%% 过多径信道 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [data_mult,h] = multipath(N_fft,data_freq_ad_cp); data_mult = awgn(data_mult,i,'measured'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 接收端 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% re_data_time = data_mult(30+1:end,:); re_data_freq = fft(re_data_time,256,1)/sqrt(256); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 生成部分傅里叶矩阵 %%%%%%%%%%%%%%%%%%%%%% eye_mtx = eye(256); dft_mtx=dftmtx(256);%（傅里也点数，1） W_p_samp=eye_mtx(pilot_pos,:)*dft_mtx;%（256*256）（导频位置，：）*dftmtx(256) X_omp_samp=eye(length(pilot_pos));%导频数*导频数 T_samp=X_omp_samp*W_p_samp; T = T_samp; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 压缩感知信道估计 %%%%%%%%%%%%%%%%%%%%%%% s = re_data_freq(pilot_pos,1); H1 =fft(h,256,1)/sqrt(256);%H [ re_y ] = CS_OMP( s,T,5);%h re_H=fft(re_y,256,1)/sqrt(256);%压缩感知信道估计重构的信道频率响应 re_data_freq = re_data_freq./re_H;%逆滤波 re_data = re_data_freq(data_pos,:);%取出数据 re_symbol_bin = qamdemod(re_data,M);%解qam映射 [num_bin,~] = biterr(data,re_symbol_bin);%误码个数 error_num_bin = error_num_bin + num_bin; diff_value=abs((re_H) -(H1)); mse = mean((diff_value./abs(H1)).^2);%获得均方差mse mse0 = mse0 + mse; end BER(n) = error_num_bin/(frame*length(data_pos)); MSE(n) = mse0/frame; n = n+1; end SNR = snr_begin:1:snr_end; figure(1);%误码率 semilogy(SNR,BER,'r--o'); grid on hold on legend('OMP'); xlabel('SNR(dB)'); ylabel('BER') figure(2);%均方差 semilogy(SNR,MSE,'r--o'); grid on hold on legend('OMP'); xlabel('SNR(dB)'); ylabel('MSE')