cpofdm.m

基于子空间的盲信道估计算法

clear all
j=sqrt(-1);
N=32,p=40;e=1000;%子载波的个数（N）,P=N+前缀，e发送的信号个数
h(1)=0.158-0.664j;
h(2)=-0.198+0.267j;
h(3)=-0.325+0.197j;
h(4)=-0.378-0.245j;
h(5)=-0.278-0.003j;
h=[h(1) h(2) h(3) h(4) h(5)]'
H=Toeplitz([h(5);0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0],[h(5) h(4) h(3) h(2) h(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]);
%76*80
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Im=eye(2);
In=eye(32);
Ip=In(25:32,:);%8*32
Tcp=[Ip' In']';%40*32   作用是用来加前缀的
G=kron(Im,Tcp);%80*64 W临时使用   发送端是两个信号的情况
A=H*G;%76*64 
 
for dB=1:5:40
%%%%%%%%%%%%%%%%%%%%%%%%%%   做傅立叶变换（生成一个N*N的用于做傅立叶变换的矩阵）
for i1=1:e
w=exp(-j*2*pi/N);
for i=1:N;
    for s=1:N;
W(i,s)=w^[(-1)*(i-1)*(s-1)];
    end 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%
X1=rand(32,1);%生成32维BPSK码
for m=1:32
       if X1(m,1)>(0.5)
        X1(m,1)=1;
    else  X1(m,1)=-1;
    end
end
X1=W*X1;
  X2=rand(32,1) ;%生成32维BPSK码
for m=1:32
       if X2(m,1)>(0.5)
        X2(m,1)=1;
    else  X2(m,1)=-1;
    end
end
 X2=W*X2; 
X=[X1;X2];%%生成64维BPSK码   64*1
Y(:,i1)=A*X;%76*1;
SNR=[];
n=randn(size(Y(:,i1)))+sqrt(-1)*randn(size(Y(:,i1)));%产生随机噪声
  n=n/norm(n)*10^(-dB/10)*norm(Y(:,i1));
  SNR=[SNR 10*log10(norm(Y(:,i1))/norm(n))];
  Y(:,i1)=Y(:,i1)+n;
end
%%%%%%%%%%%%%%%%%%%%%%%%
Ry=Y*Y'/e;%Y为76*1000   Ry为76*76
%Ry=1/K;
   [u,s,v]=svd(Ry);%奇异值分解
   Un=u(:,65:76);%76*12
   Q=0;
for i=1:12
    gi=Un(:,i).';    %1*76
    Gi=[0 0 0 0 gi;0 0 0 gi 0;0 0 gi 0 0;0 gi 0 0 0;gi 0 0 0 0];%5*80   Gi = (L+1) * N
  
    qi=(Gi*G)*(Gi*G)';%5*5
    Q=Q+qi;


end
   [U1,S1,V1]=svd(Q);          % solve equation Ah=0 by SVD
hb=U1(:,5)      

  %%%%%%%%%%%%%%%
  hb_h=mean(hb./h);            
hb1=hb/hb_h
squ_err_h(dB)=sqrt((h-hb1)'*(h-hb1))/sqrt(h'*h);
end
figure(1)
plot(squ_err_h,'-*')