OFDM中的LS算法

%OFDM Channel Estimation Based on Comb Pilot
%IFFT_bin_length: IFFT和FFT的点数
%carrier_count: 子载波个数
%bits_per_symbol: 每符号上的比特数
%symbols_per_carrier: 每桢的OFDM符号数
%Ｘ：欲发送的二进制比特流

clear all;
clc;
IFFT_bin_length=128; %IFFT和FFT的点数
carrier_count=100; %%%子载波个数
bits_per_symbol=2; %%每符号上的比特数
symbols_per_carrier=12; %%每桢的OFDM符号数
LI=7 ; %导频之间的间隔
Np=ceil(carrier_count/LI)+1;%导频数 %加1的原因：使最后一列也是导频

N_number=carrier_count*symbols_per_carrier*bits_per_symbol;
carriers=1:carrier_count+Np;

GI=8; % guard interval length
N_snr=25; % 每比特信噪比
snr=5; %信噪比间隔
% symbols_per_carrier*(IFFT_bin_length+GI)
%------------------------------------------------------------
% vector initialization


X=zeros(1,N_number);
X1=[];
X2=[];
X3=[];
X4=[];
X5=[];
X6=[];
X7=[];
Y1=[];
Y2=[];
Y3=[];
Y4=[];
Y5=[];
Y6=[];
Y7=[];
XX=zeros(1,N_number);
dif_bit=zeros(1,N_number);
dif_bit1=zeros(1,N_number);
dif_bit2=zeros(1,N_number);
dif_bit3=zeros(1,N_number);
X=randint(1,N_number);%产生二进制随即序列（非0即1）
%--------------------------------------------------------
%QPSK调制：(1 1)->pi/4;(0 1)->3*pi/4;(0 0)->-3*pi/4;(1,0)->-pi/4;

s=(X.*2-1)/sqrt(2);
sreal=s(1:2:N_number);
simage=s(2:2:N_number);
X1=sreal+j.*simage;


%---------------------------------------------------------
%产生随机导频信号
%--------------------------------------------------------
train_sym=randint(1,2*symbols_per_carrier);
t=(train_sym.*2-1)/sqrt(2);
treal=t(1:2:2*symbols_per_carrier);
timage=t(2:2:2*symbols_per_carrier);
training_symbols1=treal+j.*timage;
training_symbols2=training_symbols1.';
training_symbols=repmat(training_symbols2,1,Np);

