channel_est_ls.m

LS盲信道估计的代码程序 matble环境下运行

% subspace method for blind channel estimation
%
% Copyright: Xiaohua(Edward) Li, Assistant Professor
%            Department of Electrical and Computer Engineering
%            State University of New York at Binghamton
%            http://ucesp.ws.binghamton.edu/~xli
% June 2003
%

clear all;clc;

SNR=[5:5:40];%信躁比范围
echoF=0;   % turn on/off output display
j=sqrt(-1); %% 由于后面用到j参数，所以这里要做重新定义！
mh=[-0.049+j*0.359 0.482-j*0.569 -0.556+j*0.587 1 -0.171+j*0.061; % channel 阶数为4，共有5个点，因为从0开始计数
    0.443-j*0.0364 1 0.921-j*0.194 0.189-j*0.208 -0.087-j*0.054;
    -0.221-j*0.322 -0.199+j*0.918 1 -0.284-j*0.524 0.136-j*0.19;
    0.417+j*0.030 1 0.873+j*0.145 0.285+j*0.309 -0.049+j*0.161];

[L,M]=size(mh);%求矩阵的行数和列数
M=M-1; N=50; d=M+N; % L: antenna #. M: channel length. N: smoothing. 
h_temp=mh';% 对矩阵共扼转制
h=conj(h_temp(:));%求共扼

for k_SNR=1:length(SNR)    
for kk=1:100
    
dB=SNR(k_SNR);     
s=randsrc(1,N);%产生1行N列的随即数据1或-1

X=zeros(N-M+1,M+1,L);  % received signals 3阶零矩阵
for i=1:L,
  %x=filter(h((i-1)*(M+1)+1:i*(M+1)),1,s);
  x=conv(h((i-1)*(M+1)+1:i*(M+1)),s(1:N-M));%juanji 
  n=randn(size(x))+sqrt(-1)*randn(size(x));
  n=n/norm(n)*10^(-dB/20)*norm(x); %%% 
  x=x+n;
  for j=1:N-M
  X(j,:,i)=x(j:j+M);  %%% 只是完成转置，并无共轭；
  end
end

F=[];   %%%%%% 此处可根据BIFFT程序作相应修改加快运行速度
for i=1:L-1
  temp=[];
  for j=i+1:L
   temp1=[zeros(N-M+1,(i-1)*(M+1)),X(:,:,j),zeros(N-M+1,(j-i-1)*(M+1)),-X(:,:,i),zeros(N-M+1,(L-j)*(M+1))];%行组合
   temp=[temp;temp1] ;%列组合
  end
  F=[F;temp];
end

%%%%%%%%%%%%%% 
[U1,S1,V1]=svd(F'*F);          % solve equation Ah=0 by SVD 3个矩阵的积
hb=U1(:,L*(M+1));           % channel estimation
hhh=[];
for i=1:L
    ttt=hb((i-1)*(M+1)+1:i*(M+1));
    hhh=[hhh;ttt(end:-1:1)];
end
% hb=hhh;
%%%%% Compare channel estimation MSE
hb_h=mean(hb./h);             
hb1=hb/hb_h;
squ_err_h(k_SNR,kk)=sqrt((h-hb1)'*(h-hb1))/sqrt(h'*h); %%% 
bias(k_SNR,kk)=sum(abs(hb1-h))/(L*(M+1));%另均方差最小
qh=hb'*F'*F*hb;

if echoF,       % plot channels
  figure;
  subplot(211), te=length(h); 
  plot(1:te,real(hb1),'bo-',1:te,real(h),'r+-'), grid
  legend('Estimated','Accurate')
  title('Real part of Channel');
  
  subplot(212),
  plot(1:te,imag(hb1),'bo-',1:te,imag(h),'r+-'); grid, 
  legend('Estimated','Accurate')
  title('Imag Part of Channel'),
  xlabel(['hb/h=' num2str(hb_h)]);
end
kk 
end
k_SNR
end
figure,plot(SNR,10*log10(mean(bias'.^2)))
grid on;