close all
clear all
n=256;
nhlf=n/2;
nhlfp=nhlf+1;
nhlfm=nhlf-1;
nbi=n;
nbihlfp=nbi/2 +1;

im=zeros(1,n);
kill=1:256;
im=rectpuls(kill-60,10) + 3*rectpuls(kill-100,10);
%im(32)=4;
%im(48)=10;
%im(80)=7;

abski=abs(fft(im));
phaseki=angle(fft(im));
figure(9);
kk=ifft(abski.*exp(i*unwrap(phaseki)));
plot(abs(kk));

figure(1)
plot(im)
title('Initial Signal');
set(gca,'XLimMode','manual','XLim',[0 n-1]);
pause(3);




fim = fft(im);
sfim=fftshift(fim);
figure(2);
plot(abs(sfim));
title('FFt of Signal');
set(gca,'XLimMode','manual','XLim',[0 n-1]);
pause(3);

%bispect=zeros(nbi,nbi);
%for v=-nhlf:nhlfm
%  for u=-nhlf:nhlfm 
 %     if abs(u+v) <= nhlfm
  %           bispect(nbihlfp +v, nbihlfp +u)= sfim(nhlfp+u).*sfim(nhlfp+v) .*conj(sfim(nhlfp+u+v));
       
             
%end
 % end
%end
bispect=zeros(nbi,nbi);
for v=1:256
  for u=1:256 
       
             bispect(v,u)= sfim(u).*sfim(v) .*conj(sfim(1+mod((u+v),256)));
end
end

bispect=fftshift(bispect)
bispectud=(bispect);
figure(3);
scale=max(max(bispectud));
bispectud=bispectud/scale;
imshow(abs(bispectud));
title('Bispectrum of Signal');
pause(3);




phaseko(1)=0;
for v=2:256
    angle(bispectud(abs(v-1)+1,1))
     
    phaseko(v)=mod(phaseko(v-1)-angle(bispectud(1,v-1)),2*pi);
end


figure(10);
kk=ifft(abski.*exp(i*unwrap(phaseko)));
plot(fftshift(abs(kk)));


bispectorig=fftshift(bispect);
tricorr=ifft2(bispectorig);
figure(4);
stricorr=fftshift(tricorr);
stricorr(nbihlfp,nbihlfp)=0;

nbig=4;
temp=zeros(nbi+nbig-1,nbi+nbig-1);
for yind=1:nbig
 for xind=1:nbig
   temp(yind:yind+nbi-1,xind:xind+nbi-1)=temp(yind:yind+nbi-1,xind:xind+nbi-1)+stricorr;
end
end

imdisp=temp(1:nbi,1:nbi);
imdisp=abs(imdisp);
maxval=max(max(imdisp));
minval=min(min(imdisp));
scale=maxval-minval;
imdisp=(imdisp-minval)/scale;

imshow(imdisp);
title('Triple Correlation');
pause(3);

x=[0:nbi-1];
y=x;

figure(5)
mesh(x,y,imdisp)
title('Triple Correlation');
set(gca,'XLimMode','manual','XLim',[0 n-1]);
set(gca,'XLimMode','manual','XLim',[0 n-1]);


