clear all;
close all;
%生成一個1*64的矩陣
for N=1:100
mod_in=randint(1,64);
%子載波個數16循環前綴個數4
%QPSK調制
mod_out=zeros(1,length(mod_in)/2); 
R=reshape(mod_in,2,length(mod_in)/2);  
   T=R' ;   
   B2D=bi2de(T,'left-msb')+1;  
   Temp=[-1-j  -1+j  1-j   1+j]; 
    for i=1:length(mod_in)/2 
         mod_out(i)=Temp(B2D(i))/sqrt(2);%歸一化
      end 
%串并變換
E=reshape(mod_out,2,length(mod_out)/2);
%IFFT
f=ifft(E);
%并串
g=reshape(f,1,32);
for i=1:32
    g(i)=real(g(i));
end
a=g'*g;
if N==1
    l=a;
end
l=(a+l)/2;
end

%信道
m=[1 0.4 -0.7 0.5 0.1]; %5個
m1=zeros(1,4);
m2=zeros(1,4);
for i=1:4
    m1(i)=m(i);
end
m2(1)=m1(1);
c0=toeplitz(m1,m2);


m3=zeros(1,4);
m4=zeros(1,4);
u=5;
for i=2:5
    u=u-1;
    m4(u)=m(i);
end
m3(1)=m4(1);
c1=toeplitz(m3,m4);
c3=zeros(4,4);

H=[c0 c3 c3 c1 c3 c3 c3 c3;c1 c0 c3 c3 c3 c3 c3 c3;c3 c1 c0 c3 c3 c3 c3 c3;c3 c3 c1 c0 c3 c3 c3 c3 ;c3 c3 c3 c1 c3 c3 c3 c0;c3 c3 c3 c3 c0 c3 c3  c1;c3 c3 c3 c3 c1 c0 c3 c3;c3 c3 c3 c3 c3  c1 c0 c3;c3 c3 c3 c3 c3 c3 c1 c0];

%接收信號自相關矩陣

out1=H*l*H';
%加高斯白噪聲
snr=4;
SymCh = awgn(out1,snr,'measured');
%svd分解
[U,S,V] =svd(SymCh);
rank_l=rank(l);
b=U;
for i=1:rank_l;
    b(:,1)=[];
end
%b為噪聲子空間特征向量36*4

%生成矩陣
A=zeros(1,4);
%b的第1列生成的A矩陣
for z=1:4:33
m3=zeros(1,4);
m4=zeros(1,4);
u=5;
e=z;
for n=1:4
    u=u-1;
    m4(u)=b(e,1);
    e=e+1;
end
m3(1)=m4(1);
temp=toeplitz(m3,m4);
temp=fliplr(temp);
M=[temp;A];
if z==1
    A11=M;
elseif z==5
    A12=M;
elseif z==9
    A13=M;
elseif z==13
    A14=M;
elseif z==17
    A15=M;
elseif z==21
    A16=M;
elseif z==25
    A17=M;
elseif z==29
    A18=M;
elseif z==33
    A19=M;
end
end
 
%b的第2列生成的A矩陣
for z=1:4:33
m3=zeros(1,4);
m4=zeros(1,4);
u=5;
e=z;
for n=1:4
    u=u-1;
    m4(u)=b(e,2);
    e=e+1;
end
m3(1)=m4(1);
temp=toeplitz(m3,m4);
temp=fliplr(temp);
M=[temp;A];
if z==1
    A21=M;
elseif z==5
    A22=M;
elseif z==9
    A23=M;
elseif z==13
    A24=M;
elseif z==17
    A25=M;
elseif z==21
    A26=M;
elseif z==25
    A27=M;
elseif z==29
    A28=M;
elseif z==33
    A29=M;
end
end

%b的第3列生成的A矩陣
for z=1:4:33
m3=zeros(1,4);
m4=zeros(1,4);
u=5;
e=z;
for n=1:4
    u=u-1;
    m4(u)=b(e,3);
    e=e+1;
end
m3(1)=m4(1);
temp=toeplitz(m3,m4);
temp=fliplr(temp);
M=[temp;A];
if z==1
    A31=M;
elseif z==5
    A32=M;
elseif z==9
    A33=M;
elseif z==13
    A34=M;
elseif z==17
    A35=M;
elseif z==21
    A36=M;
elseif z==25
    A37=M;
elseif z==29
    A38=M;
elseif z==33
    A39=M;
end
end

%b的第4列生成的A矩陣
for z=1:4:33
m3=zeros(1,4);
m4=zeros(1,4);
u=5;
e=z;
for n=1:4
    u=u-1;
    m4(u)=b(e,4);
    e=e+1;
end
m3(1)=m4(1);
temp=toeplitz(m3,m4);
temp=fliplr(temp);
M=[temp;A];
if z==1
    A41=M;
elseif z==5
    A42=M;
