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;
N=1024*8;

NF=128;
Monte_Carlo_Runs=1;
SNR=30;
channel_runs=1;
realcumulant=0;
ADD_NOISE=1;
draw_fig=1;
% % % %
real_channel=1;
estimate_phase=1;
Roots_Amplitude=1.50;

%%%%% set to 1 for BPSK,0 for 4QAM
BPSK=1;

rolloff=0.11;

right_channel=ones(channel_runs,1);

K=zeros(n,n,2*NF);

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:channel_runs
    channel

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

    for SNR=[30:5: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

            %            Gernerate the input signals
            s=zeros(n,N);
            counter=iloop+80;

            if BPSK
                s=rand(n,N)-0.5;
                s = s-(mean(s'))'*ones(1,N);
                s=diag(1./std(s,0,2))*s;       %% Input Signal
                s(find(s>=0))=1;
                s(find(s<0))=-1;
            else
                for ii=1:n
                    s(ii,:)=QAM(N,2,counter,counter*11);
                    counter=counter+1234;
                end
            end

            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

            % received signals generated here

            delta=13;

            k2=0;
            k3=-13;

            % Estimate the Trispectrum
            [Tris,Tris_true]=estTris4(NF,Le,k2,k3,x,n,h,H);

            if realcumulant
                Tris=Tris_true;
                clear Tris_true;
            end

            % Estimate the Channel and the input signals

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

            error=h_est-h_long

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

            %-------------------------------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