computeispd4th.m

基于卷积信号的MIMO系统盲信号估计

function [h_est]=ComputeISPD4th(real_channel,draw_fig,NF,L_diff,k2,k3,Tris,h,H,N)

%%%%% This function calculate the estimated MIMO channel gain for given
%%%%% Trispectrum and other coefficients

[L,m,n]=size(h);

Le=L+L_diff;

%%%% FULL Frequency
r_index=[1:NF+1];

r_index_even=[1:NF/2];

%%% Initialize for the SPD

delta=-(k2+k3);
R=1;

Cum3=squeeze(Tris(mod(-R*delta,NF)+1,:,:,:,:));
Cum3ex=squeeze(Tris(mod(-(R+1)*delta,NF)+1,:,:,:,:));


% Parafac decomposition

% estA0=H(:,:,mod((R+1)*delta,NF)+1);
% estB0=conj(H(:,:,mod(R*delta,NF)+1));
% estC0=H(:,:,mod(k2,NF)+1);
% estD0=conj(H(:,:,mod(-k3,NF)+1));

%[estA0,estB0,estC0,estD0]=QALSex1(Cum3,Cum3ex,n);
[estA0,estB0,estC0,estD0]=QALS(Cum3,n);

for r=r_index-r_index(1)

    Cum3r= squeeze(Tris(mod(-R*delta-r*delta,NF)+1,:,:,:,:));

    if r==0
        F(:,:,r+1)=estA0;
    else

        UA = [];

        for iii=1:m,
            for jjj=1:m,
                UA = [UA; squeeze(Cum3r(:,iii,jjj,:)).'];
            end
        end

        F(:,:,r+1) = (pinv( krp(conj(squeeze(F(:,:,r))),krp(estC0,estD0)) )*UA).';

    end
end


PERM=allign3(estC0,H(:,:,mod(k2,NF)+1));

%...compensating permulation and scalar matrices for comparison to H_even

%compensating for constant permutation

for rr=r_index-r_index(1)
    F(:,:,rr+1)=F(:,:,rr+1)*pinv(PERM);
end

for ii=0:NF-1
    H_est1(:,:,mod((R+ii+1)*delta,NF)+1)=F(:,:,ii+1);
end

H_even=H(:,:,1:2:end);

if draw_fig

    figure; clf

    for ii=1:m
        for jj=1:n
            subplot(m,n,jj+(ii-1)*n), plot(abs(squeeze(H_even(ii,jj,:)))); hold on
            plot(abs(squeeze(H_est1(ii,jj,1:2:end))),'r');
            title('amplitude of the H_{est} and H');
        end
    end


    figure; clf
    for ii=1:m
        for jj=1:n
            subplot(m,n,jj+(ii-1)*n), plot(unwrap(angle(squeeze(H_even(ii,jj,:))))); hold on
            plot(unwrap(angle(squeeze(H_est1(ii,jj,1:2:end)))),'r');
            title('phase of H_{est} and H before solve the phase ambiguity')
        end
    end

end


% Estimate the phase
phi2=zeros(n,n);

if m>=n    %%%%% Over-determined MIMO

    phi2=diag(squeeze(F(1,:,NF+1)))*pinv(diag(squeeze(F(1,:,1))));

    phi2=diag(diag((1/NF)*angle(phi2)));

    for r=0:NF-1

        F(:,:,mod(r,NF)+1)=F(:,:,mod(r,NF)+1) * expm(-sqrt(-1)*(r)*phi2);

        H_est1(:,:,mod((R+ii+1)*delta,NF)+1)=F(:,:,ii+1);

    end

else

    [mm,jj]=max( min( abs(squeeze(F(:,:,1))),[],2 ));

    phi2=diag(squeeze(F(jj,:,NF+1)))*pinv(diag(squeeze(F(jj,:,1))));

    phi2=diag(diag((1/NF)*angle(squeeze(phi2))));

    for r=0:NF-1

        F(:,:,mod(r,NF)+1)=F(:,:,mod(r,NF)+1) * expm(-sqrt(-1)*(r)*phi2);

        H_est1(:,:,mod((R+ii+1)*delta,NF)+1)=F(:,:,ii+1);

    end
end

H_est_even=H_est1(:,:,1:2:end);

%%%%% scalar ambiguity
Scalar_even=zeros(n,n);
for ii=1:n
    Scalar_even(ii,ii)=( abs(squeeze(H_even(2,ii,r_index_even))).'*abs(squeeze(H_est_even(2,ii,r_index_even)) )) / (abs(squeeze(H_est_even(2,ii,r_index_even))).' * abs(squeeze(H_est_even(2,ii,r_index_even)) + 0.01 ));
end


%compensating for constant Scarlar
for ii=r_index_even
    H_est_even(:,:,ii)=H_est_even(:,:,ii)*Scalar_even;
end


if draw_fig
    figure; clf
    for ii=1:m
        for jj=1:n
            subplot(m,n,jj+(ii-1)*n), plot(r_index_even,abs(squeeze(H_est_even(ii,jj,r_index_even)))); hold on
            plot(abs(squeeze(H_even(ii,jj,:))),'r');
            title('amplitude after solve the scalar ambiguity');
        end
    end

    figure; clf
    for ii=1:m
        for jj=1:n
            subplot(m,n,jj+(ii-1)*n), plot(r_index_even,unwrap(angle(squeeze(H_est_even(ii,jj,r_index_even))))); hold on
            plot(unwrap(angle(squeeze(H_even(ii,jj,:)))),'r');
            title('phase after solve the non-integer phase ambiguity');
        end
    end

end

for ii=0:NF-1
    Q(:,:,ii+1)=F(:,:,mod(ii-R-1,NF)+1);
end

%%%%%% downsmpling and solve the circular shift, proper allign with the true channel

for ii=1:m
    for jj=1:n

        [H_est_even(ii,jj,:),h_est(:,ii,jj)]=relign4(NF/2, squeeze(Q(ii,jj,1:2:NF)), squeeze(h(:,ii,jj)), real_channel, L_diff, delta);

    end
end


%%%%% scalar ambiguity
Scalar_even=zeros(n,n);
for ii=1:n
    Scalar_even(ii,ii)=( squeeze(H_est_even(2,ii,r_index_even))'*squeeze(H_even(2,ii,r_index_even))) / (squeeze(H_est_even(2,ii,r_index_even))' * squeeze(H_est_even(2,ii,r_index_even)));
end


%compensating for constant Scarlar
for ii=r_index_even
    H_est_even(:,:,ii)=H_est_even(:,:,ii)*Scalar_even;
end


if draw_fig
    figure; clf
    for ii=1:m
        for jj=1:n
            subplot(m,n,jj+(ii-1)*n), plot(r_index_even,abs(squeeze(H_est_even(ii,jj,r_index_even)))); hold on
            plot(abs(squeeze(H_even(ii,jj,:))),'r');
            title('amplitude after solve the scalar ambiguity')
        end
    end

    figure; clf
    for ii=1:m
        for jj=1:n
            subplot(m,n,jj+(ii-1)*n), plot(r_index_even,unwrap(angle(squeeze(H_est_even(ii,jj,r_index_even))))); hold on
            plot(unwrap(angle(squeeze(H_even(ii,jj,:)))),'r');
            title('phase after solve the integer phase ambiguity');
        end
    end

end