% synthcoh.m
%
% srv, fm 6/5/02 

t=0:1774;
rand1=rand(1,1775);
rand2=rand(1,1775);
x=1*cos(2*pi*90/360+2*pi/21*t)+4*cos(+2*pi*00/360-2*pi/10*t)+1.5*(-1+2*rand1);
y=1*cos(2*pi*60/360+2*pi/21*t)+4*cos(-2*pi*70/360-2*pi/10*t)+1.5*(-1+2*rand2);

figure(1)
clf
plot(x,'b')
hold
plot(y,'r')
legend('series 1, x','series 2, y')
title('Raw mystery time series')
xlabel('sample number')
ylabel('amplitude')
print -djpeg raw.jpg

axis([0 100 -8 8])
grid on;
title('First 100 samples of raw mystery data series')
print -djpeg rawzoom.jpg


% start with psd/csd approach

% calculate cross/power spectra
% note that if Fs~=1, then be careful interpreting f-axis labels when plotting
fs=1;      % sampling frequency
nfft=256;  % width of subseries (# of samples)
ww=256;    % width of hanning window applied to each subseries
over=0;
[gxx,gxxc,f]=psd(x,nfft,fs,ww,over,.95);
gyy=psd(y,nfft,fs,ww,over,.95);
gxy=csd(x,y,nfft,fs,ww,over,.95);

% plot cross/power spectra

figure(2)
clf
psd(x,nfft,fs,ww,over,.95);
title('Power spectral density estimates for x')
print -djpeg psdx.jpg
clf
psd(y,nfft,fs,ww,over,.95);
title('Power spectral density estimates for y')
print -djpeg psdy.jpg
clf
csd(x,y,nfft,fs,ww,over,.95);
title('Cross power spectral density between x and y')
print -djpeg csdxy.jpg

% calculate coh and phase, then plot
cohxy=(abs(gxy).^2)./(gxx.*gyy);
phsxy= atan2(-imag(gxy),real(gxy)) * 180/pi;

figure(3)
clf
title('Coherence and phase (>95%) between x and y')
% EDOF=(length(timeseries)/windowhalfwidth)*8/3=5.33
% So, 95% confidence level = gamma_95 = 1-(0.05)^(1/(5.33-1))
edof=(8/3)*(length(x)/ww/2)
gamma_95=1-(0.05)^(1/(edof-1)) 
conf95=find(cohxy>gamma_95);
[ax,hc,hp]=plotyy(f/fs,cohxy,f(conf95)/fs,phsxy(conf95));
hold on; 
%grid on;
plot([min(f) max(f)],[gamma_95 gamma_95],'k','LineWidth',[1])
set(hc,'Color',[0 0 1],'LineStyle','-')
set(hp,'MarkerFaceColor',[0 1 0],'LineStyle','none','marker','*')
xlabel('Frequency, hr^-^1')
title('Coherence and phase (>95%) between x and y')
set(get(ax(1),'YLabel'),'string','Coherence^2');
set(get(ax(2),'YLabel'),'string','Phase');
y1tls=['0.0';'0.2';'0.4';'0.6';'0.8';'1.0'];
set(ax(1),'YTick',[0:0.2:1],'YTickLabel',y1tls);
y2tls=['-180';'-135';' -90';' -45';'   0';'  45';'  90';' 135';' 180'];
set(ax(2),'YTick',[-180:45:180],'YTickLabel',y2tls);
text(0.005,.95,sprintf('N=%g; NFFT=%g; Windowwidth=%g',length(x),nfft,ww));
text(0.005,.9,sprintf('EDOF=%g; overlap=%g',edof,over));
%print -depsc2 cohphasexy.eps2
print -djpeg cohphasexy.jpg

% now try fft approach

fx=fft(x);
fy=fft(y);

