t=1:4*1024;
T=1024;
om=2*pi/32;
x=sin(om*t);
y=3*x+10*(rand(1,4*1024)-.5);;;
x=x-mean(x);
y=y-mean(y);


nfft=512;
nf=nfft/2+1;
h=hanning(nfft);
T=nfft;
F=24*[0:nf-1]/T;

Pxx=zeros(1,nfft);
Pyy=zeros(1,nfft);
Pxy=zeros(1,nfft);
nseg=floor(length(x)/nfft);
for iseg=1:nseg;
    i1=1+(iseg-1)*nfft;
    i2=iseg*nfft;
    xx=x(i1:i2).*h';
    yy=y(i1:i2).*h';
    Xx=fft(xx,nfft);
    Yy=fft(yy,nfft);
    Xx2=Xx.*conj(Xx);
    Yy2=Yy.*conj(Yy);
    Xy2=Yy.*conj(Xx);
    Pxx = Pxx + Xx2;
    Pyy = Pyy + Yy2;
    Pxy = Pxy + Xy2;
end
   
coher=(abs(Pxy).^2)./(Pxx.*Pyy);
coher=coher(1:nf);
plot(F,coher);
clim=1-.05^(1/(nseg-1));
hold on
plot([0 12],[clim clim],'k');