blind3sub.m

一基于子空间方法的盲信道估计程序

% 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
%
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,:)].';

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));

TN=T-N+1; X=zeros(L*N,TN); SNR=[];  v=[]; % received signals
for i=1:L,
  x=filter(h((i-1)*(M+1)+1:i*(M+1)),1,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 j=1:TN,
    X((i-1)*N+1:i*N, j)=x(j+N-1:-1:j).';
  end
end 
if echoF SNR=SNR, end
ss=std(s)^2;   sv=std(v)^2;
  
%%%%%%%%%%%%%% 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:1,
   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 i=d+1:L*N,
  Vm=zeros(L*(M+1), M+N);
  for j=1:(M+1),
    for k=1:L,
      Vm((k-1)*(M+1)+j, j:(j+N-1))=U0((k-1)*N+1:k*N, i).';
    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=sqrt((h-hb1)'*(h-hb1))/sqrt(h'*h);
bias=sum(abs(hb1-h))/(L*(M+1));
qh=hb'*Q*hb;

if echoF, squ_err_h, bias, qh, end

if echoF,       % plot channels
  subplot(221), te=length(h); 
  plot(1:te,real(hb1),'bo-',1:te,real(h),'r+-'), grid
  legend('Estimated','Accurate')
  title('Real part of Channel');
  
  subplot(223),
  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