Live Script edama_09
This Live Script supports Chapter 9 of Environmental Data Analysis with MATLAB or PYTHON, Third Edition, by WIlliam Menke
The script naming system has the form edama_CC_SS, where CC is the chapter number and SS is the script number.
% edama_09_01: Covariance & corr coeff for Atlantic Rocks data set
clear all;
% read data from file
S = load('../Data/rocks.txt');
% note use of a cell-string (cellstr) to create array of variable names
name = { 'SiO2', 'TiO2', 'Al203', 'FeO', 'MgO', 'CaO', 'Na2O', 'K2O' };
% determine size of data matrix
Ns = size(S);
N = Ns(1);
M = Ns(2);
dbar = mean(S); % mean
C = cov(S); % covariance
sigma = diag(C).^0.5; % stddev
R = corrcoef(S); % correlation coefficient
% plot results
eda_draw(abs(R));
% display results; also find most correlated data
for i = [1:M]
fprintf('%s mean=%.2f stddev=%.2f\n',char(name(i)),dbar(i),sigma(i));
end
SiO2 mean=50.39 stddev=1.30 TiO2 mean=1.36 stddev=0.38 Al203 mean=15.26 stddev=1.60 FeO mean=9.86 stddev=1.32 MgO mean=7.88 stddev=3.21 CaO mean=11.64 stddev=1.33 Na2O mean=2.47 stddev=0.47 K2O mean=0.19 stddev=0.19
fprintf('\n');
icor=1;
jcor=2;
r=0;
for i = [1:M]
for j = [i+1:M-1]
fprintf('%s vs %s cov=%.2f R=%.2f',char(name(i)),char(name(j)),C(i,j),R(i,j));
if( abs(R(i,j)) > r )
r = abs(R(i,j));
icor=i;
jcor=j;
end
end
end
SiO2 vs TiO2 cov=0.18 R=0.36SiO2 vs Al203 cov=0.10 R=0.05SiO2 vs FeO cov=0.43 R=0.25SiO2 vs MgO cov=-1.87 R=-0.45SiO2 vs CaO cov=0.10 R=0.06SiO2 vs Na2O cov=0.20 R=0.33TiO2 vs Al203 cov=0.01 R=0.02TiO2 vs FeO cov=0.25 R=0.50TiO2 vs MgO cov=-0.57 R=-0.47TiO2 vs CaO cov=-0.11 R=-0.22TiO2 vs Na2O cov=0.13 R=0.73Al203 vs FeO cov=-0.82 R=-0.39Al203 vs MgO cov=-3.52 R=-0.68Al203 vs CaO cov=1.41 R=0.66Al203 vs Na2O cov=0.29 R=0.39FeO vs MgO cov=-0.96 R=-0.23FeO vs CaO cov=-0.35 R=-0.20FeO vs Na2O cov=0.05 R=0.08MgO vs CaO cov=-2.66 R=-0.62MgO vs Na2O cov=-0.79 R=-0.53CaO vs Na2O cov=-0.06 R=-0.09
figure(2);
clf;
set(gca,'LineWidth',2);
plot(S(:,icor),S(:,jcor), 'ko', 'LineWidth', 2 );
xlabel(char(name(icor)));
ylabel(char(name(jcor)));
title(sprintf('R=%f',R(icor,jcor)));
% edama_09_02: scatter plots of lagged samples
% using the Neuse river hydrograph
% in order to investigate correlation
% read data from file
D=load('../Data/neuse.txt');
% break out time t and discharge d
t=D(:,1);
d=D(:,2);
N=length(d);
figure(1);
clf;
% plot data with 1 day lag
lag=1;
subplot(1,3,1);
set(gca,'LineWidth',2);
hold on;
axis equal;
axis( [0 25000 0 25000] );
plot( d(1:N-lag), d(1+lag:N), 'ko', 'LineWidth', 2 );
xlabel('discharge');
ylabel(sprintf('discharge lagged by %d days', lag));
% plot data with 3 day lag
subplot(1,3,2);
lag=3;
set(gca,'LineWidth',2);
hold on;
axis equal;
axis( [0 25000 0 25000] );
plot( d(1:N-lag), d(1+lag:N), 'ko', 'LineWidth', 2 );
xlabel('discharge');
ylabel(sprintf('discharge lagged by %d days', lag));
% plot data with 30 day lag
subplot(1,3,3);
lag=30;
set(gca,'LineWidth',2);
hold on;
axis equal;
axis( [0 25000 0 25000] );
plot( d(1:N-lag), d(1+lag:N), 'ko', 'LineWidth', 2 );
xlabel('discharge');
ylabel(sprintf('discharge lagged by %d days', lag));
% edama_09_03, autocorrelation of the Neuse hydrograph
% load data from time
S=load('../Data/neuse.txt');
% break out time t and discharge d data
t=S(:,1);
d=S(:,2);
N=length(d);
Dt=1;
% demean the data
d=d-mean(d);
% compute the autocorrelation
a = xcorr(d); % autocorrelation
tau = Dt*[-N+1:N-1]'; % vector of time lags
Nk = [1:N, N-1:-1:1]'; % vector of overlap-lengths
a = a ./ Nk; % correction for length of timeseries
figure(1);
clf;
% plot central part of autocorrelation
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [ -30, 30, 0, 1.1*max(a) ] );
plot( tau, a, 'k-', 'LineWidth', 2 );
xlabel('lag, days');
ylabel('autocorrelation');
% plot complete autocorrelation
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [ -365.25*10, 365.25*10, -1.1*max(a), 1.1*max(a) ] );
plot( tau, a, 'k-', 'LineWidth', 2 );
xlabel('lag, days');
ylabel('autocorrelation');
% eda09_04: align two timeseries using cross-correlation
% this example uses synthetic data
% standard seet-up of time-axis
N=1000;
Dt=0.1;
t=Dt*[1:N]';
% synthetic time series, u and v
t1 = Dt*4*N/8;
t2 = Dt*5*N/8+Dt/3;
tlagtrue = t2-t1;
w0 = 2*pi*0.5;
u = cos( w0 * (t-t1) ) .* exp( -((t-t1).^2)/2 ) + random('Normal',0,0.02,N,1);
v = cos( w0 * (t-t2) ) .* exp( -((t-t2).^2)/2 ) + random('Normal',0,0.02,N,1);
% plot original timeseries
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, u+0.1, 'k-', 'LineWidth', 2);
plot( t, v-0.1, 'k-', 'LineWidth', 2);
xlabel('time')
ylabel('u(t) and v(t)');
% cross correlation of u and v
c = xcorr(u, v); % cross-correlate
NA = length(u); % length of u
NB = length(v); % length of v
Nm = max([NA,NB]); % max of two lengths
Nc = length(c); % length of cross-correlation
tau = Dt*[-Nm+1:Nm-1]'; % time axis for xcorr
% lag to nearest sample
[cmax, icmax] = max(c); % determine location of maximum
tlag = -Dt * (icmax-Nm); % time lag
% plot lagged timeseries
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, u+0.1, 'k-', 'LineWidth', 2);
plot( t-tlag, v-0.1, 'k-', 'LineWidth', 2);
xlabel('time')
ylabel('u(t) and v(t)');
% plot cross-correlation
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [-20, 20, 1.1*min(c), 1.1*max(c)] );
plot( tau, c, 'k-', 'LineWidth', 2);
plot( -[tlag,tlag]', 1.1*[min(c),max(c)]', 'k:', 'LineWidth', 1);
xlabel('time')
ylabel('cross-correlation');
fprintf('lag: true: %.4f estimated: %.4f\n', tlagtrue, tlag );
lag: true: 12.5333 estimated: 12.5000
[ tlag2, Cmax_est ] = eda_timedelay( u, v, Dt );
fprintf('lag: true: %.4f refined: %.4f\n', tlagtrue, tlag2 );
lag: true: 12.5333 refined: 12.5295
figure(2);
clf;
% plot lagged timeseries
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, u+0.1, 'k-', 'LineWidth', 2);
plot( t-tlag2, v-0.1, 'k-', 'LineWidth', 2);
xlabel('time')
ylabel('u(t) and v(t)');
% edama_09_05; align 2 timeseries using the cross correlation
% this example uses West Point ozone data
S = load('../Data/ozone_nohead.txt');
t = S(:,1); % time in days after 1:00 AM on Aug 1, 1993
ozone = S(:,2); % ozone in ppm
radiation = S(:,3); % solar radiation in W/m2
temp = S(:,4); % temperature in deg C
% time series paramaters
N=length(t);
Dt=1/24;
% plot original timeseries
figure(1);
clf;
% plot solar radiation
subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), 0, max(radiation)] );
plot( t, radiation, 'k-', 'LineWidth', 2);
xlabel('time, days');
ylabel('solar, W/m2');
% plot ozone
subplot(4,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), 0, max(ozone)] );
plot( t, ozone, 'k-', 'LineWidth', 2);
xlabel('time, days');
ylabel('ozone, ppb');
% cross correlation of u and v
c = xcorr(radiation, ozone); % cross-correlate
NA = length(radiation); % length of u
NB = length(ozone); % length of v
Nm = max([NA,NB]); % max of two lengths
Nc = length(c); % length of cross-correlation
tau = Dt*[-Nm+1:Nm-1]'; % time axis for xcorr
% lag to nearest sample
[cmax, icmax] = max(c); % determine location of maximum
tlag = -Dt * (icmax-Nm); % time lag
fprintf('lag: %.4f\n', 24*tlag);
lag: 3.0000
% plot solar radiation
subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), 0, max(radiation)] );
plot( t, radiation, 'k-', 'LineWidth', 2);
xlabel('time, days');
ylabel('solar, W/m2');
% plot lagged ozone
subplot(4,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(N), 0, max(ozone)] );
plot( t-tlag, ozone, 'k-', 'LineWidth', 2);
xlabel('time, days');
ylabel('ozone, ppb');
figure(2);
clf;
days=5;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(days*24), 0, max(radiation)] );
plot( t, radiation, 'k-', 'LineWidth', 2);
xlabel('time, days');
ylabel('solar radiation, W/m2');
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [t(1), t(days*24), 0, max(ozone)] );
plot( t, ozone, 'k-', 'LineWidth', 2);
plot( t-tlag, ozone, 'k:', 'LineWidth', 2);
title(sprintf('%.2f hour lag', 24*tlag));
xlabel('time, days');
ylabel('ozone, ppb');
fprintf('estimated lag: %.2f hours\n', 24*tlag );
estimated lag: 3.00 hours
figure(3);
clf;
set(gca,'LineWidth',2);
hold on;
% plot central part of cross-correlation function
hours=12;
axis( [-hours, hours, 0, 1.1*max(c)] );
plot( 24*tau, c, 'k-', 'LineWidth', 2);
plot( -24*[tlag,tlag]',[0, 1.1*max(c)]','k:','LineWidth',1);
xlabel('time, hours');
ylabel('cross-correlation');
% edama_09_06: filter estimation using cross-correlation
% temperature(t) = heat(t) (convolved with) impulse_response(t)
% q(t) = h(t) * g(t)
% impulse_repsonse known (analytic formula)
% temperature is data, heat is unknown
% solution via generalized least squares
% using bicongugate gradient method
% with the prior information of smoothness
% Note: symbol q for theta, temperture
clear a e H;
global a e H;
% standard set-up of t-axis
M=1024;
N=M;
Dtdays = 0.5;
tdays = Dtdays*[0:N-1]';
Dtseconds = Dtdays*3600*24;
tseconds = Dtseconds*tdays;
% materal constants
rho = 1500; % density, kg/m3
cp = 800; % heat capacity, J / kg-K
kappa = 1.25e-6; % thermal diffusivity, m2/s
z = 1; % vertical position, m
% impulse response, and alalytic formula
g = zeros(N,1);
g(1)=0.0; % manually set first value
g(2:N) = (1/(rho*cp)) * (Dtseconds/sqrt(2*pi)) .* ...
(1./sqrt(2*kappa*tseconds(2:N))) .* ...
exp( -0.5*(z^2)./ (2*kappa*tseconds(2:N)));
% true heat production
htrue = zeros(M,1);
t1 = Dtdays*M/4;
t2 = Dtdays*7*M/16;
sd = Dtdays*M/256;
Ah=10;
htrue = Ah*exp(-((tdays(1:M)-t1).^2) / (2*sd^2) ) + 0.5*Ah*exp(-((tdays(1:M)-t2).^2) / (2*sd^2) );
% predict true temperature
tmp = conv(g, htrue);
qtrue=tmp(1:N);
% create noisy simulated observations, qobs = qtrue+noise
sigmad=0.01*max(abs(qtrue));
qobs2=qtrue+random('Normal',0,sigmad,N,1);
% solve for h given qobs and g using generalized least squares
% with a prior smoothness constraint using biconjugate gradient
% autocorrelation, a, of g. Just keep positive values
al = xcorr(g);
Na = length(al);
a = al((Na+1)/2: Na);
% set up 2nd derivative matrix H, in Hm=h
e=1*max(abs(g));
K=M-2;
L=N+K;
H=spdiags([ones(K,1),-2*ones(K,1),ones(K,1)], [0,1,2], K, M);
% set up RHS h=0
h=zeros(K,1);
% st up RHS FTf and solve
[FTf] = autofunRHS(g,qobs2,e,H,h);
mest = bicg( @autofun, FTf, 1e-10, 3*L);
bicg converged at iteration 387 to a solution with relative residual 8e-11.
hest3 = mest
hest3 =
7.5049 2.6913 -0.8450 -2.5134 -2.6614 -2.0217 -1.1931 -0.4790 0.0071 0.2673
% plot simulated observations
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), 0, 1.1*max(qobs2)] );
plot(tdays(1:N),qobs2,'k-','LineWidth',2);
xlabel('time, days');
ylabel('qobs(t)');
% plot true heat production
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(M), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:M),htrue,'k-','LineWidth',2);
xlabel('time, days');
ylabel('htrue(t)');
% plot estimated heat production
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(M), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:M),hest3,'k-','LineWidth',2);
xlabel('time, days');
ylabel('hest(t)');
% eda09_07: smoothing the Neuse River hydrograph
% load data from a file
D=load('../Data/neuse.txt');
% break out time t and discharge d data
t=D(:,1);
d=D(:,2);
Dt=t(2)-t(1);
[N, i]=size(d);
% construct 3-point smoothing filter
s1 = [1, 2, 1]';
s1 = s1/sum(s1);
L1 = length(s1);
L1h = floor(L1/2);
% filter and time shift
d1 = conv(s1, d); % filter by convolution
d1 = d1(1+L1h:N+L1h); % time shift
% construct 19-point smoothing filter
s2 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1]';
s2 = s2/sum(s2);
L2 = length(s2);
L2h = floor(L2/2);
% filter and time shift
d2 = conv(s2, d); % filter by convolutio
d2 = d2(1+L2h:N+L2h); % time shift
% plot only this time range (for clarity)
L1=1;
L2=500;
% plot
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
axis( [t(L1), t(L2), 0, max(d) ] );
hold on;
plot(t(L1:L2),d(L1:L2),'k-', 'Linewidth', 2);
xlabel('time, days');
ylabel('d-obs');
subplot(3,1,2);
set(gca,'LineWidth',2);
axis( [t(L1), t(L2), 0, max(d) ] );
hold on;
plot(t(L1:L2),d1(L1:L2),'k-', 'Linewidth', 2);
xlabel('time, days');
ylabel('3-point');
subplot(3,1,3);
set(gca,'LineWidth',2);
axis( [t(L1), t(L2), 0, max(d) ] );
hold on;
plot(t(L1:L2),d2(L1:L2),'k-', 'Linewidth', 2);
xlabel('time, days');
ylabel('21-point');
% eda09_08: amplitude spectral density of uniform smoothing filter
% standard set-up of time and frequency
N=256;
Dt=0.01;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
df=fmax/(N/2);
f=df*[0:N/2,-N/2+1:-1]';
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;
% evaluate fourier transform of filter of length L
T1=3*Dt;
s1 = sinc(w(1:Nf)*T1/(2*pi));
sa1 = abs(s1);
% case 1: short T1, plot a.s.d. of sinc function
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1] );
plot( f(1:Nf), sa1, 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('|s| for L=3');
% evaluate fourier transform of filter of length L
T2=21*Dt;
s2 = sinc(w(1:Nf)*T2/(2*pi));
sa2 = abs(s2);
% case 2: long T2, plot a.s.d. of sinc function
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1] );
plot( f(1:Nf), sa2, 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('|s| for L=21');
% edama_09_09: amplitude spectrum of Normal smoothing filter
% standard set-up of time and frequency
N=256;
Dt=0.01;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
df=fmax/(N/2);
f=df*[0:N/2,-N/2+1:-1]';
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;
% evaluate fourier transform of filter, g, of length L
L1=3;
% compute variance of filter (as if it were distribution)
g1 = ones(L1,1);
g1 = g1/sum(g1);
mg1 = sum( [1:L1]' .* g1 ) ;
vg1 = sum( (([1:L1]'-mg1).^2) .* g1);
sg1 = Dt*sqrt(vg1);
% fourier transform of gausian, using analytic formula
sigma1=1/sg1;
s1 = exp(-w(1:Nf).*w(1:Nf)/(2*sigma1^2));
sa1 = abs(s1);
fprintf('length %d samples (%.2f seconds) sigma-t samples %.2f (%.22f seconds)\n', L1, Dt*L1, sg1/Dt, sg1 );
length 3 samples (0.03 seconds) sigma-t samples 0.82 (0.0081649658092772612467 seconds)
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1] );
plot( f(1:Nf), sa1, 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('|s| for L=3');
% evaluate fourier transform of filter of length L
L2=21;
% compute variance of filter (as if it were distribution)
g2 = ones(L2,1);
g2 = g2/sum(g2);
mg2 = sum( [1:L2]' .* g2 ) ;
vg2 = sum( (([1:L2]'-mg2).^2) .* g2);
sg2 = Dt*sqrt(vg2);
fprintf('length %d samples (%.2f seconds) sigma-t samples %.2f (%.2f seconds)', L2, Dt*L2, sg2/Dt, sg2 );
length 21 samples (0.21 seconds) sigma-t samples 6.06 (0.06 seconds)
% fourier transform of gausian, using analytic formula
sigma2=1/sg2;
s2 = exp(-w(1:Nf).*w(1:Nf)/(2*sigma2^2));
sa2 = abs(s2);
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1] );
plot( f(1:Nf), sa2, 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('|s| for L=21');
% edama_09_10: complex z-plane analysis of highpass filter
% standard time/frequency setup
N=32;
Nf=N/2+1;
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]';
% a monomial filter, and its one root
F = [1, -1.1]';
zp = -F(1)/F(2); % 0=F1+F2*z so z=-F1/F2
fprintf('zero at %.2f\n',zp);
zero at 0.91
% complex z-plane
Nz=128;
zmin=-2;
zmax=2;
dz=(zmax-zmin)/(Nz-1);
zr=zmin+dz*[0:Nz-1]';
zi=zmin+dz*[0:Nz-1]';
% evaluate filter polynomial everywhere in complex z-plane
S = zeros(Nz,Nz);
for n = [1:Nz]
for m = [1:Nz]
S(m,n) = abs( F(1) + F(2)*(zr(n)+complex(0,1)*zi(m)) )^2;
end
end
% plot amplitude of polynomial
figure(1);
clf;
subplot(1,2,1);
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([zmin, zmax, zmin, zmax]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(S))-min(min(S));
if( range==0 )
range=1;
end
imagesc( [zmin, zmax], [zmin, zmax], (S-min(min(S)))/range);
hold on;
xlabel('real z');
ylabel('imag z');
title('|S|^2');
% plot Fourier transform points on the unit circle
plot(0,0,'wo','LineWidth',2);
z=exp(-complex(0,1)*w*Dt);
plot(real(z),imag(z),'k+','LineWidth',2);
% plot zero of filter
plot(real(zp),imag(zp),'w*','LineWidth',2);
% plot frequency response of filter
subplot(1,2,2);
set(gca,'LineWidth',2);
hold on;
Sz = abs( F(1)+F(2)*z ) .^2;
axis( [0, f(Nf), 0, max(Sz)] );
plot(f(1:Nf),Sz(1:Nf),'k-','LineWidth',2);
plot(f(1:Nf),Sz(1:Nf),'ko','Linewidth',2);
xlabel('frequency, Hz');
ylabel('|S|^2');
% edama_09_11: complex z-plane analysis of lowpass filter
% standard time/frequency setup
N=32;
Nf=N/2+1;
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]';
% a monomial filter, and its one root
F = [1, 1.1]';
zp = -F(1)/F(2); % 0=F1+F2*z so z=-F1/F2
fprintf('zero at %.2f\n',zp);
zero at -0.91
% complex z-plane
Nz=128;
zmin=-2;
zmax=2;
dz=(zmax-zmin)/(Nz-1);
zr=zmin+dz*[0:Nz-1]';
zi=zmin+dz*[0:Nz-1]';
% evaluate filter polynomial everywhere in complex z-plane
S = zeros(Nz,Nz);
for n = [1:Nz]
for m = [1:Nz]
S(m,n) = abs( F(1) + F(2)*(zr(n)+complex(0,1)*zi(m)) )^2;
end
end
% plot amplitude of polynomial
figure(1);
clf;
subplot(1,2,1);
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([zmin, zmax, zmin, zmax]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(S))-min(min(S));
if( range==0 )
range=1;
end
imagesc( [zmin, zmax], [zmin, zmax], (S-min(min(S)))/range);
hold on;
xlabel('real z');
ylabel('imag z');
title('|S|^2');
% plot Fourier transform points on the unit circle
plot(0,0,'wo','LineWidth',2);
z=exp(-complex(0,1)*w*Dt);
plot(real(z),imag(z),'k+','LineWidth',2);
% plot zero of filter
plot(real(zp),imag(zp),'w*','LineWidth',2);
% plot frequency response of filter
subplot(1,2,2);
set(gca,'LineWidth',2);
hold on;
Sz = abs( F(1)+F(2)*z ) .^2;
axis( [0, f(Nf), 0, max(Sz)] );
plot(f(1:Nf),Sz(1:Nf),'k-','LineWidth',2);
plot(f(1:Nf),Sz(1:Nf),'ko','Linewidth',2);
xlabel('frequency, Hz');
ylabel('|S|^2');
% edama_09_12 complex z-plane analysis of narrowband filter
% standard time/frequency setup
N=32;
Nf=N/2+1;
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]';
% a monomial filter, and its one root
F1 = [1, 0.6*(1-1.1i)]';
zp1 = -F1(1)/F1(2); % 0=F1+F2*z so z=-F1/F2
F2 = conj(F1);
zp2 = -F2(1)/F2(2); % 0=F1+F2*z so z=-F1/F2
fprintf('poles at (%.2f %.2f) and (%.2f %.2f)\n',real(zp1),imag(zp1),real(zp2),imag(zp2));
poles at (-0.75 0.83) and (-0.75 -0.83)
% complex z-plane
Nz=128;
zmin=-2;
zmax=2;
dz=(zmax-zmin)/(Nz-1);
zr=zmin+dz*[0:Nz-1]';
zi=zmin+dz*[0:Nz-1]';
% evaluate filter polynomial everywhere in complex z-plane
S = zeros(Nz,Nz);
for n = [1:Nz]
for m = [1:Nz]
v1 = 1/(F1(1) + F1(2)*(zr(n)+complex(0,1)*zi(m)));
v2 = 1/(F2(1) + F2(2)*(zr(n)+complex(0,1)*zi(m)));
S(m,n) = log(abs( v1*v2 )^2);
end
end
% plot amplitude of polynomial
figure(1);
clf;
subplot(1,2,1);
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([zmin, zmax, zmin, zmax]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(S))-min(min(S));
if( range==0 )
range=1;
end
imagesc( [zmin, zmax], [zmin, zmax], (S-min(min(S)))/range);
hold on;
xlabel('real z');
ylabel('imag z');
title('|S|^2');
% plot Fourier transform points on the unit circle
plot(0,0,'wo','LineWidth',2);
z=exp(-complex(0,1)*w*Dt);
plot(real(z),imag(z),'k+','LineWidth',2);
% plot poles of filter
plot(real(zp1),imag(zp1),'w+','LineWidth',2);
plot(real(zp2),imag(zp2),'w+','LineWidth',2);
% plot frequency response of filter
subplot(1,2,2);
set(gca,'LineWidth',2);
hold on;
Sz = abs( 1./((F1(1)+F1(2)*z) .* (F2(1)+F2(2)*z)) ) .^2;
axis( [0, f(Nf), 0, max(Sz)] );
plot(f(1:Nf),Sz(1:Nf),'k-','LineWidth',2);
plot(f(1:Nf),Sz(1:Nf),'ko','Linewidth',2);
xlabel('frequency, Hz');
ylabel('|S|^2');
% edama_09_13: complex z-plane analysis of notch filter
% standard time/frequency setup
N=32;
Nf=N/2+1;
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]';
% two roots
F1 = [1, 0.9*complex(0,1)]';
zp1 = -F1(1)/F1(2);
F2 = conj(F1);
zp2 = -F2(1)/F2(2);
% two poles
F3 = [1, 0.8*complex(0,1)]';
zp3 = -F3(1)/F3(2);
F4 = conj(F3);
zp4 = -F4(1)/F4(2);
fprintf('zeros at (%.2f %.2f) and (%.2f %.2f)\n',real(zp1),imag(zp1),real(zp2),imag(zp2));
zeros at (0.00 -1.11) and (0.00 1.11)
fprintf('zeros at (%.2f %.2f) and (%.2f %.2f)\n',real(zp3),imag(zp3),real(zp4),imag(zp4));
zeros at (0.00 -1.25) and (0.00 1.25)
% complex z-plane
Nz=128;
zmin=-2;
zmax=2;
dz=(zmax-zmin)/(Nz-1);
zr=zmin+dz*[0:Nz-1]';
zi=zmin+dz*[0:Nz-1]';
% evaluate filter polynomial everywhere in complex z-plane
S = zeros(Nz,Nz);
for n = [1:Nz]
for m = [1:Nz]
v1 = (F1(1) + F1(2)*(zr(n)+complex(0,1)*zi(m)));
v2 = (F2(1) + F2(2)*(zr(n)+complex(0,1)*zi(m)));
v3 = 1/(F3(1) + F3(2)*(zr(n)+complex(0,1)*zi(m)));
v4 = 1/(F4(1) + F4(2)*(zr(n)+complex(0,1)*zi(m)));
S(m,n) = log(abs( v1*v2*v3*v4 )^2);
end
end
% plot amplitude of polynomial
figure(1);
clf;
subplot(1,2,1);
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([zmin, zmax, zmin, zmax]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(S))-min(min(S));
if( range==0 )
range=1;
end
imagesc( [zmin, zmax], [zmin, zmax], (S-min(min(S)))/range);
hold on;
xlabel('real z');
ylabel('imag z');
title('|S|^2');
% plot Fourier transform points on the unit circle
plot(0,0,'wo','LineWidth',2);
z=exp(-complex(0,1)*w*Dt);
plot(real(z),imag(z),'k+','LineWidth',2);
% plot zeros and zeros of filter
plot(real(zp1),imag(zp1),'w*','LineWidth',2);
plot(real(zp2),imag(zp2),'w*','LineWidth',2);
plot(real(zp3),imag(zp3),'w+','LineWidth',2);
plot(real(zp4),imag(zp4),'w+','LineWidth',2);
% plot frequency response of filter
subplot(1,2,2);
set(gca,'LineWidth',2);
hold on;
Sz = abs( ((F1(1)+F1(2)*z) .* (F2(1)+F2(2)*z)) ./ ((F3(1)+F3(2)*z) .* (F4(1)+F4(2)*z)) ) .^2;
axis( [0, f(Nf), 0, max(Sz)] );
plot(f(1:Nf),Sz(1:Nf),'k-','LineWidth',2);
plot(f(1:Nf),Sz(1:Nf),'ko','Linewidth',2);
xlabel('frequency, Hz');
ylabel('|S|^2');
% edama_09_14: Chebychev bandpass filter
% limits of bandpass filter, in Hz
flow=5;
fhigh=10;
% standard time/frequency setup
N=1024;
Nf=N/2+1;
Dt=0.01;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]';
% input timeseries is a spike
din=zeros(N,1);
din(100)=1;
% filter
[dout, u, v] = chebyshevfilt( din, Dt, flow, fhigh);
% fft
dinspec = abs(fft(din));
doutspec = abs(fft(dout));
% plot input
figure(1);
clf;
subplot(2,2,1);
set(gca,'LineWidth',2);
hold on;
axis( [0.75, 1.25, -1, 1 ] );
plot( t, din, 'k-', 'LineWidth', 2);
xlabel('time, s');
ylabel('input');
% plot input spectrum
subplot(2,2,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1.1*max(dinspec) ] );
plot( f(1:Nf), dinspec(1:Nf), 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('input spectrum');
% plot output
subplot(2,2,3);
set(gca,'LineWidth',2);
hold on;
axis( [0.75, 1.25, -1, 1 ] );
plot( t, dout, 'k-', 'LineWidth', 2);
xlabel('time, s');
ylabel('output');
% plot output spectrum
subplot(2,2,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1.1*max(dinspec) ] );
plot( f(1:Nf), doutspec(1:Nf), 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('output spectrum');
% roots of filters
ru = roots(flipud(u));
Nu = length(u);
Nru = length(ru);
rv = roots(flipud(v));
Nv=length(v);
Nrv = length(rv);
% complex z-plane setup
Nz=100;
zmin=-2;
zmax=2;
dz=(zmax-zmin)/(Nz-1);
zr=zmin+dz*[0:Nz-1]';
zi=zmin+dz*[0:Nz-1]';
% evaluate filter polynomials everywhere in complex z-plane
S = zeros(Nz,Nz);
for n = [1:Nz]
for m = [1:Nz]
z = zr(n)+complex(0,1)*zi(m);
uu=0.0;
zpow=1;
for k = [1:Nu]
uu = uu + u(k)*zpow;
zpow=zpow*z;
end
vv=0.0;
zpow=1;
for k = [1:Nv]
vv = vv + v(k)*zpow;
zpow=zpow*z;
end
SS = uu/vv;
S(m,n)=abs(SS);
end
end
% figure showing position of roots in the z-plane
% plot amplitude of polynomial
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
zmin=-2;
zmax=2;
axis([zmin, zmax, zmin, zmax]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(S))-min(min(S));
if( range==0 )
range=1;
end
imagesc( [zmin, zmax], [zmin, zmax], (S-min(min(S)))/range);
xlabel('real z');
ylabel('imag z');
% plot origin
plot(0,0,'ok','LineWidth',2);
z=exp(-complex(0,1)*w*Dt);
plot(real(z),imag(z),'k-','LineWidth',2);
% plot zeros and poles
for k = [1:Nru]
plot(real(ru(k)),imag(ru(k)),'k*','LineWidth',2);
end
for k = [1:Nrv]
plot(real(rv(k)),imag(rv(k)),'k+','LineWidth',2);
end
% edama_09_15 plot Reynolds Channel Water Quality dataset
% columns are days, precip, air_temp, water_temp, salinity,
% turbidity, chlorophyl
% load data from file
D = load('../Data/reynolds_interpolated.txt');
[N, Ncols] = size(D);
Dt=1;
Df=D;
flow = 1/30;
fhigh = 1/5;
% beak out columns of data
t = Df(:,1); % time, days
p = Df(:,2); % precipitation, inches
at = Df(:,3); % air temperature, deg C
wt = Df(:,4); % water temperature, deg C
s = Df(:,5); % salinity, practical salinity units
tb = Df(:,6); % turbidity, (no units given)
c = Df(:,7); % chlorophyl, micrograms per liter
% plot limits
itmin=1;
itmax=N;
tmin=t(itmin);
tmax=t(itmax);
figure(2);
clf;
subplot(6,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(p(itmin:itmax)), max(p(itmin:itmax))] );
plot( t, p, 'k-', 'LineWidth', 2 );
ylabel('precipitation');
subplot(6,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(at(itmin:itmax)), max(at(itmin:itmax))] );
plot( t, at, 'k-', 'LineWidth', 2 );
ylabel('T-air');
subplot(6,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(wt(itmin:itmax)), max(wt(itmin:itmax))] );
plot( t, wt, 'k-', 'LineWidth', 2 );
ylabel('T-water');
subplot(6,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(s(itmin:itmax)), max(s(itmin:itmax))] );
plot( t, s, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('salinity');
subplot(6,1,5);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(tb(itmin:itmax)), max(tb(itmin:itmax))] );
plot( t, tb, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('turbidity');
subplot(6,1,6);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(c(itmin:itmax)), max(c(itmin:itmax))] );
plot( t, c, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('chlorophyll');
% edama_09_16: bandpass filter Reynolds Channel Water Quality data
% load data from file
D = load('../Data/reynolds_interpolated.txt');
[N, Ncols] = size(D);
Dt=1;
Df=D;
% figure 1 (low pass filter)
fcenter = 1/365.25;
flow = 0.75*fcenter;
fhigh =1.25*fcenter;
% bandpass filter
for i = [2:Ncols]
[z, u, v] = chebyshevfilt(D(:,i), Dt, flow, fhigh);
Df(:,i) = z;
end
% break-out columns of data
t = Df(:,1); % time, days
p = Df(:,2); % precipitation, inches
at = Df(:,3); % air temperature, deg C
wt = Df(:,4); % water temperature, deg C
s = Df(:,5); % salinity, practical salinity units
tb = Df(:,6); % turbidity,
c = Df(:,7); % chlorophyll, micrograms per liter
% plot limits in samples
itmin=1;
itmax=N;
tmin=t(itmin);
tmax=t(itmax);
figure(1);
clf;
subplot(6,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(p(itmin:itmax)), max(p(itmin:itmax))] );
plot( t, p, 'k-', 'LineWidth', 2 );
ylabel('precip');
subplot(6,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(at(itmin:itmax)), max(at(itmin:itmax))] );
plot( t, at, 'k-', 'LineWidth', 2 );
ylabel('T-air');
subplot(6,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(wt(itmin:itmax)), max(wt(itmin:itmax))] );
plot( t, wt, 'k-', 'LineWidth', 2 );
ylabel('T-water');
subplot(6,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(s(itmin:itmax)), max(s(itmin:itmax))] );
plot( t, s, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('salinity');
subplot(6,1,5);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(tb(itmin:itmax)), max(tb(itmin:itmax))] );
plot( t, tb, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('turbidity');
subplot(6,1,6);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(c(itmin:itmax)), max(c(itmin:itmax))] );
plot( t, c, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('chlorophyl');
% figure 2 (high pass filter)
figure(2);
clf;
fcenter = 1/5;
flow = 0.75*fcenter;
fhigh =1.25*fcenter;
% bandpass filter
for i = [2:Ncols]
[z, u, v] = chebyshevfilt(D(:,i), Dt, flow, fhigh);
Df(:,i) = z;
end
% break out columns of data
t = Df(:,1); % time, days
p = Df(:,2); % precipitation, inches
at = Df(:,3); % air temperature, deg C
wt = Df(:,4); % water temperature, deg C
s = Df(:,5); % salinity, practical salinity units
tb = Df(:,6); % turbidity,
c = Df(:,7); % chlorophyll, micrograms per liter
% plot limits in samples
itmin=400;
itmax=600;
tmin=t(itmin);
tmax=t(itmax);
subplot(6,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(p(itmin:itmax)), max(p(itmin:itmax))] );
plot( t, p, 'k-', 'LineWidth', 2 );
ylabel('precip');
subplot(6,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(at(itmin:itmax)), max(at(itmin:itmax))] );
plot( t, at, 'k-', 'LineWidth', 2 );
ylabel('T-air');
subplot(6,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(wt(itmin:itmax)), max(wt(itmin:itmax))] );
plot( t, wt, 'k-', 'LineWidth', 2 );
ylabel('T-water');
subplot(6,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(s(itmin:itmax)), max(s(itmin:itmax))] );
plot( t, s, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('salinity');
subplot(6,1,5);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(tb(itmin:itmax)), max(tb(itmin:itmax))] );
plot( t, tb, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('turbidity');
subplot(6,1,6);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, min(c(itmin:itmax)), max(c(itmin:itmax))] );
plot( t, c, 'k-', 'LineWidth', 2 );
xlabel('time, days');
ylabel('chlorophyl');
% 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));
pass-band: 0.187500 to 0.312500
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
fprintf('Calculated directly from the band-passed time series:\n');
Calculated directly from the band-passed time series:
fprintf(' Power in u %f.2\n',Pu);
Power in u 109.512161.2
fprintf(' Power in v %f.2\n',Pv);
Power in v 84.136364.2
fprintf(' Zero-lag cross correlation %.2f\n',x0);
Zero-lag cross correlation 39.74
fprintf(' Curly C squared %.2f\n', x0^2/(Pu*Pv));
Curly C squared 0.17
fprintf('\n');
% 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
fprintf('Calculated from the cross spectral density\n');
Calculated from the cross spectral density
fprintf(' Pover in u %.2f\n',usua);
Pover in u 109.51
fprintf(' Pover in v %.2f\n',vsva);
Pover in v 84.14
fprintf(' Zero-lag cross correlation %.2f\n',usvar);
Zero-lag cross correlation 39.74
fprintf(' Curly C squared %.2f\n', curlyc2);
Curly C squared 0.17
fprintf(' Coherence %,2f\n', coherence);
Coherence
% edama_09_18: coherence of Reynolds Channel dataset
% set to 1 to add results from mscohere() function
% (which are plotted in red). Set to 0 to omit them.
DOMSCOHERE = 1;
% this is an example of a cell-string (cellstr) variable, that is,
% an array of character strings its used here to hold the name of the variables
names = {'time' 'precipitation' 'air-temp' 'water-temp' 'salinity' 'turbidity' 'chlorophyll'};
% load data from file
% columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl
D = load('../Data/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.
% figure our various lengths needed by mscohere()
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
% edama_09_19: effect of boxcar window on power spectral density
% standard time/frequency definitions
N=512;
Nf=N/2+1;
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]';
% boxcar window function, W
W = zeros(N,1);
W(1:N/8)=1;
Wtmax = tmax/8;
% exemplary timeseries is a sinusoid
w0 = 40*Dw;
d = sin( w0 * t );
% windowed timeseries
Wd = W .* d;
% fourier transform
dfft=fft(d);
Wfft=fft(W);
Wdfft=fft(Wd);
% power spectral density
ds=(Dt^2/tmax)*(abs(dfft).^2);
Ws=(Dt^2/Wtmax)*(abs(Wfft).^2);
Wds=(Dt^2/Wtmax)*(abs(Wdfft).^2);
% time domain plot
figure(1);
clf;
subplot(3,2,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
xlabel('time t, s');
ylabel('d(t)');
subplot(3,2,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,W,'k-','LineWidth',2);
xlabel('time t, s');
ylabel('W(t)');
subplot(3,2,5);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,Wd,'k-','LineWidth',2);
xlabel('time t, s');
ylabel('W(t)*d(t)');
% frequency domain plot
subplot(3,2,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(ds))] );
plot(f(1:Nf),sqrt(ds(1:Nf)),'k-','LineWidth',2);
xlabel('frequency f, Hz');
ylabel('ASD of d');
subplot(3,2,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(Ws))] );
plot(f(1:Nf),sqrt(Ws(1:Nf)),'k-','LineWidth',2);
xlabel('frequency f, Hz');
ylabel('ASD of W');
subplot(3,2,6);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(Wds))] );
plot(f(1:Nf),sqrt(Wds(1:Nf)),'k-','LineWidth',2);
xlabel('frequency f, Hz');
ylabel('ASD of Wd)');
% eda09_20: effect of Hamming window on power spectral density
% standard time/frequency definitions
N=512;
Nf=N/2+1;
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]';
% Hamming window function, W
W = zeros(N,1);
Nw = N/8;
Wtmax = tmax/8;
W(1:Nw)=0.54-0.46*cos(2*pi*[0:Nw-1]'/(Nw-1));
% timeseriee
w0 = 40*Dw;
d = sin( w0 * t );
% windowed timeseries
Wd = W .* d;
% fourier transform
dfft=fft(d);
Wfft=fft(W);
Wdfft=fft(Wd);
% power spectral density
ds=(Dt^2/tmax)*(abs(dfft).^2);
Ws=(Dt^2/Wtmax)*(abs(Wfft).^2);
Wds=(Dt^2/Wtmax)*(abs(Wdfft).^2);
% time domain plot
figure(1);
clf;
subplot(3,2,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
xlabel('time t, s');
ylabel('d(t)');
subplot(3,2,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,W,'k-','LineWidth',2);
xlabel('time t, s');
ylabel('W(t)');
subplot(3,2,5);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,Wd,'k-','LineWidth',2);
xlabel('time t, s');
ylabel('W(t)*d(t)');
% frequency domain plot
subplot(3,2,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(ds))] );
plot(f(1:Nf),sqrt(ds(1:Nf)),'k-','LineWidth',2);
xlabel('frequency f, Hz');
ylabel('ASD of d');
subplot(3,2,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(Ws))] );
plot(f(1:Nf),sqrt(Ws(1:Nf)),'k-','LineWidth',2);
xlabel('frequency f, Hz');
ylabel('ASD of W');
subplot(3,2,6);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(Wds))] );
plot(f(1:Nf),sqrt(Wds(1:Nf)),'k-','LineWidth',2);
xlabel('frequency f, Hz');
ylabel('ASD of Wd)');
% eda09_21
% MultiTaper analysis
% standard time/frequency set up
N=64;
Nf=N/2+1;
Dt=1;
t=Dt*[0:N-1]';
tmax = Dt*(N-1);
fmax=1/(2*Dt);
Df=fmax/(N/2);
Dw=2*pi*Df;
f=Df*[0:N/2,-N/2+1:-1]';
% choose bandwidth to be +/-3 Dw
f0 = 2*Df;
w0 = 2*pi*f0;
% matrix, M, of sinc functions
M = zeros(N,N);
for m = [1:N]
for n = [1:N]
M(m,n) = 2*w0*sinc( w0*(m-n)*Dt / pi);
end
end
% solve algebric eigenvalue problem
% eig sorts in ascending order
[V,D] = eig(M);
lambda=diag(D);
fprintf('6 largest eigenvalues:\n');
6 largest eigenvalues:
for m = [N-5:N]
fprintf(' %d: %.2f\n', m, lambda(m) );
end
59: 0.27 60: 1.72 61: 4.54 62: 6.03 63: 6.27 64: 6.28
% four largest eigenvalues
w1 = V(:,N);
w2 = V(:,N-1);
w3 = V(:,N-2);
w4 = V(:,N-3);
% power spectral density of window functions
s1 = ((Dt^2)/tmax) * abs(fft(w1)).^2;
s2 = ((Dt^2)/tmax) * abs(fft(w2)).^2;
s3 = ((Dt^2)/tmax) * abs(fft(w3)).^2;
s4 = ((Dt^2)/tmax) * abs(fft(w4)).^2;
% plot window functions
figure(1);
clf;
winmin = -0.25;
winmax = -winmin;
subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w1, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w1(t)');
subplot(4,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w2, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w2(t)');
subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w3, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w3(t)');
subplot(4,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w4, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w4(t)');
% plot amplitude spectral density of window functions
figure(2);
clf;
subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s1))] );
plot( f(1:Nf), sqrt(s1(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s1))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w1(f)|');
subplot(4,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s2))] );
plot( f(1:Nf), sqrt(s2(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s2))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w2(f)|');
subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s3))] );
plot( f(1:Nf), sqrt(s3(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s3))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w3(f)|');
subplot(4,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s4))] );
plot( f(1:Nf), sqrt(s4(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s4))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w1(f)|');
% new, longer time interval, with same sampling interval
N2=512;
Nf2=N2/2+1;
t2=Dt*[0:N2-1]';
tmax2 = Dt*(N2-1);
fmax=1/(2*Dt);
Df2=fmax/(N2/2);
Dw2=2*pi*Df2;
f2=Df2*[0:N2/2,-N2/2+1:-1]';
% create time series and window it
ww= 2*pi*fmax/3.3;
d=cos(ww*t2);
dw0 = [[1:N], 0*[1:N2-N] ]'.*d;
dw1 = [w1', 0*[1:N2-N] ]'.*d;
dw2 = [w2', 0*[1:N2-N] ]'.*d;
dw3 = [w3', 0*[1:N2-N] ]'.*d;
dw4 = [w4', 0*[1:N2-N] ]'.*d;
% plot windowed time series
figure(3);
clf;
subplot(5,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(d), max(d) ]);
plot( t2, d, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('d(t)');
subplot(5,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(d), max(d) ]);
plot( t2, dw0, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('d(t)');
subplot(5,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(dw1), max(dw1) ]);
plot( t2, dw1, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w1(t)d(t)');
subplot(5,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(dw2), max(dw2) ]);
plot( t2, dw2, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w2(t)d(t)');
subplot(5,1,5);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(dw3), max(dw3) ]);
plot( t2, dw3, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w3(t)d(t)');
% power spectral density of time series
ds = ((Dt^2)/tmax) * abs(fft(d)).^2;
% power spectral density of windowed time series
dw0s = ((Dt^2)/tmax) * abs(fft(dw0)).^2;
dw1s = ((Dt^2)/tmax) * abs(fft(dw1)).^2;
dw2s = ((Dt^2)/tmax) * abs(fft(dw2)).^2;
dw3s = ((Dt^2)/tmax) * abs(fft(dw3)).^2;
dw4s = ((Dt^2)/tmax) * abs(fft(dw4)).^2;
% average power spectral densities
dws = (1/3)*( dw1s + dw2s + dw3s );
% plot power spectral density of timeseries
figure(4);
clf;
subplot(6,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(ds))] );
plot( f2(1:Nf2), sqrt(ds(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d(f)|');
subplot(6,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw0s))] );
plot( f2(1:Nf2), sqrt(dw0s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d0(f)|');
subplot(6,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw1s))] );
plot( f2(1:Nf2), sqrt(dw1s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d1(f)|');
subplot(6,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw2s))] );
plot( f2(1:Nf2), sqrt(dw2s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d2(f)|');
subplot(6,1,5);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw3s))] );
plot( f2(1:Nf2), sqrt(dw3s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d3(f)|');
subplot(6,1,6);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dws))] );
plot( f2(1:Nf2), sqrt(dws(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d(f)|avg');
u = [1, 0, 0, 0, 0]';
v = [1, 0, 0]';
% cross correlation
c = xcorr(u, v); % cross-correlate
NA = length(u); % length of u
NB = length(v); % length of v
Nm = max([NA,NB]); % max of two lengths
Nc = length(c); % length of cross-correlation
tau = Dt*[-Nm+1:Nm-1]'; % time axis for xcorr
c'
ans =
0 0 0 0 1 0 0 0 0
tau'
ans =
-0.4000 -0.3000 -0.2000 -0.1000 0 0.1000 0.2000 0.3000 0.4000
NA
NA = 5
NB
NB = 3
Nc
Nc = 9
% lag to nearest sample
[cmax, icmax] = max(c); % determine location of maximum
tlag = -Dt * (icmax-Nm); % time lag
tlag
tlag = 0