%OFDM系統(tǒng)基于子空間的盲信道估計(jì)

clear all
close all
clc
Mt=1;Mr=1;%the number of transmit and received antennas 
N=48;   %the number of subcarrier
nos=200;%number of symbol
P=12;   %length of CP
Q=N+P;  %the totle length of an OFDM symbol
J=3;    %collect J consecutively received OFDM symbols
L=4;    %the order of channel
Nm=1;



%構(gòu)造W矩陣
j=sqrt(-1);
W1=zeros(N,N);
for i1=1:N
    for j1=1:N
        W1(j1,i1)=exp(j*2*pi*(N-i1)*(j1-1)/N)/sqrt(N);
    end
end
W=[W1,W1(:,1:P)].';
Ij=eye(J);
Imt=eye(Mt);
W2=kron(Ij,W);
W3=kron(W2,Imt);
W4=conj(W)*(W.');
A=kron(Ij,W4);


%------------------------------------------------------------%

%-----------------------go through the channel-----------------%
mh=[-0.049+j*0.359  0.482-j*0.569  -0.556+j*0.587 1              -0.171+j*0.061];  %h11
    %0.443-j*0.0364  1              0.921-j*0.194  0.189-j*0.208  -0.087-j*0.054;  %h12
    %-0.221-j*0.322  -0.199+j*0.918 1              -0.284-j*0.524 0.136-j*0.19;    %h13
    %0.417+j*0.030   1              0.873+j*0.145  0.285+j*0.309  -0.049+j*0.161]; %h14

H=mh(:);%將mh變?yōu)榱邢蛄?
%構(gòu)造H1
H1=zeros((J*Q-L)*Mr,J*Q*Mt);

H1(1:Mr,:)=[H(1:Mr,:) H(Mr+1:2*Mr,:) H(2*Mr+1:3*Mr,:) H(3*Mr+1:4*Mr,:) H(L*Mr+1:(L+1)*Mr,:) zeros(Mr,(J*Q-L-1)*Mt)];
for i=1:J*Q-L-1
    htmp=H1((i-1)*Mr+1:i*Mr,:);
    H1(i*Mr+1:(i+1)*Mr,:)=[htmp(:,(J*Q-1)*Mt+1:J*Q*Mt),htmp(:,1:(J*Q-1)*Mt)];
end
%------------------------進(jìn)行子空間估計(jì)--------------------------------%

NRMSE=zeros(1,9);SNR=[];a=[];


SNR1=0:5:40;
for snr=0:5:40
     Herr=0;
    for i1=1:Nm
        Rr=zeros((J*Q-L)*Mr,(J*Q-L)*Mr);
        %------------------產(chǎn)生OFDM符號(hào)-----------------------%
        b0=round(rand(1,2*N*nos));
        %QPSK modulation
        b=zeros(1,N*nos);
        b1=reshape(b0,2,N*nos).';
        b2=bi2de(b1,2,'left-msb');
        b3=[1+j 1-j -1+j -1-j];
        b=b3(b2+1)/sqrt(2);

        %串并變換
        d1=b(:);
        
        for i=1:nos-J+1
            d4(:,i)=d1((i-1)*N*Mt+1:(i+J-1)*N*Mt);%連取J個(gè)符號(hào)
        end
        s=W3*d4;
        %信號(hào)經(jīng)過信道
        r1=H1*s;
        
        n=rand(size(r1))+sqrt(-1)*rand(size(r1));
        n=n/norm(n)*norm(r1)*10^(-snr/20);
        SNR=[SNR 20*log10(norm(r1)/norm(n))];
        r=r1+n;
        
        %對(duì)Rr進(jìn)行平滑

        Rr=r*r'/(nos-J+1); 
        %特征值分解
        [U,S,V]=svd(Rr);
        %[U,V]=eigs(Rr,(J*Q-L)*Mr-J*N*Mt,'sm');
        %構(gòu)造C矩陣
        C=zeros((L+1)*Mr,(L+1)*Mr);
        %for k=1:(J*Q-L)*Mr-J*N*Mt
        for k=J*N*Mt+1:(J*Q-L)*Mr
            uk=U(:,k);
            vk1=reshape(uk,Mr,J*Q-L);
            vk=zeros((L+1)*Mr,J*Q);
            vk(1:Mr,:)=[vk1 zeros(Mr,L)];
            for i3=1:L
                vktmp=vk((i3-1)*Mr+1:i3*Mr,:);
                vk(i3*Mr+1:(i3+1)*Mr,:)=[vktmp(:,J*Q) vktmp(:,1:J*Q-1)];
            end
            C=C+vk*A*vk';
        end
        [UC SC VC]=svd(C);
        %[UC VC]=eigs(C,Mt,'sm');
        hr0=UC(:,(L+1)*Mr-Mt+1:(L+1)*Mr);
        a=lscov(hr0,H);
        a1=mean(H./hr0);
        %a=lscov(hr0,H);
        Hr=hr0*a;
        Herr1=H-Hr;
        
        for i4=1:Mt
            Herr=Herr+(norm(Herr1(:,i4)))^2/(norm(H(:,i4)))^2;
        end
    end
    NRMSE(snr/5+1)=sqrt(Herr/(Nm*Mt*Mr*(L+1)));
end
semilogy(SNR1,NRMSE,'-*r')
xlabel('SNR/dB');ylabel('NRMSE');grid on;