% eda09_18
% read Reynolds Channel Water Quality dataset
% columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl
% compute coherence and plot it

% this is an example of a cellstr variable, that is, an array of character strings
% its used here to hold the name of the variables

clear all;

% set to 1 to add results from mscohere() function
% (which are plotted in red).  Set to 0 to omit them.
DOMSCOHERE = 1;

names = {'time' 'precipitation' 'air-temp' 'water-temp' 'salinity' 'turbidity' 'chlorophyll'};

% load data
D = load('reynolds_interpolated.txt');
[N, Ncols] = size(D);
N=2*floor(N/2);
Nf = N/2+1;
Dt=1;
fny = 1/(2*Dt);
Df = fny/(N/2);
f = Df*[0:Nf-1];

% bandwidth factors
lowside = 0.75;
highside =1.25;

ifig=1;
figure(1);
clf;

% loop over all pairs of columns, computing coherence for each
for icol = [2:Ncols]
for jcol = [icol+1:Ncols]

c = zeros(Nf,1); % coherence

% compute cross spectral density and power spectral density
% no need to normalize, since all normalizations cancel
u = fft( D(1:N,icol) ); % fft
v = fft( D(1:N,jcol) );
u = u(1:Nf); % delete negative frequencies
v = v(1:Nf);
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

% average over band
usva = zeros(Nf,1);
usua = zeros(Nf,1);
vsva = zeros(Nf,1); 
for i = 1:Nf
    fi = Df*(i-1); % center frequency
    flow = lowside*fi;
    fhigh = highside*fi;
    ilow = floor(flow/Df)+1;
    ihigh = floor(fhigh/Df)+1;
    if( ihigh==ilow)
        ihigh=ilow+1;
    end
    if (ilow < 1)
        ilow=1;
    elseif (ilow > Nf)
        ilow=Nf;
    end
    if (ihigh < 1)
        ihigh=1;
    elseif (ihigh > Nf)
        ihigh=Nf;
    end
    usva = mean( usv(ilow:ihigh) );
    usua = mean( usu(ilow:ihigh) );
    vsva = mean( vsv(ilow:ihigh) );
    c(i) = (conj(usva)*usva) / (usua * vsva) ;
end

% plot coherence (black)
subplot(3,5,ifig);
set(gca,'LineWidth',2);
hold on;
axis( [0, 1/5, 0, 1.1] );
plot( f, c, 'k-', 'LineWidth', 2 );
xlabel('frequency, cycles per day');
ylabel('coherence');
title(sprintf('%s and %s',char(names(icol)),char(names(jcol))));

if( DOMSCOHERE == 1)
% result using mscohere() (plotted in red)
% the black and red traced will not agree
% perfectly, because the bandwidth parameters
% cannot be made to match, expecially at low
% frequencies.
N2 = 2;
while (N2 < N ) % find smallest power of 2 greater than N
    N2 = 2*N2;
end
N3 = floor(N2/16);
N4 = 3;
[Cmatlab, fmatlab] = mscohere( D(1:N,icol), D(1:N,jcol), hamming(N3), N3-N4, N2 );
plot( (fny/pi)*fmatlab, Cmatlab, 'r-', 'LineWidth', 2 );
end
ifig=ifig+1;

end
end

figure(6);


