function [h_est]=ComputeISPD3rd(real_channel,draw_fig,NF,L_diff,Bis,h,H,N,delta,R)

%%%%% This function calculate the estimated MIMO channel gain for given
%%%%% Bispectrum 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

Cum0 = squeeze(Bis(mod(-R*delta,NF)+1,mod(-delta,NF)+1,:,:,:));
Cum1 = squeeze(Bis(mod(-R*delta-delta,NF)+1,mod(-delta,NF)+1,:,:,:));

%Parafac decomposition
[estA0,estB0,estC0,fit,it]=comfac(Cum0,n);

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

estA0=estA0*pinv(PERM);
estB0=estB0*pinv(PERM);
estC0=estC0*pinv(PERM);

for r=r_index-r_index(1)

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

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

        UA = [];

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

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


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

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')
        end
    end

end

%%% solve for phase ambiguity.........

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_even(2,ii,r_index_even))'*squeeze(H_est_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