% 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 displaydB=15; T=1000;  % SNR, sample amountL=4; M=4; N=5; d=M+N; % L: antenna #. M: channel length. N: smoothing. j=sqrt(-1);           % d: equalization delaymh=[-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 symbolss=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 signalsfor 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).';  endend if echoF SNR=SNR, endss=std(s)^2;   sv=std(v)^2;  %%%%%%%%%%%%%% subspace method beginRx=X*X'/TN;                 % calculate correlation matrix[U0,S0,V0]=svd(Rx);         % SVD to find null subspacefor i=L*N:-1:1,   if S0(i-1,i-1)-S0(i,i)>S0(i,i), break; end  endi=d+1;                      % check rank of null subspace%i=rank(S0)+1;if echoF, d=i-1, else d=i-1; end   % display ranksigma=0;for i=i:L*N, sigma=sigma+S0(i,i); end   % remove noisesigma=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 SVDhb=U1(:,L*(M+1));           % channel estimation%%%%%%% Compare channel estimation MSEhb_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, endif 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 resultsH=zeros(L*N, M+N);  %% channel matrixfor j=1:N,   for k=1:L,      H((k-1)*N+j, j:(j+M))=hb1((k-1)*(M+1)+1:k*(M+1)).';   endendY=H'*U0(:,1:d)*inv(S0(1:d, 1:d)-sigma*eye(d))*U0(:,1:d)'*X; % zero-forcing equalizergd=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],    endif echoF,  subplot(222), plot(s,'ro'), grid, title('Transmitted Symbols')  subplot(224), plot(Y(round(d/2),:),'ro'), grid,   title('Estimated Symbols')end