firbpsknocp.m

应用SIMO技术来构造高矩阵

clear all
close all
L=4; M=4; N=10; d=M+N; %L天线，M信道长度（阶数），N子载波数，d发送信号长度
j=sqrt(-1);    
T=1000;

h1=[-0.049+j*0.359 0.482-j*0.569 -0.556+j*0.587 1 -0.171+j*0.061].'; % channel
  h2=[0.443-j*0.0364 1 0.921-j*0.194 0.189-j*0.208 -0.087-j*0.054].';
  h3=[-0.221-j*0.322 -0.199+j*0.918 1 -0.284-j*0.524 0.136-j*0.19].';
  h4=[0.417+j*0.030 1 0.873+j*0.145 0.285+j*0.309 -0.049+j*0.161].';
h=[h1;h2;h3;h4];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=exp(-j*2*pi/d); %W=14*14
for i=1:d;
    for s=1:d;
W(i,s)=w^[(-1)*(i-1)*(s-1)];
    end 
end 
In=eye(N);
Ip=In(2*N-d+1:N,:);
Tcp=[Ip' In']';%40*36 

for dB=1:25

H1=toeplitz([h1(5);0;0;0;0;0;0;0;0;0],[h1(5) h1(4) h1(3) h1(2) h1(1) 0 0 0 0 0 0 0 0 0 ]);%n*n+m   %10*14
H2=toeplitz([h2(5);0;0;0;0;0;0;0;0;0],[h2(5) h2(4) h2(3) h2(2) h2(1) 0 0 0 0 0 0 0 0 0 ]);
H3=toeplitz([h3(5);0;0;0;0;0;0;0;0;0],[h3(5) h3(4) h3(3) h3(2) h3(1) 0 0 0 0 0 0 0 0 0 ]);
H4=toeplitz([h4(5);0;0;0;0;0;0;0;0;0],[h4(5) h4(4) h4(3) h4(2) h4(1) 0 0 0 0 0 0 0 0 0 ]);
H=[H1;H2;H3;H4]; %40*14

for i=1:T
    
 %s=sign(rand(1,d)-0.5)+2*sign(rand(1,d)-0.5);  % 16 QAM symbols
%s=s+sqrt(-1)*(sign(rand(1,d)-0.5)+2*sign(rand(1,d)-0.5));  
t=rand(1,d) ;%生成32维BPSK码
    for m=1:d
       if t(1,m)>(0.5)
          t(1,m)=1;
       else  t(1,m)=-1;
       end
    end 

Y(:,i)=H*W*(t.');%40*1

SNR=[];
n=randn(size(Y(:,i)))+sqrt(-1)*randn(size(Y(:,i)));%产生随机噪声
  n=n/norm(n)*10^(-dB/10)*norm(Y(:,i));
  SNR=[SNR 10*log10(norm(Y(:,i))/norm(n))];
  Y(:,i)=Y(:,i)+n;
end %Y=40*1000
%%%%%%%%%%%%%%%%%%%%%%%%
Ry=Y*Y'/T; %40*40
%Ry=1/K;
   [u,s,v]=svd(Ry);%奇异值分解
  
    Q=zeros(L*(M+1), L*(M+1));             % Construct matrix A (in Q)   %%%%%%%%20*20
for i=d+1:L*N,%15:40
  Vm=zeros(L*(M+1), M+N);%20*14
  for j2=1:(M+1),%1:5
    for k=1:L,%1:4
      Vm((k-1)*(M+1)+j2, 6-j2:(5-j2+N))=u((k-1)*N+1:k*N, i).';%u=40*40   构造成了和H矩阵模型结构相同的矩阵
    end
  end
  Q=Q+Vm*Vm';
end

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

  %%%%%%%%%%%%%%%
  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,'-*')