
%% AA *** continuous vs. discrete ***
close all
figure(1)
Fs = 60
dt = 1/Fs;
N = 1000;
ts = 1:N
ts = dt.*ts;
T = ts(N);
fl = 1/T;
fun2 = sin(2*pi*15*fl.*ts);
plot(ts,fun2,'k-','LineWidth',2)
grid on
hold on
axis([0 3 -1.5 1.5])
ylabel('amplitude (arb. units)')
title([num2str(15*fl),' Hz sine wave'])
xlabel('time (s)-->')
ylabel('(arbitrary units)')
R = input('press enter to continue:');
plot(ts,fun2,'r*','MarkerSize',12)
title(' is really a DISCRETE function')
R = input('press enter to continue')
cla
plot(ts,fun2,'r*','MarkerSize',12)
axis([0.8 1.5 -0.9 1])

%% *** BB this is foldover, sample rate is 60 (Hz), two
%     sin waves, one with f=5 and one with f=60+5 look
%     identical when plotted at 60 Hz (60 samples per second). ***
close all
clear i
figure(1)
ts = 0:1/60:1;
func1 = 2.*sin(2*pi*5.*ts);
func2 = 2.*sin(2*pi*65.*ts);
plot(ts,real(func1),'b:')
hold on
plot(ts,real(func1),'b*')
plot(ts,real(func2),'rs')
axis([0 1 -3 3])
grid on
xlabel('time (s)')
ylabel('amplitude')
title (' 5 Hz (blue) and 65 Hz (red) sampled at 60 Hz')


%close all
%% CC *** another look at foldover. ***
figure(2)
ts = 0:1/60:(1024-1)/60;
func1 = sin(2*pi*5.*ts);
func2 = sin(2*pi*65.*ts);
xfc1 = fft(func1,512);
xfc2 = fft(func2,512);
xmfc1 = abs(xfc1(1:256).^2);
xmfc2 = abs(xfc2(1:256).^2);
fn = 60*(1:512/2)/512;
fn2 = [0 fn(1:(512/2-1))];
subplot(2,1,1)
bar(fn2,xmfc1);
%axis([0 0.125 0 200000])
ylabel('magnitude of fft coefficient 5 Hz')
subplot(2,1,2)
bar(fn2,xmfc2);
%axis([0 0.125 0 200000])
ylabel('magnitude of fft coefficient 65 Hz')
xlabel('frequency (Hz)')

%% DD) *** an example to show orthogonality ***
close all
figure(1)
Fs = 60
dt = 1/Fs;
N = 1000;
ts = 1:N
ts = dt.*ts;
T = ts(N);
fl = 1/T;
fun1 = sin(2*pi*10*fl.*ts);
fun2 = sin(2*pi*15*fl.*ts);
subplot(2,1,1)
plot(ts,fun1,'b*','MarkerSize',2)
grid on
hold on
plot(ts,fun2,'rd','MarkerSize',2)
axis([0 20 -2 2])
ylabel('amplitude (arb. units)')
title([num2str(10*fl),' (10\nu_s) and ',num2str(15*fl),' (15\nu_s) Hz sine waves'])

cross = fun1.*fun2;
sum(cross)
subplot(2,1,2)
plot(ts,cross,'m*')
axis([0 20 -2 2])
grid on
ylabel('one sinusoid x other (arb. units)')
xlabel(' time (s)')
title(' An example of ORTHOGONALITY, point by point multiplication of red and blue, above')
text(4,1.5,[' the SUM of this, over an integer number of waveforms, is : ',num2str(sum(cross),'%0.3g')]);
% num2str(sum(cross),'1.3g')
mesg = [' the SUM of fun1 * fun1 is: ',num2str(sum(fun1.*fun1))]
mesg = [' the SUM of fun2 * fun2 is: ',num2str(sum(fun2.*fun2))]
%% EE) *** an example to show non-orthogonality ***
figure(2)
Fs = 60
dt = 1/Fs;
N = 1000;
ts = 1:N
ts = dt.*ts;
T = ts(N);
fl = 1/T;
fun1 = sin(2*pi*10.3*fl.*ts);
fun2 = sin(2*pi*15*fl.*ts);
subplot(2,1,1)
plot(ts,fun1,'LineWidth',2)
grid on
hold on
plot(ts,fun2,'r:','LineWidth',2)
axis([0 20 -2 2])
ylabel('amplitude (arb. units)')
title([num2str(10.3*fl),' (10.3\nu_s) and ',num2str(15*fl),' (15\nu_s) Hz sine waves'])

cross = fun1.*fun2;
sum(cross)
subplot(2,1,2)
plot(ts,cross,'m*')
axis([0 20 -2 2])
grid on
ylabel('one sinusoid multiplied by the other (arb. units)')
xlabel(' time (s)')
title(' An example of NON-ORTHOGONALITY')
text(4,1.5,[' the SUM of this, over an integer # of waveforms, is: ',num2str(sum(cross),'%0.3g')]);
% num2str(sum(cross),'0.3g')
text(4,-1.4,[' this is why we do FFTs ONLY for integer multiples of some '])
text(4,-1.75,[' fundamental frequency (to ensure ORTHOGONALITY)']);

