Live Script edama_06
This Live Script supports Chapter 6 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_06_01: cosine wave
clear all; % clears all previously-defined variables
% I recommend that this be done at the top of every LibeScript
% standard set-up of t-axis
N=100;
Dt=1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
t = tmin + Dt*[0:N-1]';
% comoute d(t)
t0 = 10.0; % time delay of cosine
T = 50; % period of cosine
C= 2.0; % amplitude of cosine
d = C*cos( 2*pi*(t-t0)/T ); % d(t)
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -3, 3] );
plot(t, d, 'k-', 'LineWidth', 2);
plot([tmin, tmax]', [0, 0]', 'k-', 'LineWidth', 2);
xlabel('time, t');
ylabel('d(t)');
title('cosine example');
% edama_06_02: plot columns of G
% (not normalized - see eda06_13 for properly normalized version)
N=32
M=N;
% standard set-up of x-axis
Dt=1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
t = tmin + Dt*[0:N-1]';
% standard set-up of Fourier parameters
% note: N required to be an even integer
Nf=N/2+1; % number of non-negative frequencies
fmax = 1/(2*Dt); % Nyquist frequency
Df = fmax/(N/2); % frequency increment
f = Df*[0:Nf-1]'; % non-negative frequency axis
Nw=Nf; % number of angular frequencies
wmax = 2*pi*fmax; % Nyquist angukar frequency
Dw = wmax/(N/2); % angular frequency incremnt
w = Dw*[0:Nw-1]'; % non-negative angular frequency axis
% set up G
G=zeros(N,M); % first set G to zeroes
G(:,1)=1; % zero frequency column
% interior M/2-1 columns
for i = [1:M/2-1] % loop over frequencies
j = 2*i; % index of cosine column
k = j+1; % index of sine column
G(:,j)=cos(w(i+1).*t); % compute cosine column
G(:,k)=sin(w(i+1).*t); % compute sine column
end
G(:,M)=cos(w(Nw).*t); % nyquist column
% plot spectral density
figure(1)
clf;
% plot columns of G vs time
j=1;
for i = [1:5,M]
subplot(1,6,j);
j=j+1;
set(gca,'LineWidth',2);
hold on;
set(gca,'XAxisLocation','top');
axis( [-2, 2, tmin, tmax] );
axis ij;
plot( G(:,i), t, 'k-','LineWidth',2);
plot( G(:,i), t, 'k.','LineWidth',2);
xlabel(sprintf('col %d',i));
if (i==1)
ylabel('time, s');
end
end
% edama_06_03: example of aliasing
% standard set-up of t-axis
N=50;
Dt=1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
t = tmin + Dt*[0:N-1]';
% standard set-up of Fourier parameters
% note: N required to be an even integer
Nf=N/2+1; % number of non-negative frequencies
fmax = 1/(2*Dt); % Nyquist frequency
Df = fmax/(N/2); % frequency increment
f = Df*[0:Nf-1]'; % non-negative frequency axis
Nw=Nf; % number of angular frequencies
wmax = 2*pi*fmax; % Nyquist angular frequency
Dw = wmax/(N/2); % angular frequency increment
w = Dw*[0:Nw-1]'; % non-negative angular frequency axis
% low frequency oscillation
w1= 2*Dw;
w1
d1 = cos( w1*t );
% high frequency oscillation
w2= (2+N)*Dw;
w2
d2 = cos( w2*t );
% set up higher sample rate time
N=500;
Dt=0.1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
t3 = tmin + Dt*[0:N-1]';
d3 = cos( w2*t3 );
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax/2, -2, 2] );
plot(t, d1, 'k-','LineWidth',2);
plot(t, d1, 'ko','LineWidth',2);
plot([tmin, tmax]', [0, 0]', 'k-');
xlabel('time, t');
ylabel('d1(t)');
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax/2, -2, 2] );
plot(t, d2, 'k-','LineWidth',2);
plot(t, d2, 'ko','LineWidth',2);
plot(t3, d3, 'k:','LineWidth',2');
plot([tmin, tmax]', [0, 0]', 'k-','LineWidth',2);
xlabel('time, t');
ylabel('d2(t)');
% edama_06_04: another example of aliasing
% standard set-up of t-axis
N=50;
Dt=1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
t = tmin + Dt*[0:N-1]';
% standard set-up of Fourier parameters
% note: N required to be an even integer
Nf=N/2+1; % number of non-negative frequencies
fmax = 1/(2*Dt); % Nyquist frequency
Df = fmax/(N/2); % frequency increment
f = Df*[0:Nf-1]'; % non-negative frequency axis
Nw=Nf; % number of angular frequencies
wmax = 2*pi*fmax; % Nyquist angular frequency
Dw = wmax/(N/2); % angular frequency increment
w = Dw*[0:Nw-1]'; % non-negative angular frequency axis
% low frequency oscillation
w1= -2*Dw;
w1
d1 = sin( w1*t );
% high frequency oscillation
w2= (N-2)*Dw;
w2
d2 = sin( w2*t );
% set up higher sample rate time
N=500;
Dt=0.1;
tmin=0.0;
tmax = tmin+Dt*(N-1);
t3 = tmin + Dt*[0:N-1]';
d3 = cos( w2*t3 );
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax/2, -2, 2] );
plot(t, d1, 'k-','LineWidth',2);
plot(t, d1, 'ko','LineWidth',2);
plot([tmin, tmax]', [0, 0]', 'k-');
xlabel('time, t');
ylabel('d1(t)');
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax/2, -2, 2] );
plot(t, d2, 'k-','LineWidth',2);
plot(t, d2, 'ko','LineWidth',2);
plot(t3, d3, 'k:','LineWidth',2');
plot([tmin, tmax]', [0, 0]', 'k-','LineWidth',2);
xlabel('time, t');
ylabel('d2(t)');
% edama_06_05: a.s.d. of the Neuse River hydrograph
% (not normalized - see edapy_06_14 for properly normalized version)
% load data from file
D = load('../Data/neuse.txt');
Dt = D(2,1)-D(1,1);
d = D(:,2);
N=length(d);
% round off to even number of points
N=floor(N/2)*2;
M=N;
d=d(1:N);
% standard set-up of t-axis
tmin=0.0;
tmax = tmin+Dt*(N-1);
t = tmin + Dt*[0:N-1]';
% standard set-up of Fourier parameters
% note: N required to be an even integer
Nf=N/2+1; % number of non-negative frequencies
fmax = 1/(2*Dt); % Nyquist frequency
Df = fmax/(N/2); % frequency increment
f = Df*[0:Nf-1]'; % non-negative frequency axis
Nw=Nf; % number of angular frequencies
wmax = 2*pi*fmax; % Nyquist angukar frequency
Dw = wmax/(N/2); % angular frequency incremnt
w = Dw*[0:Nw-1]'; % non-negative angular frequency axis
% set up data kernel G
G=zeros(N,M); % initialize to zero
G(:,1)=ones(N,1); % zero frequency column
% interior M/2-1 columns
for i = [1:M/2-1] % loop over frequencies
j = 2*i; % index of cosine colunn
k = j+1; % index of sine column
G(:,j)=cos(w(i+1).*t); % evaluate cosine
G(:,k)=sin(w(i+1).*t); % evaluate sine
end
G(:,M)=cos(w(Nw).*t); % nyquist column
% solve least squares problem using
% analytic inverse of GTG
gtgi = 2*ones(M,1)/N;
gtgi(1)=1/N;
gtgi(M)=1/N;
mest = gtgi .* (G'*d);
% check results
dpre=G*mest;
E=(d-dpre)'*(d-dpre);
fprintf('Error in d->mest->dnew calculation: %f\n', E );
% compute power spectral density s
s=zeros(Nw,1);
s(:,1)=sqrt(mest(1)^2); % zero frequency
for i = [1:M/2-1] % interior points
j = 2*i;
k = j+1;
s(i+1) = sqrt(mest(j)^2 + mest(k)^2);
end
s(Nw) = sqrt(mest(M)^2); % Nyquist frequency
% plot spectral density
figure(1)
clf;
% plot frequency
subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
plot(f(2:Nf),s(2:Nf),'k-','LineWidth',2);
axis( [0, 0.5, 0,2000 ] );
xlabel('frequency, cycles per day');
ylabel('spectrum');
% plot log in frequency
subplot(4,1,2);
lw = 2.2;
set(gca,'LineWidth',lw,'Xscale','log','Yscale','lin'); hold on;
plot(f(2:Nf),s(2:Nf),'k-','LineWidth',2);
axis( [.001, 0.5, 0,2000 ] );
xlabel('frequency, cycles per day');
ylabel('spectrum');
% plot period
subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, 1000, 0, max(s(2:Nf)) ] );
plot(1./f(2:Nf),s(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('spectrum');
% log plot in period
subplot(4,1,4);
lw = 2.2;
set(gca,'LineWidth',lw,'Xscale','log','Yscale','lin');
hold on;
axis([10,1000,0,2000]);
hold on;
plot(1./f(2:Nf),s(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('spectrum');
% edama_06_06: a.s.d. of the Neuse River hydrograph
% using fft command
% (unnormalized spectral density - see eda06_13 for proper normalization)
% load data from file
D = load('../Data/neuse.txt');
Dt = D(2,1)-D(1,1);
d = D(:,2);
N=length(d);
% round off to even number of points
N=floor(N/2)*2;
d=d(1:N);
% generic time/frequency set up
M=N;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
% generic frequency setup
fmax=1/(2.0*Dt); % nyquist frequency
df=fmax/(N/2); % frequency sampling
Nf=N/2+1; % number of non-negagive frequencies
% f is full freqeuncy vector:
% f = [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]
f=df*[0:N/2,-N/2+1:-1]';
dw=2*pi*df; % angular frequency sampling
w=2*pi*f; % full angular frequency vector
Nw=Nf; % number of non-negative angular f's
fpos = df*[0:N/2]'; % non-negative frequencies
wpos = 2*pi*fpos; % non-negative angular frequencies
% Compute fourier coefficients
mest = fft(d);
% compute amplitude spectral density
s=abs(mest(1:Nw)); % a.s.d.
% plot spectrum
figure(1)
clf;
% plot frequency
subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
plot(f(2:Nf),s(2:Nf),'k-','LineWidth',2);
axis( [0, 0.5, 0, 2000*N/2 ] );
xlabel('frequency, cycles per day');
ylabel('spectrum');
% plot log in frequency
subplot(4,1,2);
lw = 2.2;
set(gca,'LineWidth',lw,'Xscale','log','Yscale','lin'); hold on;
plot(f(2:Nf),s(2:Nf),'k-','LineWidth',2);
axis( [.001, 0.5, 0,2000*N/2 ] );
xlabel('frequency, cycles per day');
ylabel('spectrum');
% plot period
subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, 1000, 0, max(s(2:Nf)) ] );
plot(1./f(2:Nf),s(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('spectrum');
% log plot in period
subplot(4,1,4);
lw = 2.2;
set(gca,'LineWidth',lw,'Xscale','log','Yscale','lin');
hold on;
axis([10,1000,0,2000*N/2]);
hold on;
plot(1./f(2:Nf),s(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('spectrum');
% compute original data from fourier coefficients
% Inverse FFT
dnew = ifft(mest);
E = (d-dnew)'*(d-dnew);
fprintf('Error in d->mest->dnew calculation: %f\n', E );
% edama_06_07: amplitude spectral density of Normal curves
% standard set up of time/frequency axes
N=256;
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]';
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;
% list of sqrt(variances)
s = [1.25, 2.5, 5, 10, 20, 40]';
Ns = length(s);
D = zeros(N,Ns);
Ds = zeros(Nf,Ns);
for j = [1:Ns]
t1 = tmax/2;
s2 = s(j)^2;
D(:,j) = exp( -((t-t1).^2) / (2*s2) );
tmp = fft(D(:,j));
Ds(:,j)=abs(tmp(1:Nf));
end
eda_draw(' ',D(:,1),D(:,2),D(:,3),D(:,4),D(:,5),D(:,6),' ',' ',' ',Ds(:,1),Ds(:,2),Ds(:,3),Ds(:,4),Ds(:,5),Ds(:,6));
% edama_06_08: Fourier transform of a spike
% standard time/frequency set-up
N=256;
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]';
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;
% spike function
d = zeros(N,1);
d(1)=1.0;
% fourier transform
dt = Dt*fft(d);
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [-tmax/10, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
plot(t,d,'k.','LineWidth',2);
xlabel('time, t');
ylabel('d(t)');
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [-fmax, fmax, 0, 2] );
plot(f,real(dt),'k-','LineWidth',2);
plot(f,real(dt),'k.','LineWidth',2);
xlabel('frequency, f');
ylabel('d(f)');
% edama_06_09: area under a curve using FFT
% standard time/frequency set-up
N=256;
Dt=1;
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;
% exeplary function d(t)
w0 = 2*pi*fmax/20;
d = cos(w0*t).*exp(-0.25*w0*t);
area1 = Dt*sum(d);
area1
% fourier transform; note normalization factor of Dt
dt=Dt*fft(d);
% area by fft
area2 = real(dt(1));
area2
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
title(sprintf('sum area %f fft area2 %f', area1, area2'));
xlabel('time, t');
ylabel('d(t)');
% edama_06_10: time shift using FFT
% standard time/frquency set-up
N=256;
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]';
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;
% exepplary function
w0 = 2*pi*fmax/20;
d = cos(w0*t).*exp(-0.25*w0*t);
% time shift
t0 = t(16);
ds=ifft(exp(-complex(0,1)*w*t0).*fft(d));
% Note: I recommend that you use complex(0,1)
% and not i. Although MATLAB initializes i to
% complex(0,1), it can easily be overwritten to
% some other value (e.g. a for-loop counter)
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
xlabel('time, t');
ylabel('d(t)');
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,real(ds),'k-','LineWidth',2);
xlabel('time, t');
ylabel('d_shifted(t)');
% edama_06_11: derivative using FFT
% standard time/frequency set-up
N=256;
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]';
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;
% exemplary function
t0 = tmax/2;
s = 50;
d = exp( -((t-t0).^2) / (2*s^2) );
% derivative using finite differences
dddt1 = Dt*(d(2:N)-d(1:N-1));
% derivative using fft
dddt=ifft(complex(0,1)*w.*fft(d));
% Note: I recommend that you use complex(0,1)
% and not i. Although MATLAB initializes i to
% complex(0,1), it can easily be overwritten to
% some other value (e.g. a for-loop counter)
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
xlabel('time, t');
ylabel('d(t)');
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -0.02, 0.02] );
plot(t(1:N-1),dddt1,'k-','LineWidth',2);
xlabel('time, t');
ylabel('dd/dt(t)');
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -0.02, 0.02] );
plot(t,real(dddt),'k-','LineWidth',2);
xlabel('time, t');
ylabel('dd/dt(t)');
% edama_06_12: integral using FFT
% standard time/frequency set-up
N=256;
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]';
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;
% exemplary function
t1 = tmax/4;
t2 = 3*tmax/4;
s = 25;
d = exp( -((t-t1).^2) / (2*s^2) ) - exp( -((t-t2).^2) / (2*s^2) );
% integral using Riemann sum
integral1 = Dt*cumsum(d);
% take Fourier transform and divide by iw
% except set zero frequency to zero
wr=1.0./(complex(0,1)*w(2:N));
integral=ifft(fft(d).*[0,wr']');
% Note: I recommend that you use complex(0,1)
% and not i. Although MATLAB initializes i to
% complex(0,1), it can easily be overwritten to
% some other value (e.g. a for-loop counter)
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
xlabel('time, t');
ylabel('d(t)');
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -100, 150] );
plot(t,integral1,'k-','LineWidth',2);
xlabel('time, t');
ylabel('integral');
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -100, 150] );
plot(t,real(integral),'k-','LineWidth',2);
xlabel('time, t');
ylabel('integral');
% edama_06_13: Crib Sheet 6.2
% makes exemplary time series d(t),
% amplitude spectral density s(f) pairs for
% standard set-up for time
N=128; % number of data
Dt = 0.5; % ampling interval
t = Dt*[0:N-1]'; % time vector
% standard set-up for frequency
T=N*Dt; % length of time window
fmax=1/(2*Dt); % Nyquist frequency
Df =fmax/(N/2); % frequency sampling
Nf=N/2+1; % number of non-negative frquencies
f=Df*[0:N/2,-N/2+1:-1]'; % frequncy vector
% first of two figures
figure(1);
clf;
sp=1;
% 1A. narrow function
sigmad = 1.0;
d = exp(-((t-T/2).^2)/(2*(sigmad^2))) + random('Normal',0,0.1,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, 0, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
% 1B. zero area wavelet
sigmad1 = 2.5;
sigmad2 = 1.0;
d = exp(-((t-T/2).^2)/(2*(sigmad1^2)));
d = d.*cos((t-T/2)/sigmad2);
d = d + random('Normal',0,0.1,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
% 1C. wide function
sigmad = 7.0;
d = exp(-((t-T/2).^2)/(2*(sigmad^2))) + random('Normal',0,0.1,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, 0, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
% 1D. triangle
d = [zeros(N/4,1); (0:N/4-1)'/(N/4); (N/4-1:-1:0)'/(N/4); zeros(N/4,1)];
d = d + random('Normal',0.0,0.005,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
figure(2);
clf;
sp=1;
% in all examples
% d: data time series
% dtilde: Fourier transform
% s2: power spectral density
% s: amplitude spectral density
% 2.A approximately sinusoidal
w0 = 0.25*2*pi*(0.95+0.1*[0:N-1]'/(N-1));
d = sin( w0 .* t ) + random('Normal',0,0.2,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
% 2.B beating
w1 = 0.25*2*pi;
w2 = w1/2;
d = sin( w1*t ) + 0.5*sin( w2*t ) + random('Normal',0,0.2,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
% 2.C chirp
w0 = 0.25*2*pi*(2*[0:N-1]'/(N-1));
d = sin( w0 .* t ) + random('Normal',0,0.05,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
% 2.D sawtooth
d = mod(t,N/10)/(N/10)-0.5;
d = d + random('Normal',0.0,0.005,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);
subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);
% edama_06_14: power spectral density
% comparison of variances between time series and its power spectral density
% random time series with variance sigmad^2
N=256;
Dt=0.01;
t=Dt*[0:N-1];
sigmad1=10;
sigmad1sq=sigmad1^2;
d = random('Normal',0,sigmad1,N,1);
fprintf('true variance %.2f\n', sigmad1sq);
% variance estimated from time series
sigmad2sq = (d'*d)/N;
fprintf('variance estimated from time series %.2f\n', sigmad2sq);
% Part 1: sines and cosines approach
% standard set-up of frequency axis
Nf=N/2+1;
fmax = 1/(2*Dt);
Df = fmax/(N/2);
f = Df*[0:Nf-1]';
Nw=Nf;
wmax = 2*pi*fmax;
Dw = wmax/(N/2);
w = Dw*[0:Nw-1]';
% set up G for least squares fit of sines and cosines
M=N;
G=zeros(N,M);
G(:,1)=1;
for i = [1:M/2-1]
j = 2*i;
k = j+1;
G(:,j)=cos(w(i+1).*t);
G(:,k)=sin(w(i+1).*t);
end
G(:,M)=cos(w(Nw).*t);
% solve least squares problem with analytic inv(GTG)
gtgi = 2*ones(M,1)/N;
gtgi(1)=1/N;
gtgi(M)=1/N;
mest = gtgi .* (G'*d);
sigmad3sq = 0.5*mest'*mest;
% the next line corrects for the zero-frequency and Nyquist frequncy
% (though since the correction is small, it is often omitted)
sigmad3sq = sigmad3sq + 0.5*mest(1)*mest(1) + 0.5*mest(N)*mest(N);
fprintf('variance estimated from sine/cosines %.2f\n', sigmad3sq);
% Part 2. fft() approach
% standard set-up of frequency axis
N=256;
Dt=0.1;
tmax=Dt*N;
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;
% discrete fourier transform
C=fft(d);
% variance estimated from fft coefficients
sigmad4sq = (2/(N^2)) * C(1:Nf)'*C(1:Nf);
sigmad4sq = sigmad4sq - (1/(N^2)) * (abs(C(1))^2 + abs(C(Nf))^2);
fprintf('variance estimated from fft coefficients %.2f\n', sigmad4sq);
% apprximation to integral transform
dbar=Dt*C;
% variance estimated from fourier transform
sigmad5sq = (2/(Dt^2*N^2)) * dbar(1:Nf)'*dbar(1:Nf);
% the next line corrects for the zero-frequency and Nyquist frequncy
% (though since the correction is small, it is often omitted)
sigmad5sq = sigmad5sq - (1/(Dt^2*N^2))*(abs(dbar(1))^2+abs(dbar(Nf))^2);
fprintf('variance estimated from the fourier transform %.2f\n', sigmad5sq);
% integral from 0 to fny of |dbar(f)|2
theintegral = Df * dbar(1:Nf)'*dbar(1:Nf);
sigmad5sq = (2/(Df*Dt^2*N^2)) * theintegral;
% the next line corrects for the zero-frequency and Nyquist frequncy
% (though since the correction is small, it is often omitted)
sigmad5sq = sigmad5sq -(1/(Df*Dt^2*N^2))*Df*(abs(dbar(1))^2+abs(dbar(Nf))^2);
fprintf('variance estimated from frequency integral %.2f\n', sigmad5sq);
% definition of power spectral density
thepsd = (2/tmax) * abs(dbar(1:Nf)).^2;
% integral of power spectral density
P = Df * sum(thepsd);
sigmad6sq = P;
% the next line corrects for the zero-frequency and Nyquist frequncy
% (though since the correction is small, it is often omitted)
sigmad6sq = sigmad6sq - 0.5*Df*(thepsd(1)+thepsd(Nf));
fprintf('variance estimated from power spectral density %.2f\n', sigmad6sq);
% maximum, minimum
v = [sigmad2sq,sigmad3sq,sigmad4sq,sigmad5sq,sigmad6sq]';
vmax = max(v);
vmin=min(v);
fprintf('maximum difference in estimated variances: %e\n',vmax-vmin);
% comparison of coefficients
fprintf('check agreement between sines/cosine and fft at 2nd frequency:\n');
fprintf('A2 B2 %.2f %.2f (2/N)real(C) -(2/N)imag(C) %.2f %.2f', mest(2), mest(3), (2/N)*real(C(2)), -(2/N)*imag(C(2)) );
% edama_06_15: p.s.d. of the Neuse River hydrograph
% read data from file
D = load('../Data/neuse.txt');
Dt = D(2,1)-D(1,1);
d = D(:,2);
N=length(d);
% round off to even number of points
N=floor(N/2)*2;
d=d(1:N);
% generic time/frequency set up
t=Dt*[0:N-1]';
tmax = Dt*(N-1);
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]';
Nf=N/2+1;
% compute power spectral sensity, s2
dbar=Dt*fft(d); % fourier transform
T = N*Dt; % duration of time series
s2 = (2/T) * abs(dbar(1:Nf)).^2; % p.s.d.
% plot data
figure(1)
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis([0,tmax,0,max(d)]);
plot(t,d,'k-','LineWidth',2);
xlabel('time, days');
ylabel('discharge, cfs');
% plot power spectral density (skip zero fequency)
% note that plot is only of low frequencies
figure(2);
clf;
subplot(1,3,1);
set(gca,'LineWidth',2);
pbaspect([1 2 2]);
hold on;
axis([0,0.5,0,max(s2(2:Nf))]);
plot(f(2:Nf),s2(2:Nf),'k-','LineWidth',2);
xlabel('frequency, cycles per day');
ylabel('PSD, (cfs)^2 per cycle/day');
title('linear in frequency')
subplot(1,3,2);
set(gca,'LineWidth',2);
pbaspect([1 2 2]);
hold on;
axis([0,1000,0,max(s2(2:Nf))]);
plot(1./f(2:Nf),s2(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('PSD, (cfs)^2 per cycle/day');
title('linear in period')
subplot(1,3,3);
set(gca,'LineWidth',2);
set(gca,'Xscale','log');
pbaspect([1 2 2]);
hold on;
axis([0.001,1,0,max(s2(2:Nf))]);
plot(f(2:Nf),s2(2:Nf),'k-','LineWidth',2);
xlabel('frequency, cycles per day');
ylabel('PSD, (cfs)^2 per cycle/day');
title('logarithmic in frequency')
% power in the data and the PSD (which should be the same)
% note that the mean of the data zero frequency are being excluded
v1 = std(d)^2;
v2 = Df*sum(s2(2:Nf)); % note omits 0,Nyquist corrction
fprintf('variance estimated from time series %e\n',v1);
fprintf('variance estimated from poser spectral density %e\n',v2);
% eda06_15: make a random noise dataset
% Note that he noise is saved to mynoise.txt, not noise.txt
% to avoid overwriting the standard distribution. So change
% the name if you want to overwrite noise.txt
clear all;
% standard set-up of t-axis
N=1024;
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
% make vector d of random noise
dbar=0.0; % mean
sigmad=1.0; % sqrt(variance)
d = random('Normal',dbar, sigmad, N, 1);
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0.0,tmax,-3.0*sigmad,3.0*sigmad]);
plot(t,d,'k-','LineWidth',1);
plot(t,d,'k.','LineWidth',1);
plot([0,tmax]',[0,0]','k-','LineWidth',1);
% write file
D=zeros(N,2); % concatenate (t,d) ino matrix D
D(:,1)=t;
D(:,2)=d;
dlmwrite('../Data/mynoise.txt',D,'\t');