Live Script edama_07
This Live Script supports Chapter 7 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_07_01: least squares estimation of filter input
% example using impulse response of a heat-producing layer
% with solution by damped least squares
% Note: symbol q for theta, temperture
% This is the non-sparse implementation.
% time setup for g (length M)
M=1024;
N=M;
Dtdays = 0.5;
tdays = Dtdays*[0:N-1]';
Dtseconds = Dtdays*3600*24;
tseconds = Dtseconds*tdays;
max(tdays)
% materal constants
rho = 1500; % kg/m3
cp = 800; % J / kg-K
kappa = 1.25e-6; % m2/s
z = 1; % m
% impulse response, which is a complicated analytic formula
g = zeros(M,1);
g(1)=0.0; % manually set first value
g(2:M) = (1/(rho*cp)) * (Dtseconds/sqrt(2*pi)) .* ...
(1./sqrt(2*kappa*tseconds(2:M))) .* ...
exp( -0.5*(z^2)./ (2*kappa*tseconds(2:M)));
% true model parameters, heat production by convolving g and htrue
htrue = zeros(N,1);
t1 = Dtdays*N/4;
t2 = Dtdays*7*N/16;
sd = Dtdays*N/256;
Ah=10;
htrue = Ah*exp(-((tdays(1:N)-t1).^2) / (2*sd^2) ) + 0.5*Ah*exp(-((tdays(1:N)-t2).^2) / (2*sd^2) );
% predict true temperature
tmp = conv(g, htrue);
qtrue=tmp(1:N); % keep only first N values
% plot impulse response
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(M), 0, 1.1*max(g)] );
plot(tdays(1:M),g,'k-','LineWidth',2);
xlabel('time, days');
ylabel('g(t)');
% plot true heat production
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),htrue,'k-','LineWidth',2);
xlabel('time, days');
ylabel('htrue(t)');
% plot true temperature
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), 0, 1.1*max(qtrue)] );
plot(tdays(1:N),qtrue,'k-','LineWidth',2);
xlabel('time, days');
ylabel('qtrue1(t)');
% Toeplitz matrix version of impulse response
G = toeplitz([g', zeros(N-M,1)']', [g(1), zeros(1,M-1)] );
% temperature = G times heat_production
qtrue2 = G*htrue;
% create simulated observations, qobs = qtrue+noise
sigmad=0.001*max(abs(qtrue)); % noise is 0.1% amplitude of signal
qobs=qtrue+random('Normal',0,sigmad,N,1);
% plot simulated observations
figure(2);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), 0, 1.1*max(qobs)] );
plot(tdays(1:N),qobs,'k-','LineWidth',2);
xlabel('time, days');
ylabel('qobs(t)');
% solve for h given qobs and g using damped least squares
GTG=G'*G; % G-transpose times G
GTGmax = max(max(abs(GTG))); % maximum value of G
e2=1e-4*GTGmax; % damping is a fraction of GTG
hest1=(GTG+e2*eye(M,M))\(G'*qobs);
% plot true heat production
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),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(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),hest1,'k-','LineWidth',2);
xlabel('time, days');
ylabel('hest(t)');
% edama_07_02: estimation of filter input, with smoothing
% using impulse repsonse of a heat-producing layer
% with solution via generalized least squares
% with the prior information of smoothness
% Note: symbol q for theta, temperture
% This is the non-sparse implementation.
% generic time setup, except make two time arrays.
% one in seconds for calulations, one in days for plotting
M=1024;
N=M;
Dtdays = 0.5;
tdays = Dtdays*[0:N-1]';
Dtseconds = Dtdays*3600*24;
tseconds = Dtseconds*tdays;
% materal constants
rho = 1500; % kg/m3
cp = 800; % J / kg-K
kappa = 1.25e-6; % m2/s
z = 1; % m
% impulse response, which is a analysic formula that's fairly complicated
g = zeros(M,1);
g(1)=0.0; % manually set first value
g(2:M) = (1/(rho*cp)) * (Dtseconds/sqrt(2*pi)) .* ...
(1./sqrt(2*kappa*tseconds(2:M))) .* ...
exp( -0.5*(z^2)./ (2*kappa*tseconds(2:M)));
% true model parameters, heat production
htrue = zeros(N,1);
t1 = Dtdays*N/4; % time of heat pulse 1
t2 = Dtdays*7*N/16; % time of heat pulse 2
sd = Dtdays*N/256; % width of pulses
Ah = 10; % amplitde of pulses
% pulses have a gaussian shape
htrue = Ah*exp(-((tdays(1:N)-t1).^2) / (2*sd^2) ) + 0.5*Ah*exp(-((tdays(1:N)-t2).^2) / (2*sd^2) );
% predict true temperature
tmp = conv( g, htrue);
qtrue=tmp(1:N);
% matrix version of impulse response
G = toeplitz([g', zeros(N-M,1)']', [g(1), zeros(1,M-1)] );
qtrue2 = G*htrue;
% create noisy simulated observations, qobs = qtrue+noise
sigmad=0.01*max(abs(qtrue));
qobs2=qtrue+random('Normal',0,sigmad,N,1);
% 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)');
% solve for h given qobs and g using damped least squared
GTG=G'*G;
GTGmax = max(max(abs(GTG))); % maximum value of GTG
e2=1e-4*GTGmax; % damping is fraction of GTG
% damped least squares solution
hest2=(GTG+e2*eye(M,M))\(G'*qobs2);
% plot true heat production
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),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(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),hest2,'k-','LineWidth',2);
xlabel('time, days');
ylabel('hest(t)');
% solve for h given qobs and g using genralized least squares
% with a prior smoothness constraint
e=10*max(max(abs(G))); % damping = sigma_h/sigma_d;
F=zeros(N+M-2,M);
f=zeros(N+M-2,1);
F(1:N,1:M)=G;
f(1:N)=qobs2;
H=toeplitz( [1, zeros(1,M-3)]', ...
[1, -2, 1, zeros(1,M-3) ] );
F(N+1:N+M-2,:)=e*H;
f(N+1:N+M-2)=0;
hest3=(F'*F)\(F'*f);
% plot simulated observations
figure(2);
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(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),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(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),hest3,'k-','LineWidth',2);
xlabel('time, days');
ylabel('hest(t)');
% edama_07_03: sparse, least sq estimation of filter input
% using impulse repsonse of a heat-producing layer
% with solution via generalized least squares
% with the prior information of smoothness
% Note: symbol q for theta, temperture
% This is the sparse implementation.
clear all;
clear e g H;
global e g H;
% time setup for g (length M)
M=1024;
N=M;
Dtdays = 0.5;
tdays = Dtdays*[0:N-1]';
Dtseconds = Dtdays*3600*24;
tseconds = Dtseconds*tdays;
% materal constants
rho = 1500; % density in kg/m3
cp = 800; % heat capacity in J / kg-K
kappa = 1.25e-6; % thermal diffusivity in m2/s
z = 1; % depth of observation in m
% impulse response, an analysic formula that's fairly complicated
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)));
% note that g is a global variable
% true heat production, sum of two pulses
htrue = zeros(M,1);
t1 = Dtdays*M/4; % time of pulse 1
t2 = Dtdays*7*M/16; % time of pulse 2
sd = Dtdays*M/256; % width of pulses
Ah=10; % amplitude of pulses
% sum of two gaussian pulses
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 by colvolving htrue with impulse response
tmp = conv(g, htrue);
qtrue=tmp(1:N);
% create noisy simulated observations, qobs = qtrue+noise
sigmad=0.01*max(abs(qtrue)); % fraction of signal
qobs=qtrue+random('Normal',0,sigmad,N,1);
% for comparison purposes, solve for h given qobs and G using damped least squares
% build Toeplitz matrix G equivalent to g
% matrix version of impulse response
G = toeplitz([g', zeros(N-M,1)']', [g(1), zeros(1,M-1)] );
% predict true temperature from true heat prodcution
%qtrue2 = G*htrue;
% solve with damped least squares
GTG=G'*G;
GTGmax = max(max(abs(GTG)));
e2=1e-4*GTGmax;
hest2=(GTG+e2*eye(M,M))\(G'*qobs);
% solve for h given qobs and g using genralized least squares
% with a prior smoothness constraint using biconjugate gradient
e=10*max(abs(g)); % damping as a fraction of g
% Implement prior information of smoothness (2nd derivative = 0).
% Note the factor of 1/(Dt^2) in the derivative is omitted because
% the r.h.s. is zero.
% H is Toeplitz, so use spdiags() to create a sparse version of it.
% The [0,1,2] means the non-zero diagonals are the main diagonal
% (named 0) and twp immediateley to the right of it (the 1 and 2)
% The [ones(K,1), -2*ones(K,1), ones(K,1)] are the corresponding
% values of the diagonals. The sparse matrix is K by M.
K=M-2;
L=N+K;
H = spdiags([ones(M,1), -2*ones(M,1), ones(M,1)],[0,1,2],K,M);
h=zeros(K,1);
% set up RHS and solve
FTf = filterfunRHS(g,qobs,e,H,h,M);
mest=bicg(@filterfun,FTf,1e-10,3*L);
hest3=mest;
% 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(qobs)] );
plot(tdays(1:N),qobs,'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),hest2,'k-','LineWidth',2);
xlabel('time, days');
ylabel('hest(t)');
% plot simulated observations
figure(2);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), 0, 1.1*max(qobs)] );
plot(tdays(1:N),qobs,'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)');
% eda07_04: filter estimation problem
% example of estimating a filter r that predicts river
% discharge d from precipitation data g; that is, g*r=d
clear all;
global e g H;
% in this example, the filter length is chosen to be
% M<<N, to embody the idea that the response of the river
% (meaning the length M of the filter r) is short
M=25; % length of filter r
% load merged discharge and precipitation data
% column 1: time in inter count of days starting 1/1/2002
% column 2: Hudson River discharge in cfs
% column 3: Albany NY precipitation in inches
D = load('../Data/precip_discharge_data.txt');
% break out into separate variables, metric conversion
% time in days
t = D(:,1);
% discharge in cubic meters per second
d = D(:,2)/35.3146;
% precipitation in millimeters per day
g = D(:,3)*25.4;
% in this example, we use only a 101 day portion of the
% dataset, mainly so that when we plot it, storm events
% are clearly visible
N=101;
istart = 900;
t = t(istart:istart+N-1);
d = d(istart:istart+N-1);
g = g(istart:istart+N-1);
tmin = min(t);
tmax = max(t);
% the discharge and precipitation data
figure(1);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',3);
axis( [tmin, tmax, -20, 60 ] );
plot( t, g, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('precipitation (mm/day)');
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1000, 3000 ] );
plot( t, d, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('discharge (m3/s)');
% Given the filter equation g*r=d, solve for r given d and g
% using genralized least squares
% Note that the convolution relationship g*r=d is equivalent
% to the matrix equation G r = d where G is the (g*) operator,
% Since we are using filterfun(), we do not need to construct
% G. We do, however, need to implement the prior information
% equation Hr=h.
% This code can toggle between the prior information of r
% being small, its first derivative being small and its second
% derivative being small, depending on whether DERIVATIVE is 0,
% 1 or 2. (However, the results turn out not to be especially
% sensitive to the choice).
DERIVATIVE=1;
if( DERIVATIVE==2 )
% The prior information equation H*r = h is implemented
% for the smoothness condition of the second derivative
% being close to zero
e=2*max(abs(g)); % the coefficient was manually tuned
% to produce a reasonably smooth filter
K=M-2;
L=N+K;
H = spdiags([ones(M,1), -2*ones(M,1), ones(M,1)],[0,1,2],K,M);
h=zeros(K,1);
elseif (DERIVATIVE==1)
% The prior information equation H*r = h is implemented
% for the smoothness condition of the first derivative
% being close to zero
e=4.0*max(abs(g)); % the coefficient was manually tuned
% to produce a reasonably smooth filter
K=M-1;
L=N+K;
H = spdiags( [ones(M,1), -ones(M,1)], [0,1], K, M );
h=zeros(K,1);
else % DERRIVATIVE == 0
% The prior information equation H*r = h is implemented
% for the smallness condition of the filter is close to zero
e=0.1*max(abs(g)); % the coefficient was manually tuned
% to produce a reasonably small filter
K=M;
L=N+K;
H = spdiags([ones(M,1)],[0],K,M);
h=zeros(K,1);
end
% set up RHS and solve
FTf = filterfunRHS(g,d,e,H,h,M);
r=bicg(@filterfun,FTf,1e-12,3*L);
% plot the filter r
figure(2);
clf;
set(gca,'LineWidth',3);
hold on;
axis( [0, M, -5, 10 ] );
plot( [1:M]', r, 'k-', 'LineWidth', 3 );
plot( [1,M]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('filter r');
% plot the predicted discharge
dpre = conv(g,r);
dpre = dpre(1:N);
figure(1);
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1000, 3000 ] );
plot( t, dpre, 'k-', 'LineWidth', 2 );
% edama_07_05: prediction error filter (pef)
% this examples uses the for Neuse Riverr Hydrograph
% this code uses generalized least squares to implement
% the prior information that pef(1)=(-1). Two solutions
% are calculated, one via standard matrix algebra,
% (FTF)\(FTf) and the other via the biconjugate gradient method
clear all;
global e g H;
% read in data
D = load('../Data/neuse.txt');
% standard set-up of t-axis
t=D(:,1);
d=D(:,2);
N = length(d);
Dt = t(2)-t(1);
tmin = 0.0;
t = tmin + Dt*[0:N-1]';
tmax = t(N);
% sizes of various arrays
M=100; % length of pef is user tunable
K=1;
L=N+K;
% for a prediction error filter, pef, the data is the filter, g
g = d;
% prior information that first element of the pef is (-1)
e=(1000)*max(abs(g)); % large weight
H=zeros(K,M);
H(1,1)=1.0;
h=zeros(1,1);
h(1)=(-1.0);
% G needed only for standard colution
G=zeros(N,M);
G=toeplitz(g, [g(1), zeros(1,M-1)] );
% F needed only for standard solutiom
F=zeros(L,M);
F(1:N,:)=G;
F(N+1:N+K,:)=e*H;
% f needed only for standard solution
f=zeros(L,1);
f(1:N)=zeros(N,1);
f(N+1:N+K)=e*h;
% pef by standard matrix algebra
pef1=(F'*F)\(F'*f);
% FT f needed only for biconjugate gradient method
% FTf = GT 0 + e HT e h = e^2 [1 0 0 ... 0]T [-1 0 0 ... 0] = -e^2
FTf = zeros(M,1);
FTf(1) = -(e^2);
% solve for predictione error filter using generalized least squares
% constraint that pef(1)=(-1) implemented with generalized least squares
pef2=bicg(@filterfun,FTf,1e-12,3*L);
% prediction error
perror = conv(pef2,g);
% plot pef
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [ 0, Dt*(M-1), min(pef1), max(pef1)] );
plot( t(1:M), pef1, 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'pef(t) (std)' );
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [ 0, Dt*(M-1), min(pef1), max(pef1)] );
plot( t(1:M), pef2, 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'pef(t)( bicg)' );
% plot hydrograph data and prediction error for first year
figure(2);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
Np = 365;
axis( [tmin, tmin+Dt*Np, -5000, 15000]' );
hold on;
plot( t(1:Np), d(1:Np), 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'd(t)' );
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmin+Dt*Np, -5000, 15000]' );
plot( t(1:Np), perror(1:Np), 'k-', 'LineWidth', 2 );
xlabel( 'time t, days' );
ylabel( 'e(t)' );
% edama_07_06: Script supporting CribSheet 7.2
% Least squares filter estimation of filter m
% obeying g*m=d
clear all;
global e g H;
% standard set-up of time axis, t
N=101;
Dt = 1;
t = Dt*[0:N-1]';
tmin = t(1);
tmax = t(end);
% make some synthetic data g
g = random('exp',1,N,1);
% true filter is decaying exponential
M=10;
c=5.0*Dt;
mtrue=0.3*exp(-t(1:M)/c);
% true data d = g (convolve) m, truncated to N
dtrue = conv(g,mtrue);
dtrue = dtrue(1:N);
% observed data is true data plus noise
sigmad = 0.1;
d = dtrue + random('Normal',0,sigmad,N,1);
% prior infoirmation of smoothness (2nd dev = 0)
K = M-2;
L = N+K;
e = 0.001; % damping = sigma_d/sigma_h
col = ones(M,1)/(Dt^2);
H = spdiags([col,-2*col,col],[-1,0,1],K,M);
h = zeros(M-2,1);
% build RHS and solve
FTf = filterfunRHS(g,d,e,H,h,M);
mest=bicg(@filterfun,FTf,1e-12,3*L);
% predicted data
dpre = conv(g,mest);
dpre = dpre(1:N);
e = d-dpre;
E=e'*e;
% plot g, d
figure(1);
clf;
subplot(2,1,1)
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -5, 5] );
plot( t, g, 'k-', 'LineWidth', 2 );
xlabel('t');
ylabel('g');
subplot(2,1,2)
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -5, 5] );
plot( t, d, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('d');
% plot predicted data
subplot(2,1,2)
set(gca,'LineWidth',2);
hold on;
plot( t, dpre, '-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );
% plot filter m
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [tmin, t(M), 0, 1] );
plot( t(1:M), mtrue, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('m');
% plot estimated filtes
figure(2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, t(M), 0.0, 0.5] );
plot( t(1:M), mest, '-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );
% edama_07_07: inverse filter via z-transforms
% use z-transforms to:
% check if a filter is minumum phase
% represent it as a cascade of length-2 filters
% compute the inverse filter
% define a filter
g = [6, 5.5, 5.3, 5, 4, 3, 2]';
N = length(g);
gN = g(N);
% find roots of g
r = roots(flipud(g));
% check size of roots
badroots=0;
for j = [1:N-1]
if( abs(r(j)) <= 1 )
badroots=badroots+1;
end
end
fprintf('%d roots not outside the unit circle\n', badroots);
% rebuild filter from its roots
gp = zeros(1,1);
gp(1)=gN;
for j = [1:N-1]
gp = conv(gp,[-r(j),1]');
end
% there should be no imiginary part
gp = real(gp);
% check that it worked
E1 = (g-gp)'*(g-gp);
fprintf('error of filter rebuilt from its roots: %.2f',E1);
% construct inverse filter, gi, of length Ni
Ni = 50;
gi = zeros(Ni,1);
gi(1)=1/gN;
% filter cascade, one filter per loop
for j = [1:N-1]
% construct inverse filter of a length-2 filter
tmp = zeros(Ni,1);
tmp(1) = 1;
for k = [2:Ni]
tmp(k)=tmp(k-1)/r(j);
end
tmp = -tmp/r(j);
gi = conv(gi,tmp);
gi=gi(1:Ni);
end
% delete imaginary part (which should be zero)
gi = real(gi);
% test quality of inverse filter
gig=conv(g,gi);
gig=gig(1:Ni);
d = zeros(Ni,1);
d(1)=1;
E2 = (gig-d)'*(gig-d);
fprintf('error of inverse filter: %.2f\n',E2);
% plot results
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, Ni, min(g), max(g)] );
plot( [1:Ni]', [g',zeros(Ni-N,1)']', 'k-', 'LineWidth', 2);
xlabel('element j');
ylabel('g(j)');
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, Ni, min(gi), max(gi)] );
plot( [1:Ni]', gi, 'k-', 'LineWidth', 2);
xlabel('element j');
ylabel('ginv(j)');
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, Ni, -1.2, 1.2] );
plot( [1:Ni]', gig, 'k-', 'LineWidth', 2);
xlabel('element j');
ylabel('[g*ginv](j)');
% edama_07_08: recursive smoothing filter, by for-loops
% example of recursive smoothing filter
% this version used for-loops, for demonstration purposes
% the method in edama_07_09 is prefered
% standard set-up of t-axis
N=101;
Dt=1;
t=Dt*[0:N-1]';
% create two synthetic datasets
% h1, a unit spike
% h2, random noise with zero mean and unit variance
h1=zeros(N,1);
h1(floor(N/2))=1;
h2=random('Normal',0,1,N,1);
% apply recursive filter to h1 to yield q1
q1=zeros(N,1);
q1(1)=0.5*h1(1);
for j=[2:N]
q1(j)=0.5*(h1(j)+q1(j-1));
end
% apply recursive filter to h2 to yield q2
q2=zeros(N,1);
q2(1)=0.5*h2(1);
for j=[2:N]
q2(j)=0.5*(h2(j)+q2(j-1));
end
% plot
figure(1);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',2);
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, q1, 'k-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );
plot( t, h1, 'k.', 'LineWidth', 2 );
xlabel('time, t');
ylabel('h1(t) and q1(t)');
subplot(2,1,2);
hold on;
set(gca,'LineWidth',2);
axis( [t(1), t(N), -5, 5] );
plot( t, q2, 'k-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );
plot( t, h2, 'k.', 'LineWidth', 2 );
xlabel('time, t');
ylabel('h2(t) and q2(t)');
fprintf('sum of elements of the unit spike: %.2f', sum(h1));
fprintf('sum of elements of the filtered spike: %.2f', sum(q1));
% edama_07_09: recursive smoothing filter, by filter()
% example of recursive smoothing filter
% this version uses the filter() function
% standard set-up of time axis, t
N=101;
Dt=1;
t=Dt*[0:N-1]';
% create two synthetic datasets
% h1, a unit spike
% h2, random noise with zero mean and unit variance
h1=zeros(N,1);
h1(floor(N/2))=1;
h2=random('Normal',0,1,N,1);
u=[0.5,0]'; % create component filters
v=[1.0,-0.5];
q1 = filter(u, v, h1); % filter h1 to q1
q1 = q1(1:N); % truncate to length N
q2 = filter(u, v, h2); % filter h2 to q2
q2 = q2(1:N); % truncate to length N
% plot
figure(1);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',2);
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, q1, 'k-', 'LineWidth', 2', 'Color', [0.8,0.8,0.8] );
plot( t, h1, 'k.', 'LineWidth', 2 );
xlabel('time, t');
ylabel('h1(t) and q1(t)');
subplot(2,1,2);
hold on;
set(gca,'LineWidth',2);
axis( [t(1), t(N), -5, 5] );
plot( t, q2, 'k-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );
plot( t, h2, 'k.', 'LineWidth', 2 );
xlabel('time, t');
ylabel('h2(t) and q2(t)');
disp(sprintf('sum of elements of the unit spike: %f', sum(h1)));
disp(sprintf('sum of elements of the filtered spike: %f', sum(q1)));