mesg = [' the SUM of fun1 * fun1 is: ',num2str(sum(fun1.*fun1))]
mesg = [' the SUM of fun2 * fun2 is: ',num2str(sum(fun2.*fun2))]

%% FF) *** FFT example ***
Fs = 1000;                    % Sampling frequency
dt = 1/Fs;                     % Sample time
N = 1000;                     % Length of signal
t = (1:N).*dt;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
xx = sin(2*pi*120*t);
x = 1 + 0.7*sin(2*pi*50*t+pi/2) + 1.0*sin(2*pi*120*t); 
y = x + 2*randn(size(t));     % Sinusoids plus noise
% *** what is scale of noise?  Note that randn produces numbers from -1:1,
% and that the amplitude of the combined signal is about 1.7.  S/N ratio is
% then 1.7/2.  That's pretty noisy, as you can see. ***
figure(3)
subplot(2,1,1)
plot(Fs*t(1:50),xx(1:50),'r.')
hold on
plot(Fs*t(1:50),x(1:50),'b*')
plot(Fs*t(1:50),x(1:50),'b:')
title('x(t) = 120 Hz signal (red) and sum of 50 & 120 Hz sines (x,blue)')
grid on
subplot(2,1,2)
plot(Fs*t(1:50),y(1:50),'b*')
hold on
plot(Fs*t(1:50),y(1:50),'b:')
title('y(t) = Signal Corrupted with Zero-Mean Random Noise (y)')
xlabel('time (milliseconds)')
grid on


figure(4)
NFFT = 2^nextpow2(N); % Next power of 2 from length of y
%NFFT = 1024
% *** BEWARE differences between NFFT and L ***
X = fft(x,NFFT)/N;
Y = fft(y,NFFT)/N;
f = Fs/2*linspace(0,1,NFFT/2);

% Plot single-sided amplitude spectrum.
subplot(2,1,1)
bar(f,2*abs(X(1:NFFT/2))) 
title('Single-Sided Amplitude Spectrum of x(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
grid on
text(200,0.8,'note that, even with perfect sinusoidal input, signal is ')
text(200,0.7,'smeared into bars beside main FC bar')
subplot(2,1,2)
bar(f,2*abs(Y(1:NFFT/2))) 
title('Single-Sided Amplitude Spectrum of y(t) (not bad, eh?)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
grid on
text(200,0.8,'noise is obvious, note also some mis-estimation of FCs for ')
text(200,0.6,'two input sinusoids (look closely at 50 Hz signal) ')

%% GG) *** now try 'reconstructing' some tides from FFTs ***
close all
load 390_ffts
figure(5)
plot(ghr./(24*365),Sttide)
xlabel('year from 1998')
ylabel('meters')
title( 'Seattle NOAA hourly tide gauge data-- 1998')
axis([0 1 -1 4])
grid on
hold on
ns = 1;
ne = 1024;
nfft = 1024;

plot(ghr(ns:ne)./(24*365),Sttide(ns:ne),'r:')

% *** first, look at the first 1024 data points ***
X = fft(Sttide(ns:ne),nfft);
xmag = abs(X(1:nfft/2).^2);
figure(6)
subplot(3,1,1)
plot((ns:ne),Sttide(ns:ne));
axis([0 1200 -1 5])
ylabel('wl (m)')
xlabel(' hours')
title(' original data, first 1024 points')
grid on

fn = (1:nfft/2)/nfft;
fn2 = [0 fn(1:(nfft/2-1))];
pn2 = 1./fn2;
pn2(1) = 1500;
figure(7)
bar(pn2,xmag/nfft);
%axis([0 0.125 0 200000])
axis([0 60 0 200])
ylabel('magnitude of fft coefficient (cm)')
xlabel('period (hours/cycle) (DC omitted)')
hold on

%% HH) *** first reconstruct first 1024 points ***
% figure(1)
% *** find FC coeff magnitudes > 2000)
ipt = find(xmag>1000);
isz = size(ipt,1)
iarray = zeros(1,nfft);
t = ns:ne;
for k=1:isz
%  k = 2;
  xp = t.*-2i*pi.*(ipt(k)-1)./nfft;
  denom = nfft/2;
  if (k==1), denom = nfft;, end
  hk = X(ipt(k)).*exp(xp)./denom;
