% Laboratory in Oceanography - Data and Methods
% SMAST, UMass Dartmouth
%
% Examples from signal processing toolbox:
%   fft basics, Fourier spectra, variance preserving form, filtering 
%
% Written by Miles A. Sundermeyer, 4/2/09

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clc

printflag = 1;

keyboard
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do some simple ffts to see what it looks like in Matlab

% Consider the following idealized time series
x = 1 + cos(2*pi*[0:8]/9)

% The Fourier transform is computed in Matlab as follows
X = fft(x, length(x)); 			% forward fft
xnew = ifft(X);			% inverse fft

[x' fft(x)' xnew']

%ans =
%   2.0000    8.0000             	2.0000          
%   1.7071    4.0000 + 0.0000i   	1.7071          
%   1.0000    0.0000 - 0.0000i   	1.0000          
%   0.2929    0      - 0.0000i		0.2929          
%        0    0                		0          
%   0.2929    0      + 0.0000i		0.2929          
%   1.0000    0.0000 + 0.0000i   	1.0000          
%   1.7071    4.0000 - 0.0000i   	1.7071          

% Note:
% Imaginary parts are all zero - no sine componenet
% First fft value is freq (k-1) = 0, and cos(0) = 1 => fft gives 8 (num data points) * (mean(x))
% Second & seventh fft values have same real value, and represent cosine variability with 8 points,
%  i.e., freq of 2pi/8.  Compute 2*X2/N = amp of cosine variability in orig signal.
% Other terms are zero since zero energy at other freqs.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add a sine component and repeat
x = 1 + cos(2*pi*[0:7]/8) -2*sin(4*pi*[0:7]/8)

% The Fourier transform is computed in Matlab as follows
X = fft(x); 			% forward fft
xnew = ifft(X);			% inverse fft

[x' fft(x).' xnew']

% Note:
% Input signal changed, but output is similar
% X3 = 8i, X7 = -8i ... Xn and XN+2-n are complex conjugates
% Imag parts of X2 and X7 => sine w/ freq 2*2pi/N has amp 2*Y3/8 = 2.

% Also note: In general, frequencies represented by fft are:  2*pi(k-1)/N, k = 0:(N/2)
% 	zero freq (mean), 
%	2*pi*(1/N) (lowest) ... 2*pi*((N/2 - 1)/N) (highest = Nyquist freq)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear X		% need to re-use X for harmonic analysis below

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run through an example using stage data from Muddy Creek, Chatham, MA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set some tidal frequencies (roughly)
M2 = 0.140519e-3;		% (rad/s)
S2 = 0.145444e-3;
N2 = 0.137880e-3;
O1 = 0.675977e-4;
K1 = 0.729212e-4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thisfnm = 'MuddyCreekUpper15min';
disp([' Loading data from file: ',thisfnm]);
eval(['load ',thisfnm])

if(0)		% create synthetic time series in place of actual depth data
  N = length(depth);
  depth = mean(depth) + 0.1*rand(N,1) ...
  		+ 0.1*sum(sin(86400*jday*[M2 S2 N2 O1 K1])')';
end

jdaylims = [min(jday) max(jday)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start with some summary plots
figure(1)
clf

subplot(5,1,1)
plot(jday,depth,'.-')
%xlabel('Yearday');
ylabel('Depth (m)');
grid on
set(gca,'xlim',jdaylims);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(5,1,2)
plot(jday,temp,'.-')
%xlabel('Yearday');
ylabel(['T (',setstr(176),'C)']);
grid on
set(gca,'xlim',jdaylims);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(5,1,3)
plot(jday,salin,'.')
%xlabel('Yearday');
ylabel('S (ppt)');
grid on
set(gca,'xlim',jdaylims);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(5,1,4)
plot(jday,Chla,'.-')
%xlabel('Yearday');
ylabel('Chla (ug/L)');
grid on
set(gca,'xlim',jdaylims);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(5,1,5)
plot(jday,DO,'.-')
xlabel('Yearday');
ylabel('DO (mg/L)');
grid on
set(gca,'xlim',jdaylims);

ax = axis;
axkeep = ax;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if printflag
  print -djpeg signal_proc_1.jpg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now wish to understand decomposition into sines and cosines
% consider just depth(t)
figure(2)
clf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replot depth data
subplot(3,1,1)
plot(jday,depth,'.-')
hold on
xlabel('Yearday');
ylabel('Depth (m)');
title('Stage (m)')
grid on
set(gca,'xlim',jdaylims);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now do a regression for tides
tfreqHz = [M2 S2 N2 O1 K1]/(2*pi);	% (rad/s)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We will do detiding using the following model, y = Xb
%
% | u_1 |   | sin(wt_1) cos(wt_1) 1 |
% | u_2 |   | sin(wt_2) cos(wt_2) 1 |
% | u_3 |   | sin(wt_3) cos(wt_3) 1 |  | a |
% |  .  | = |     .         .     . |  | b |
% |  .  |   |     .         .     . |  | M |
% |  .  |   |     .         .     . |
% |     |   |                       |
%
%   y	  =		X		b

M2sin = sin(86400*M2*jday);
M2cos = cos(86400*M2*jday);
S2sin = sin(86400*S2*jday);
S2cos = cos(86400*S2*jday);
N2sin = sin(86400*N2*jday);
N2cos = cos(86400*N2*jday);
O1sin = sin(86400*O1*jday);
O1cos = cos(86400*O1*jday);
K1sin = sin(86400*K1*jday);
K1cos = cos(86400*K1*jday);
meanones = ones(size(jday));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create model array, X, for detiding
% M2 & mean only
X{1} = [M2sin M2cos meanones];
% M2, S2, mean
X{2} = [M2sin M2cos S2sin S2cos meanones];
% M2, S2, N2, mean
X{3} = [M2sin M2cos S2sin S2cos N2sin N2cos meanones];
% M2, S2, O1, mean
X{4} = [M2sin M2cos S2sin S2cos N2sin N2cos O1sin O1cos meanones];
% M2, S2, N2, O1, K1, mean
X{5} = [M2sin M2cos S2sin S2cos N2sin N2cos O1sin O1cos K1sin K1cos meanones];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the regression for successively more complex models
for n=1:5
  [B,BINT,R,RINT,STATS]=regress(depth,X{n});
  if exist('h')
    delete(h);
  end
  h = plot(jday,X{n}*B,'r-','linewidth',2)
  %disp(' Press any key to continue ...')
  %pause
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In last iteration of this, have M2, S2, N2, O1, K1 consituents in B
An = sqrt(sum([B(1:2:10)'; B(2:2:10)'].^2));
  
subplot(3,2,3)
loglog(86400*tfreqHz,An.^2,'r*');
axis([1e-2 1e2 1e-10 1e5]);
xlabel('frequency (cpd)')
ylabel('A_n^2 (m^2)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now suppose I create a model that has all possible frequencies between 
% Nyquist, 1/(2*deltat) and 1/range(jday)
npts = length(jday)/2;
% create array of freqs in (rad/s)
newfreqs = 2*pi / 86400 * linspace(1/range(jday), 1/(2*diff(jday(1:2))), npts);

% create new X array with all these
Xnew = [cos(86400*jday*newfreqs) sin(86400*jday*newfreqs) meanones];

if(0)
  [B,BINT,R,RINT,STATS]=regress(depth,Xnew);
  save regress_spectrum B BINT R RINT STATS
else
  load regress_spectrum
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% see what this looks like
subplot(3,1,1)
hnew = plot(jday,Xnew*B,'m-','linewidth',2)

% also plot as power spectrum
subplot(3,2,3)
hold on
Annew = sqrt(sum([B(1:length(newfreqs))'; B(length(newfreqs)+1:end-1)'].^2));
plot(86400/(2*pi)*newfreqs,Annew.^2,'m-');
grid on
title('Regress ''spectrum''')

% Note: the difference (vertical offset) between above 'An_new' spectrum 
% and Matlab's Fourier spectra is that in Matlab, the fft is multiplied by N.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the spectrum the "old fashioned" fft way
N = length(depth);
Xn = fft(depth,N);
Pxx = Xn.*conj(Xn)/N;
freqs = linspace(1/range(jday), 1/(2*diff(jday(1:2))), N/2);
% get rid of negative freq parts of Pxx
Pxx(N/2+1:end) = [];
% and double power of pos freq parts
Pxx(2:end) = 2*Pxx(2:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,2,4)
loglog(freqs,Pxx);
xlabel('frequency (cpd)')
ylabel('Pxx (m^2)')
title('''By Hand'' FFT spectrum')
axis([1e-2 1e2 1e-10 1e5]);
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% verify Parseval's Theorem while we're here
p = [sum(depth.^2) sum(Pxx)]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now use Matlab's 'spectrum'			% Note: this function is now obsolete in Matlab
Fs = 1/(86400*diff(jday(1:2)));			% sample frequency (1/s)
[P,F] = spectrum(depth,length(depth),0,[],Fs,0.95);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,2,5)
specplot(P,Fs*86400);
axis([1e-2 1e2 1e-10 1e5]);
grid on

hold on
%plot(F,F.*P(:,1),'b.')         % variance preserving form
plot(F*86400,P(:,1),'b')         % non-variance preserving form

set(gca,'xscale','log')
xlabel('frequency (cpd)')
ylabel('Pxx (m^2)')
title('Matlab''s Power spectrum')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now use Matlab's periodogram
WINDOW = ones(size(depth));
NFFT = length(depth);
[Pxx2,F] = periodogram(depth-mean(depth),WINDOW,NFFT,Fs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(3,2,6)
loglog(86400*F,Pxx2/N); 
axis([1e-2 1e2 1e-10 1e5]);
grid on

xlabel('frequency (cpd)')
ylabel('Pxx (m^2)')
title('Matlab ''Periodogram'' spectrum')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the freq's of the basic tidal constituents on the spectra
for nn=1:4 
  subplot(3,2,2+nn)
  hold on
  ax=axis;
  plot(M2*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(S2*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(N2*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(O1*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(K1*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
end 

if printflag
  orient tall
  print -djpeg signal_proc_2.jpg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now create a second plot of two of same spectra, but in variance preserving form
figure(3)
clf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot 'by hand' spectrum in same form as above
subplot(2,2,1)
loglog(freqs,Pxx);
xlabel('frequency (cpd)')
ylabel('Pxx (m^2)')
title('''By Hand'' FFT spectrum')
axis([1e-2 1e2 1e-10 1e5])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now same spectrum in variance preserving form
subplot(2,2,2)
semilogx(freqs,freqs'.*Pxx);
xlabel('frequency (cpd)')
ylabel('f * Pxx (m^2/day)')
title('''By Hand'' FFT spectrum - Variance Preserving')
axis([1e-2 1e2 0 1e-0])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the freq's of the basic tidal constituents on the spectra
for nn=1:2 
  subplot(2,2,nn)
  hold on
  ax=axis;
  plot(M2*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(S2*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(N2*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(O1*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
  plot(K1*24*3600/(2*pi)*[1 1],ax(3:4),'r:')
end 

if printflag
  print -djpeg signal_proc_3.jpg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstrate filtering
figure(4)
clf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
loglog(freqs,Pxx);
xlabel('frequency (cpd)')
ylabel('Pxx (m^2)')
title('''By Hand'' FFT spectrum')
axis([1e-2 1e2 1e-10 1e5]);
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define a cutoff in this frequency series
cutoffind = find(F <= 1/86400*0.5);
cutoffind = max(cutoffind)
order = 8;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create a Butterworth filter 'by hand'
G02 = 1;
G2 = G02./(1+(F/F(cutoffind)).^(2*order));

% plot this over the spectrum
hold on
loglog(86400*F,G2,'r');

if printflag
  print -djpeg signal_proc_4.jpg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use Matlab to make a low-pass Butterworth digital filter 
% for filtering out the tides
Wn = cutoffind/length(F);
[b,a] = butter(order,Wn);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% and freqz to look at the filter just created
[H,W] = freqz(b,a,length(F));
Wcpd = W/pi * (1/(2*diff(jday(1:2))));
H2 = H.*conj(H);

loglog(Wcpd,H2,'m.');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% can also use fvtool to look at the filter just created
Hd = dfilt.df1(b,a); 		% Discrete-time filter (DFILT) object
h2 = fvtool(Hd);  

% make the fvtool plot look like mine
set(gca,'xscale','log')

% plot the filter by hand over the same 'decibel' scale
hold on
loglog(W/pi,10*log10(H2),'m.');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implement the filter by hand using convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Had the following from above:
% N = length(depth);
% Xn = fft(depth,N);
% Pxx = Xn.*conj(Xn)/N;
% freqs = linspace(1/range(jday), 1/(2*diff(jday(1:2))), N/2);
% % get rid of negative freq parts of Pxx
% Pxx(N/2+1:end) = [];
% % and double power of pos freq parts
% Pxx(2:end) = 2*Pxx(2:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% expand the filter, H, out to its full symmetric form
H_full = abs([H; flipud(H(2:end-1))]);

% take product (in Fourier space) with Xn
Xn_filt = Xn.*H_full;

% take the inverse fft to get the filtered time domain signal
x_filt = ifft(Xn_filt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prep this for plotting in Fourier space
% get rid of negative freq parts of Xn_filt
Qxx = Xn_filt.*conj(Xn_filt)/N;
Qxx(N/2+1:end) = [];

% and double power of pos freq parts
Qxx(2:end) = 2*Qxx(2:end);

figure(4)
loglog(freqs,Qxx,'color',[0 0.5 0],'marker','.'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% implement the filter in time domain using Matlab's filter (or filtfilt)
x_filt_single = filter(b,a,depth);
x_filt_too = filtfilt(b,a,depth);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot this
figure(6)
clf
subplot(2,1,1)
plot(jday,depth)
hold on
plot(jday,x_filt,'g.','linewidth',2)
plot(jday,x_filt_single,'r-','linewidth',1)
plot(jday,x_filt_too,'m-','linewidth',2)
xlabel('Yearday');
ylabel('Depth (m)');
title('Muddy Creek - Raw and Low-Pass Filtered Stage');
grid on
set(gca,'xlim',jdaylims);

h = legend('raw data', 'convolution filter', '''filter.m''', '''filtfilt.m''','location','SouthEast')
set(h,'position',[0.6634 0.5017 0.2321 0.2002])

if printflag
  print -djpeg signal_proc_6.jpg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following is my best understanding of the latest functionality in Matlab
% for computing spectra - similar to obsolete spectrum.m function, using a 
% Hanning window as default

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the spectrum of 'depth' time series using Welch's method

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create a Welch spectral estimator
seglen = N/4;
overlap = 50;						% percent overlap for segments
h = spectrum.welch('Hann',seglen,overlap)		% create a Welch spectral estimator

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the spectrum
hpsd = psd(h,depth,'Fs',Fs);
Pw = hpsd.Data;
Fw = hpsd.Frequencies;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the original 'by hand' version, plus the windowed one
figure(7)
clf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
loglog(freqs,Pxx);
xlabel('frequency (cpd)')
ylabel('Pxx (m^2)')
axis([1e-2 1e2 1e-10 1e5]);
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on
loglog(86400*Fw,Pw/N,'m-','linewidth',3);

title('''By Hand'' FFT spectrum with Hanning windowed spectrum')

if printflag
  print -djpeg signal_proc_7.jpg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% look at this Hanning window in real space
hanning_wind = hanning(seglen);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(8)
clf

subplot(2,1,1)
plot(jday,depth)
hold on
xlabel('Yearday');
ylabel('Depth (m)');
grid on
set(gca,'xlim',jdaylims);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot Hanning window with above overlap
dt = diff(jday(1:2));
x = min(jday) + dt*[0:seglen-1]';		% first hanning window
counter=0;
while(1)
  counter = counter+1;
  if counter>1
    x = x + dt*(seglen*overlap/100);
  end
  if max(x)>max(jday)
    break
  end
  plot(x,hanning_wind,'m-','linewidth',3);
end

title('Muddy Creek Upper stage, with Hanning windows');

if printflag
  print -djpeg signal_proc_8.jpg
end