sub3.m

基于子空间的MIMO-OFDM系统的盲信道估计

% 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
%
%1发4收
clear all
close all
clc
echoF=1;   % turn on/off output display
dB=15; T=1000;  % SNR, sample amount
L=4; M=4; N=5; d=M+N; % L: antenna #. M: channel length. N: smoothing. 
j=sqrt(-1);           % d: equalization delay
mh=[-0.049+j*0.359 0.482-j*0.569 -0.556+j*0.587 1 -0.171+j*0.061; % channel
    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];
h=[mh(1,:) mh(2,:) mh(3,:) mh(4,:)].';
%h1=ones(20,1);
s=sign(rand(1,T)-0.5);%+2*sign(rand(1,T)-0.5);  % 16 QAM symbols
s=s+sqrt(-1)*(sign(rand(1,T)-0.5));%+2*sign(rand(1,T)-0.5));

s=ifft(s);

TN=T-N+1; X=zeros(L*N,TN); SNR=[];  v=[]; % received signals
Nm=20;%进行10次仿真
squ_err_h=zeros(1,25);
 bias=zeros(1,25);
for dB=1:25
 
 for i3=1:Nm %进行Nm次仿真
for i1=1:L,
  x=filter(h((i1-1)*(M+1)+1:i1*(M+1)),1,s);%将信号s通过滤波器滤波
  %x=s;
  n=randn(size(x))+sqrt(-1)*randn(size(x));%产生随机噪声
  n=n/norm(n)*10^(-dB/20)*norm(x);
  SNR=[SNR 20*log10(norm(x)/norm(n))];%信噪比
  x=x+n;   v=[v n];%接收信号及接收噪声
  for j1=1:TN,
    X((i1-1)*N+1:i1*N, j1)=x(j1+N-1:-1:j1).';
  end
end 
if echoF SNR=SNR; end
ss=std(s)^2;   sv=std(v)^2;%std标准差，ss即为方差
  
%%%%%%%%%%%%%% subspace method begin
Rx=X*X'/TN;                 % calculate correlation matrix
[U0,S0,V0]=svd(Rx);         % SVD to find null subspace
for i=L*N:-1:2,
   if S0(i-1,i-1)-S0(i,i)>S0(i,i), break; end  
end
i=d+1;                      % check rank of null subspace
%i=rank(S0)+1;
if echoF, d=i-1, else d=i-1; end   % display rank
sigma=0;
for i=i:L*N, sigma=sigma+S0(i,i); end   % remove noise
sigma=sigma/(L*N-d);%噪声子空间的平均特征值

Q=zeros(L*(M+1), L*(M+1));             % Construct matrix A (in Q)
for i2=d+1:L*N,
  Vm=zeros(L*(M+1), M+N);
  for j2=1:(M+1),
    for k=1:L,
      Vm((k-1)*(M+1)+j2, j2:(j2+N-1))=U0((k-1)*N+1:k*N, i2).';
    end
  end
  Q=Q+Vm*Vm';
end

[U1,S1,V1]=svd(Q);          % solve equation Ah=0 by SVD
hb=U1(:,L*(M+1));           % channel estimation

%%%%%%% Compare channel estimation MSE
hb_h=mean(hb./h);     %        
hb1=hb/hb_h;          %消
squ_err_h(1,dB)=squ_err_h(1,dB)+sqrt((h-hb1)'*(h-hb1))/sqrt(h'*h);
bias(1,dB)=bias(1,dB)+sum(abs(hb1-h))/(L*(M+1));
qh=hb'*Q*hb;
 end
squ_err_h(1,dB)=squ_err_h(1,dB)/(Nm*L*(M+1));%均方误差
bias(1,dB)=bias(1,dB)/Nm;%平均偏差
end

figure(1)
semilogy(squ_err_h,'-*');xlabel('SNR');ylabel('MSE');title('The MSE of the Channel Estimation');grid
figure(2)
plot(bias,'-*r');xlabel('SNR');ylabel('bias');title('The bias of the Channel Estimation');grid

if echoF, squ_err_h, bias, qh, end
figure(3)
if echoF,       % plot channels
  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

%%%% plot equalization results
%H=zeros(L*N, M+N);  %% channel matrix
%for j=1:N,
%   for k=1:L,
%      H((k-1)*N+j, j:(j+M))=hb1((k-1)*(M+1)+1:k*(M+1)).';
%   end
%end
%Y=H'*U0(:,1:d)*inv(S0(1:d, 1:d)-sigma*eye(d))*U0(:,1:d)'*X; % zero-forcing equalizer

%gd=H'*U0(:,1:d)*inv(S0(1:d, 1:d)-sigma*eye(d))*U0(:,1:d)';
%gd=gd(round(d/2), :).';
%   fh=zeros(M+N,1);
%   for j=1:L
%       fh=fh+conv(h((j-1)*(M+1)+1:j*(M+1)), gd((j-1)*(N)+1:j*(N)));
%   end
%   ISI=[(fh'*fh-max(abs(fh))^2)/max(abs(fh))^2];
%         dmax=find(max(abs(fh))==abs(fh));
%         fh1=fh.'/fh(dmax); F1=gd.'/fh(dmax);
%         MSE=ss*(fh1*fh1'-1)+sv*(F1*F1');
%   if echoF, abs(fh.')/max(abs(fh)),  ISI_MSE=[ISI MSE],    end

%if echoF,
%  subplot(222), plot(s,'ro'), grid, title('Transmitted Symbols')
%  subplot(224), plot(Y(round(d/2),:),'ro'), grid, 
%  title('Estimated Symbols')
%end