% eda09_17
% cohernce calculation (single frequency band only)
% as a check of the formula for C

% time definitions
N=256;
Dt=1;
t=Dt*[0:N-1]';

% random time series
rts1=random('Normal', 0, 1, 256, 1 );
rts2=random('Normal', 0, 1, 256, 1 );
rts3=random('Normal', 0, 1, 256, 1 );

% create partially correlated time series
uorig = rts1 + rts2;
vorig = rts1 + rts3;

% frequency definition
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]'; 
Nf=N/2+1;

% frequency and bandwith for calculation
f0 = fmax/2;
lowside = 0.75;
highside =1.25;
flow = f0*lowside;
fhigh= f0*highside;
ilow = floor(flow/Df)+1;
ihigh = floor(fhigh/Df)+1;
disp(sprintf('pass-band: %f to %f', flow, fhigh));
disp(' ');

% box car filter matching the parameters above
bp = zeros(N,1);
bp(ilow:ihigh)=1.0;
bp(N-ihigh+2:N-ilow+2)=1.0;

% bandpass filter the data
ufilt = ifft( bp .* fft(uorig) );
vfilt = ifft( bp .* fft(vorig) );

% calculation of power and zero-lag cross-correlation
Pu = sum( ufilt.^2 );
Pv = sum( vfilt.^2 );
x = xcorr( ufilt, vfilt );
x0 = x(N);

% plot time series
figure(1);
clf;
set(gca,'LineWidth',2);
axis( [t(1), t(N), -10, 10 ]);
hold on;
plot( t, ufilt+3, 'k-', 'LineWidth', 2 );
plot( t, vfilt-3, 'k-', 'LineWidth', 2 );
xlabel('time');
ylabel('u(t) and v(t)');

% display results
disp('Calculated directly from the band-passed time series:');
disp(sprintf('    Power in u %f',Pu));
disp(sprintf('    Power in v %f',Pv));
disp(sprintf('    Zero-lag cross correlation %f',x0));
disp(sprintf('    Curly C squared %f', x0^2/(Pu*Pv)));
disp(' ');

% Fourier transform
u = fft( uorig ); % fft
v = fft( vorig );
u = u(1:Nf); % delete negative frequencies
v = v(1:Nf);

% spectral densities
usv = conj(u) .* v; % cross spectral density of u and v
usu = conj(u) .* u; % power spectral density of u
vsv = conj(v) .* v; % power spectral density of v

% band-averaged quantities
usva = sum( (usv(ilow:ihigh)) ) / (Nf-1);
usvar = sum( real( (usv(ilow:ihigh) ) ) )/ (Nf-1);
usua = sum( usu(ilow:ihigh) )/ (Nf-1);
vsva = sum( vsv(ilow:ihigh) )/ (Nf-1);
curlyc2 = (conj(usvar)*usvar) / (usua * vsva);
coherence = (conj(usva)*usva) / (usua * vsva);

% display quantities
disp('Calculated from the cross spectral density');
disp(sprintf('    Pover in u %f',usua));
disp(sprintf('    Pover in v %f',vsva));
disp(sprintf('    Zero-lag cross correlation %f',usvar));
disp(sprintf('    Curly C squared %f', curlyc2));
disp(sprintf('    Coherence %f', coherence));
