%       Calculate coherence function
%       Written by Marinna Martini
%
%       Returns the vector:
%       [H,PAxy,PAxx,PAyy,f]=coh(x,y,Nfft,dt);
%               where
%					 Nfft is the fft block size
%					 dt if the time between samples
%               H is the coherence,
%               PAxy is the averaged Cross Spectral Density
%               PAxx is the averaged Auto Spectral Density of x
%               PAyy is the averaged Auto Spectral Density of y
%               f is the corresponding matrix of frequencies
%
function [I,PAxy,PAxx,PAyy,freq]=coh(x,y,Nfft,dt);

n = max(size(x));
m = max(size(y));

if n~=m,
        error('*** Vectors are different lengths!');
        return;
end

xm=mean(x);
ym=mean(y);
x=x-xm;
y=y-ym;

if ~exist('Nfft')==1,
   N = input('Enter an even number of blocks to average:  ');
else
   N = Nfft;
end
if ~exist('dt') == 1,
   dt = input('Enter the time interval between samples:  ');
end
%fprintf('OK Please Wait....\n');
l = round(n/N);       % length of each block

% divide into blocks
for i=0:(N-1),
        a = i*l+1;
        b = (i+1)*l;
        for j=a:b,
                xx(i+1,j-l*i)=x(j);
                yy(i+1,j-l*i)=y(j);
        end
end

%  Get the properly normalized fourier coefficients
for i=1:N,
        X(i,:) = (2/l)*fft(xx(i,:)); % get the Fourier coefficients
        Y(i,:) = (2/l)*fft(yy(i,:)); % for each block
end

% Figure out the corresponding frequency matrix
df = 1/(l*dt);  % NOTE: use length of averaged block
Fn = 1/(2*dt);  % Nyquist frequency
freq = 0:df:Fn; % matrix of harmonics
fm = max(size(freq));

% Calculate the Cross Spectral Density
for i=1:N,
        Pxy(i,:) = .5*X(i,:).*conj(Y(i,:));        % do for each block
end
for i=1:l,
        PAxy(i) = sum(Pxy(:,i))/N;      % average the blocks
end
CAxy = real(PAxy);        % Co-spectra
QAxy = imag(PAxy);        % Quadrature spectra
PAxy = (1/df)*PAxy;       % Normalize

% plot the cross spectra
clg
subplot(223)
axis;
loglog(freq(2:fm),PAxy(2:fm));
title('Averaged Cross Spectra')


% Calculate the Auto Spectral Density
for i=1:N,
        Pxx(i,:) = .5*X(i,:).*conj(X(i,:));        % do for each block
        Pyy(i,:) = .5*Y(i,:).*conj(Y(i,:));        % do for each block
end
for i=1:l,
        PAxx(i) = sum(Pxx(:,i))/N;      % average the blocks
        PAyy(i) = sum(Pyy(:,i))/N;      % average the blocks
end
PAxx = (1/df)*PAxx;         % Normalize
PAyy = (1/df)*PAyy;
fprintf('\nThe variance cov(x) = %f', cov(x));
fprintf('\nThe variance sum(PAxx(0:Fn)*df) = %f\n', sum(PAxx(1:fm-1)*df));
fprintf('\nThe variance cov(y) = %f', cov(y));
fprintf('\nThe variance sum(PAyy(0:Fn)*df) = %f\n', sum(PAyy(1:fm-1)*df));
fprintf('\nHit any key to continue....');
pause;

% plot the auto spectra
subplot(221)
%axis([0 Fn 0 max(PAxx)]);
loglog(freq(2:fm),PAxx(2:fm))
title('Auto Spectra #1')
subplot(222)
%axis([0 Fn 0 max(PAyy)]);
loglog(freq(2:fm),PAyy(2:fm))
title('Auto Spectra #2')

% Calculate coherence function
I = (PAxy.*conj(PAxy))./(PAxx.*PAyy);
axis([0 Fn 0 1]);
subplot(224)
plot(freq(2:fm),I(2:fm))
title('Coherence')