%  figure(2)
%  bar(fn2(ipt(k)),xmag(ipt(k)),colormap(cool))
  iarray = iarray + hk;
  figure(6)
  subplot(3,1,2)
  plot((nfft:-1:1),real(iarray));
  axis([0 1200 -1 5])
  ylabel('wl (m)')
  title(' synthesized data')
  grid on

  subplot(3,1,3)
  plot((1:nfft),real(hk));
  axis([0 1200 -1 1])
  ylabel('wl (m)')
  title(' next cosine wave added')
  grid on
  
  min(abs(hk))
  max(abs(hk))
  playerobj = audioplayer(8.*[hk,hk,hk,hk,hk,hk],2048);
  play(playerobj)
   %sound(8.*abs([hk,hk,hk,hk,hk,hk]),4096)

%  hold on
  r=input('hit CR for next?')
  end

%% I) *** now somewhere deeper in the data, for window of 4096 points
ns = 2048;
ne = 2048+4096-1;
nfft = 4096;
figure(5)
plot(ghr(ns:ne)./(24*365),Sttide(ns:ne),'g:')

% *** first, look at the first 1024 data points ***
X = fft(Sttide(ns:ne),nfft);
xmag = abs(X(1:nfft/2).^2);
figure(8)
subplot(3,1,1)
plot((ns:ne),Sttide(ns:ne));
axis([ns-20 ne+20 -1 5])
ylabel('wl (m)')
xlabel(' hours')
title(' original data')
grid on

fn = (1:nfft/2)/nfft;
fn2 = [0 fn(1:(nfft/2-1))];
pn2 = 1./fn2;
pn2(1) = 1500;
figure(9)
bar(pn2,xmag/nfft);
%axis([0 0.125 0 200000])
axis([0 60 0 200])
ylabel('magnitude of fft coefficient')
xlabel('period (hours/cycle) (DC omitted)')
hold on

% figure(1)
% *** find FC coeff magnitudes > 2000)
ipt = find(xmag>1000);
isz = size(ipt,1)
iarray = zeros(1,nfft);
t = ns:ne;
for k=1:isz
%  k = 2;
  xp = t.*-2i*pi.*(ipt(k)-1)./nfft;
  denom = nfft/2;
  if (k==1), denom = nfft;, end
  hk = X(ipt(k)).*exp(xp)./denom;
%  figure(2)
%  bar(fn2(ipt(k)),xmag(ipt(k)),colormap(cool))
  iarray = iarray + hk;
  figure(8)
  subplot(3,1,2)
  plot((nfft:-1:1),real(iarray));
  axis([ns-20 ne+20 -1 5])
  ylabel('wl (m)')
  title(' synthesized data')
  grid on

  subplot(3,1,3)
  plot((1:nfft),real(hk));
  axis([ns-20 ne+20 -1 1])
  ylabel('wl (m)')
  title(' next cosine wave added')
  grid on

  % ***This bit will play sped-up cosine so students can see and hear coeffs
  %   as they are added to synthesized data. ***
  min(abs(hk))
  max(abs(hk))
  playerobj = audioplayer(8.*[hk,hk,hk,hk,hk,hk],2048);
  play(playerobj)
   %sound(8.*abs([hk,hk,hk,hk,hk,hk]),4096)

%  hold on
  r=input('hit CR for next?')
  end

%%  J) *** an example about window lengths with Seattle hourly tide data ***
% *** shows how increasing the NFFT increases the frequency resolution ***
% (as in delta_f = fs/2*NFFT)
close all
load 390_ffts
figure(1)
plot(ghr./(24*365),Sttide)
xlabel('year from 1998')
ylabel('meters')
title( 'Seattle')
axis([0 1 -1 4])
grid on

figure(2)
Hs = spectrum.welch;
set(Hs,'WindowName','Blackman');
set(Hs,'SegmentLength',256);
hpsd = psd(Hs,Sttide,'Fs',1);
F = hpsd.Frequencies;
P = hpsd.Data;
% [P,F] = spectrum.welch(Sttide,256,0,blackman(256),1);
bar(F,P)
xlabel('frequency (cycles/hour)')
ylabel('cabs of FFT coeff')
title(' spectrum for N=256')
axis([0 0.125 0 500])
grid on

figure(3)
Hs = spectrum.welch;
set(Hs,'WindowName','Blackman');
set(Hs,'SegmentLength',512);
hpsd = psd(Hs,Sttide,'Fs',1);
F = hpsd.Frequencies;
P = hpsd.Data;
%[P,F] = spectrum.welch(Sttide,512,0,blackman(512),1);
bar(F,P)
xlabel('frequeny (cycles/hour)')
ylabel('cabs of FFT coeff')
title(' spectrum for N=512')
axis([0 0.125 0 500])
grid on

figure(4)
Hs = spectrum.welch;
set(Hs,'WindowName','Blackman');
set(Hs,'SegmentLength',1024);
hpsd = psd(Hs,Sttide,'Fs',1);
F = hpsd.Frequencies;
P = hpsd.Data;
%[P,F] = spectrum.welch(Sttide,1024,0,blackman(1024),1);
bar(F,P)
xlabel('frequeny (cycles/hour)')
ylabel('cabs of FFT coeff')
title(' spectrum for N=1024')
axis([0 0.125 0 500])
grid on