elseif z==9
    A43=M;
elseif z==13
    A44=M;
elseif z==17
    A45=M;
elseif z==21
    A46=M;
elseif z==25
    A47=M;
elseif z==29
    A48=M;
elseif z==33
    A49=M;
end
end

%b的第1列生成的B矩陣
for z=1:4:33
m1=zeros(1,4);
m2=zeros(1,4);
u=1;
e=z;
for n=1:4
    m1(u)=b(e,1);
    u=u+1;
    e=e+1;
end
m2(1)=m1(1);
temp=toeplitz(m1,m2);
temp=fliplr(temp);
M=[A;temp];
if z==1
    B11=M;
elseif z==5
    B12=M;
elseif z==9
    B13=M;
elseif z==13
    B14=M;
elseif z==17
    B15=M;
elseif z==21
    B16=M;
elseif z==25
    B17=M;
elseif z==29
    B18=M;
elseif z==33
    B19=M;
end
end

%b的第2列生成的B矩陣
for z=1:4:33
m1=zeros(1,4);
m2=zeros(1,4);
u=1;
e=z;
for n=1:4
    m1(u)=b(e,2);
    u=u+1;
    e=e+1;
end
m2(1)=m1(1);
temp=toeplitz(m1,m2);
temp=fliplr(temp);
M=[A;temp];
if z==1
    B21=M;
elseif z==5
    B22=M;
elseif z==9
    B23=M;
elseif z==13
    B24=M;
elseif z==17
    B25=M;
elseif z==21
    B26=M;
elseif z==25
    B27=M;
elseif z==29
    B28=M;
elseif z==33
    B29=M;
end
end
%b的第3列生成的B矩陣
for z=1:4:33
m1=zeros(1,4);
m2=zeros(1,4);
u=1;
e=z;
for n=1:4
    m1(u)=b(e,3);
    u=u+1;
    e=e+1;
end
m2(1)=m1(1);
temp=toeplitz(m1,m2);
temp=fliplr(temp);
M=[A;temp];
if z==1
    B31=M;
elseif z==5
    B32=M;
elseif z==9
    B33=M;
elseif z==13
    B34=M;
elseif z==17
    B35=M;
elseif z==21
    B36=M;
elseif z==25
    B37=M;
elseif z==29
    B38=M;
elseif z==33
    B39=M;
end
end
%b的第4列生成的B矩陣
for z=1:4:33
m1=zeros(1,4);
m2=zeros(1,4);
u=1;
e=z;
for n=1:4
    m1(u)=b(e,4);
    u=u+1;
    e=e+1;
end
m2(1)=m1(1);
temp=toeplitz(m1,m2);
temp=fliplr(temp);
M=[A;temp];
if z==1
    B41=M;
elseif z==5
    B42=M;
elseif z==9
    B43=M;
elseif z==13
    B44=M;
elseif z==17
    B45=M;
elseif z==21
    B46=M;
elseif z==25
    B47=M;
elseif z==29
    B48=M;
elseif z==33
    B49=M;
end
end
%生成C矩陣
C11=A11+B12;    
C12=A12+B13;
C13=A13+B14;
C14=A14+B15;
C15=A15+B16;
C16=A16+B17;
C17=A17+B18;
C18=A18+B19;

C21=A21+B22;    
C22=A22+B23;
C23=A23+B24;
C24=A24+B25;
C25=A25+B26;
C26=A26+B27;
C27=A27+B28;
C28=A28+B29;

C31=A31+B32;    
C32=A32+B33;
C33=A33+B34;
C34=A34+B35;
C35=A35+B36;
C36=A36+B37;
C37=A37+B38;
C38=A38+B39;   

C41=A41+B42;    
C42=A42+B43;
C43=A43+B44;
C44=A44+B45;
C45=A45+B46;
C46=A46+B47;
C47=A47+B48;
C48=A48+B49;
%生成的G矩陣
G1=[C11,C12,C13,C14+B11,C15+A19,C16,C17,C18];
G2=[C21,C22,C23,C24+B21,C25+A29,C26,C27,C28];
G3=[C31,C32,C33,C34+B31,C35+A39,C36,C37,C38];
G4=[C41,C42,C43,C44+B41,C45+A49,C46,C47,C48];
%計算Q矩陣
Q=G1*G1'+G2*G2'+G3*G3'+G4*G4';

%生成代價矩陣
ub=[1,1,1,1,1];
lb=[-1,-1,-1,-1,-1];
%計算h在范數為1的情況下使h'*Q*h為最小情況下的h
x = quadprog(Q,[],[],[],[],[],lb,ub) ;
%計算均方誤差
p=0;
for i=1:5
    t=(m(i)-x(i))^2;
    p=p+t;
end
q=0
for i=1:5
    r=m(i)^2;
    q=q+r;
end
aaaa=p/q