clear all;

close all;

MATLAB_VERSION_NUMBER=version;   % Because the implementation of function "xcorr" are different
% in MATLAB 5 and MATLAB 6, this program is designed in MATLAB5,
% a small modification is needed to run well in MATLAB6.

L = 4;                         % Channel length, including h(0).
m=2;
n=2;

L_diff=0;

Le=L+L_diff;

NF=128;
Monte_Carlo_Runs=10;
SNR=30;
N=1024*16;
channel_runs=1;
realcumulant=0;
ADD_NOISE=1;
draw_fig=1;
% % % %
real_channel=1;
Roots_Amplitude=1.50;

rolloff=0.11;

right_channel=ones(channel_runs,1);

F=zeros(m,n,2*NF);

H_est1=zeros(m,n,NF);

H_est_even=zeros(m,n,NF/2);

timSPD=zeros(Monte_Carlo_Runs);


for channel=1:70
    channel

    ONMSE=[];
    ONMSEFD=[];
    BER=[];

    for SNR=[30:6:30]

        h=zeros(L,m,n);
        I=[-round(L/2)+1:round(L/2)];

        cond_delta=100;

        randn('state',(channel+80));

        while cond_delta>40

            for ii=1:m
                for jj=1:n
                    h(:,ii,jj)=randn*rcos((I-2)*0.25,rolloff,1.0)+randn*rcos((I-1)*0.25,rolloff,1.0);
                    h(:,ii,jj)=h(:,ii,jj)./max(abs(h(:,ii,jj)));
                end
            end

            for ii=1:m
                for jj=1:n
                    H(ii,jj,:)=fft(h(:,ii,jj),NF);
                end
            end
            cond_delta=cond(H(:,:,NF));
        end

        h_long=zeros(Le,m,n);

        h_long(1:L,:,:)=h;

        for ii=1:m
            for jj=1:n
                denom(ii,jj)=sum((squeeze(h(:,ii,jj))).^2);
                denomF(ii,jj)=sum(abs(H(ii,jj,1:2:NF)).^2);
            end
        end

        NMSE=zeros(m,n,Monte_Carlo_Runs);
        NMSEFD=zeros(m,n,Monte_Carlo_Runs);

        for iloop=1:Monte_Carlo_Runs

            SNR
            iloop

            counter=iloop+80;

            rand('seed',counter)
            s = -log(rand(n,N));      %%% Single side Exponential distributed real signal

            s = s-(mean(s'))'*ones(1,N);
            s=diag(1./std(s,0,2))*s;       %% Input Signal
            x = zeros(m,N);          %%% Observed mixture , received signal
            for ii=1:m
                for jj=1:n
                    x(ii,:) = x(ii,:) + filter(h(:,ii,jj),1,s(jj,:));
                end
            end
            x = x-(mean(x'))'*ones(1,N);

            for ii=1:m
                std_x(ii)=std(x(ii,:));
            end

            if ADD_NOISE
                for ii=1:m
                    rand('state',sum(100*clock)*rand);
                    noise=randn(1,N);
                    noise=noise-mean(noise);
                    power_coefficient = (10^(SNR/20))*std(noise)/std(x(ii,:));
                    x(ii,:) = x(ii,:)+noise/power_coefficient;
                end
            end


            for ii=1:m
                x(ii,:) = std_x(ii)*x(ii,:)/std(x(ii,:));
            end

            %%% The received signal are generated by now.

            % Estimate the Bispectrum
            [pow, pow_true]=estCorr(NF,H,h,iloop+80,SNR,L_diff);
            [Bis, Bis_true, Bis_true1]=estCum1(NF,Le,x,h,H);

            if realcumulant
                Bis=Bis_true;
                pow=pow_true;
                clear Bis_true;
            end

            delta=11;
            R=1;

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

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

            error=h_est-h_long


            %-------------------------------NMSEFD PART-------------------

            for ii=1:m
                for jj=1:n
                    for kk=1:NF/2
                        tmp=(1/denomF(ii,jj)) *(abs(H_est_even(ii,jj,kk))-abs(H_even(ii,jj,kk)))^2;
                        NMSEFD(ii,jj,iloop)=NMSEFD(ii,jj,iloop)+tmp;
                    end

                    for ll=1:Le
                        tmp=(1/denom(ii,jj))*abs(h_est(ll,ii,jj)-h_long(ll,ii,jj))^2;
                        NMSE(ii,jj,iloop)=NMSE(ii,jj,iloop)+tmp;

                    end
                end
            end
            ONMSEFD_mc(iloop)=(1/(n*m)) * sum(sum(squeeze(NMSEFD(:,:,iloop))));    %% average over Ni and No
            ONMSE_mc(iloop)=(1/(n*m)) * sum(sum(squeeze(NMSE(:,:,iloop))));    %% average over Ni and No
        end     %%%% end of MC
        ONMSEFD=[ONMSEFD,sum(ONMSEFD_mc)/Monte_Carlo_Runs];
        ONMSE=[ONMSE,sum(ONMSE_mc)/Monte_Carlo_Runs];
    end    %%%% end of SNR
    ONMSEFD_total(:,channel)=ONMSEFD;
    ONMSE_total(:,channel)=ONMSE;
end   %%%% end of channel loops