%disp(training_symbols)
pilot=1:LI+1:carrier_count+Np;
if length(pilot)~=Np
pilot=[pilot,carrier_count+Np];
end
%--------------------------------------------------------
%串并转换
X2=reshape(X1,carrier_count,symbols_per_carrier).';
%---------------------------------------------------------
%插入导频
signal=1:carrier_count+Np;
signal(pilot)=[];
X3(:,pilot)=training_symbols;
X3(:,signal)=X2;
%X3=cat(1,training_symbols,X2);
IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);
IFFT_modulation(:,carriers)=X3;
%IFFT_modulation(:,conjugate_carriers)=conj(X3);
X4=ifft(IFFT_modulation,IFFT_bin_length,2);
%X5=X4.';
%加保护间隔（循环前缀）
for k=1:symbols_per_carrier;
for i=1:IFFT_bin_length;
X6(k,i+GI)=X4(k,i);
end
for i=1:GI;
X6(k,i)=X4(k,i+IFFT_bin_length-GI);
end
end
%---------------------------------------------------------
%并串转换
X7=reshape(X6.',1,symbols_per_carrier*(IFFT_bin_length+GI));
% f=10

0; r=pi/6;
% for n=1:length(X7)
% % a(n)=cos((2*pi*n)/8);
% a(n)=exp((i*((2*pi*f*n)/(8*length(X7))+r)));
% %c=X7.*a(n);
% end
% b=real(X7.*a);
%
% m1=xcorr(b(:,1:136));
% m2=xcorr(b(:,137:272));
%a=zeros(1000,271);
%for i=1:1000
% b=xcorr(X7(:,1:136));
% c=xcorr(X7(:,137:272));
% a(i,:)=b;
%end
%---------------------------------------------------------
%信道模型:带多普勒频移的瑞利衰落信道

fd=100; %多普勒频移
r=6; %多径数
a=[0.123 0.3 0.4 0.5 0.7 0.8]; %多径的幅度
d=[2 3 4 5 9 13]; %各径的延迟
T=1; %系统采样周期
th=[90 0 72 144 216 288]*pi./180;%相移
h=ones(1,carrier_count);
hh=[];
for k=1:r
%deta=[zeros(1,d(k)-1),1,zeros(1,carrier_count-d(k))];
h1=a(k)*exp(j*((2*pi*T*fd*d(k)/carrier_count)));
%h1=a(k)*exp(j*((2*pi*T*fd*d(k)/carrier_count)));
hh=[hh,h1];
end
h(d+1)=hh;
%noise=randn(1,length(X7))+j.*randn(1,length(X7));

%--------------------------------------------------------

channel1=zeros(size(X7));
channel1(1+d(1):length(X7))=hh(1)*X7(1:length(X7)-d(1));
channel2=zeros(size(X7));
channel2(1+d(2):length(X7))=hh(2)*X7(1:length(X7)-d(2));
channel3=zeros(size(X7));
channel3(1+d(3):length(X7))=hh(3)*X7(1:length(X7)-d(3));
channel4=zeros(size(X7));
channel4(1+d(4):length(X7))=hh(4)*X7(1:length(X7)-d(4));
channel5=zeros(size(X7));
channel5(1+d(5):length(X7))=hh(5)*X7(1:length(X7)-d(5));
channel6=zeros(size(X7));
channel6(1+d(6):length(X7))=hh(6)*X7(1:length(X7)-d(6));

%---------------------------------------------------------------
Tx_data=X7+channel1+channel2+channel3+channel4;
%---------------------------------------------------------------


%---------------------------------------------------------------
%----------------------------------------------------------------
%加高斯白噪声
Error_ber=[];%误比特率
Error_ber1=[];
Error_ber2=[];%误比特率
Error_ber3=[];
Error_ber11=[];
%Error_ser=[];%误符号率
for snr_db=0:snr:N_snr
%
% code_power=0;
% code_power=[norm(Tx_data)]^2/(length(Tx_data));%信号的符号功率
% %bit_power=var(Tx_data);
% bit_power=code_power/bits_per_symbol;%比特功率
% noise_power=10*log10((bit_power/(10^(snr_db/10))));%噪声功率
% noise=wgn(1,length(Tx_data),noise_power,'complex');%产生GAUSS白噪声信号
vn=randn(1,symbols_per_carrier*(IFFT_bin_length+GI))+sqrt(-1)*randn(1,symbols_per_carrier*(IFFT_bin_length+GI));
vn=vn/norm(vn)*10^(-snr_db/20)*norm(Tx_data);
Y7=Tx_data+vn;

%-------------------------------------------------------
%串并变换
Y6=reshape(Y7,IFFT_bin_length+GI,symbols_per_carrier).';
Y6_1=reshape(Tx_data,IFFT_bin_length+GI,symbols_per_carrier).';
%去保护间隔
for k=1:symbols_per_carrier;
for i=1:IFFT_bin_length;
Y5(k,i)=Y6(k,i+GI);

Y5_1(k,i)=Y6_1(k,i+GI);
end
end
Y4=fft(Y5,IFFT_bin_length,2);
Y3=Y4(:,carriers);
Y4_1=fft(Y5_1,IFFT_bin_length,2);
Y3_1=Y4_1(:,carriers);
%-------------------------------------------------------------
%LS信道估计
H=[];
Y2=Y3(:,signal);
Y2_1=Y3_1(:,signal);
Rx_training_symbols=Y3(:,pilot);
Rx_training_symbols_1=Y3_1(:,pilot);
Rx_training_symbols0=reshape(Rx_training_symbols,symbols_per_carrier*Np,1);
Rx_training_symbols0_1=reshape(Rx_training_symbols_1,symbols_per_carrier*Np,1);
training_symbol0=reshape(training_symbols,1,symbols_per_carrier*Np);
training_symbol0_1=reshape(training_symbols,1,symbols_per_carrier*Np);
training_symbol1=diag(training_symbol0);
training_symbol1_1=diag(training_symbol0_1);
%disp(training_symbols)
training_symbol2=inv(training_symbol1);
Hls=training_symbol2*Rx_training_symbols0;
Hls1=reshape(Hls,symbols_per_carrier,Np);
HLs=[];
HLs1=[];
if ceil(carrier_count/LI)==carrier_count/LI
for k=1:Np-1
HLs2=[];
for t=1:LI
HLs1(:,1)=(Hls1(:,k+1)-Hls1(:,k))*(t-1)./LI+Hls1(:,k);
HLs2=[HLs2,HLs1];
end
HLs=[HLs,HLs2];
end
else
for k=1:Np-2
HLs2=[];
for t=1:LI
HLs1(:,1)=(Hls1(:,k+1)-Hls1(:,k))*(t-1)./LI+Hls1(:,k);
HLs2=[HLs2,HLs1];
end
HLs=[HLs,HLs2];
end
HLs3=[];
for t=1:mod(carrier_count,LI)
HLs1(:,1)=(Hls1(:,Np)-Hls1(:,Np-1))*(t-1)./LI+Hls1(:,Np-1);
HLs3=[HLs3,HLs1];
end;
HLs=[HLs,HLs3];
end
%Hls1=Hls.';
training_symbol2_1=inv(training_symbol1_1);
Hls_1=training_symbol2_1*Rx_training_symbols0_1;
Hls1_1=reshape(Hls_1,symbols_per_carrier,Np);
HLs_1=[];
HLs1_1=[];
if ceil(carrier_count/LI)==carrier_count/LI
for k=1:Np-1
HLs2_1=[];
for t=1:LI
HLs1_1(:,1)=(Hls1_1(:,k+1)-Hls1(:,k))*(t-1)./LI+Hls1(:,k);
HLs2_1=[HLs2_1,HLs1_1];
end
HLs_1=[HLs_1,HLs2_1];
end
else
for k=1:Np-2
HLs2_1=[];
for t=1:LI
HLs1_1(:,1)=(Hls1_1(:,k+1)-Hls1_1(:,k))*(t-1)./LI+Hls1_1(:,k);
HLs2_1=[HLs2_1,HLs1_1];
end
HLs_1=[HLs_1,HLs2_1];
end
HLs3_1=[];
for t=1:mod(carrier_count,LI)
HLs1_1(:,1)=(Hls1_1(:,Np)-Hls1_1(:,Np-1))*(t-1)./LI+Hls1_1(:,Np-1);
HLs3_1=[HLs3_1,HLs1_1];
end;
HLs_1=[HLs_1,HLs3_1];
end
%H=repmat(Hls1,symbols_per_carrier,1);%将导频扩展成symbols_per_carrier*carrier_count矩阵
Y1=Y2./HLs;
H_2=ifft(HLs,100,2);
H_3=zeros(12,100);
H_3(:,1:8)=H_2(:,1:8);
H_4=fft(H_3,100,2);

% H_2_1=ifft(HLs_1,IFFT_bin_length,2);
Y11=Y2./H_4;
%Y11=Y2./


%-------------------------------------------------------------------

%-------------------------------

------------------------------
%并串变换
YY=reshape(Y2.',1,N_number/bits_per_symbol);
YY1=reshape(Y1.',1,N_number/bits_per_symbol);
YY11=reshape(Y11.',1,N_number/bits_per_symbol);

%------------------------------------------------------------
%QPSK解调
y_real=sign(real(YY));
y_image=sign(imag(YY));
y_re=y_real./sqrt(2);
y_im=y_image./sqrt(2);
y_real1=sign(real(YY1));
y_image1=sign(imag(YY1));
y_re1=y_real1./sqrt(2);
y_im1=y_image1./sqrt(2);

y_real11=sign(real(YY11));
y_image11=sign(imag(YY11));
y_re11=y_real11./sqrt(2);
y_im11=y_image11./sqrt(2);

r00=[];
r01=[];
r10=[];
r11=[];

r110=[];
r111=[];

for k=1:length(y_real);
r00=[r00,[y_real(k),y_image(k)]];
end;
for k=1:length(y_real1);
r10=[r10,[y_real1(k),y_image1(k)]];
end;
for k=1:length(y_real11);
r110=[r10,[y_real11(k),y_image11(k)]];
end;




for k=1:length(y_re);
r01=[r01,[y_re(k),y_im(k)]];
end;
for k=1:length(y_re1);
r11=[r11,[y_re1(k),y_im1(k)]];
end;
for k=1:length(y_re11);
r111=[r111,[y_re11(k),y_im1(k)]];
end;


% XX(find(r01>0))=1;
%-------------------------------------------------------------
%计算在不同信噪比下的误比特率并作图

dif_bit=s-r01;
dif_bit1=s-r11;
dif_bit11=s-r111;

ber_snr=0; %纪录误比特数
for k=1:N_number;
if dif_bit(k)~=0;
ber_snr=ber_snr+1;
end
end;
ber_snr1=0; %纪录误比特数
for k=1:N_number;
if dif_bit1(k)~=0;
ber_snr1=ber_snr1+1;
end
end


ber_snr11=0; %纪录误比特数
for k=1:N_number;
if dif_bit11(k)~=0;
ber_snr11=ber_snr11+1;
end
end

Error_ber=[Error_ber,ber_snr];
Error_ber1=[Error_ber1,ber_snr1];
Error_ber11=[Error_ber11,ber_snr11];
%hls=(sum(HLs,1))/12;
%mse=abs(hls-h)/abs(h);
end

BER=zeros(1,length(0:snr:N_snr));
BER1=zeros(1,length(0:snr:N_snr));
BER11=zeros(1,length(0:snr:N_snr));

BER=Error_ber./N_number;
BER1=Error_ber1./N_number;
BER11=Error_ber11./N_number;
%-------------------------------------------------------------
%-------------------------------------------------------------
i=0:snr:N_snr;
% semilogy(i,BER,'-*r');
%
% hold on;
semilogy(i,BER1,'-og');
hold on;
semilogy(i,BER11,'-*b');
hold on;
grid on;
xlabel('SNR（dB）');

ylabel('BER');
% legend('No Channel Estimation','LS Channel Estimation','DFT Channel Estimation');
legend('LS Channel Estimation','DFT Channel Estimation');
title('LS和DFT算法误码率曲线图');
hold off