In [5]:
%pylab inline
import os
import sys
import math
import datetime
import scipy.fftpack as spfft
import pandas as pd
#--default matplotlib
rcParams['mathtext.default'] = 'regular'
rcParams['legend.fontsize'] = 8
rcParams['axes.labelsize'] = 8
rcParams['xtick.labelsize'] = 8
rcParams['ytick.labelsize'] = 8
rcParams['figure.subplot.wspace'] = 0.35
rcParams['figure.subplot.hspace'] = 0.35
Create synthetic signals
In [6]:
def create_signal(harmonics,t,vp):
nsamples = t.shape
v = zeros(nsamples, float)
for idx,[c,period,a,theta] in enumerate(harmonics):
period *= 60. * 60.
fr = 1. / period
omega = 2. * pi * fr
theta *= pi / 180.
te = vp[0]
phase = vp[1] * pi / 180.
v += te * a * cos(omega*t + theta - phase)
return v
In [7]:
# [name, period, amplitude in m, phase in degrees]
harmonics = [['M2',12.4206012,0.298,256.0],['K1',23.93447213,0.031,187.9]]
ndays = 14
dmin = 6
total_sec = ndays * 24 * 60 * 60
dsec = dmin * 60
nsamples = total_sec / dsec
t = linspace(0,total_sec,nsamples)
# [tidal efficiency, phase shift]
v_parms = array([[1.,0.],
[0.75,45.],
[0.25,180.]])
signals = []
for vp in v_parms:
v = create_signal(harmonics,t,vp)
plot(t/86400., v, linewidth=1.0)
signals.append(v)
xlim(0,t[-1]/86400.)
ylabel(r'GMT Z$_0$ (MLLW), meters')
xlabel('Time, days')
show()
Now lets see if we can recover the correct amplitude and phase of the original signal and the tidal efficiency and phase shift of the secondary signals
In [8]:
def find_nearest_idx(array,value):
idx = (abs(array-value)).argmin()
return idx
def calculate_fft(station,t,v,harmonics,base=None):
nsamples = t.size
scale = 1./float(nsamples/2)
t = linspace(0,total_sec,nsamples)
FFT = spfft.fft(v)
freqs = spfft.fftfreq(nsamples, t[1]-t[0])
#--calculate fft of base
if base is not None:
bf = spfft.fft(base)
#--
alpha = []
dtau = []
for idx,[c,period,a,theta] in enumerate(harmonics):
period *= 60. * 60.
fr = 1. / period
omega = 2. * pi * fr
idx_fr = find_nearest_idx(freqs,fr)
if base is None:
alpha.append(abs(FFT[idx_fr]*scale))
dtau.append(0.)
else:
alpha.append(abs(FFT[idx_fr])/abs(bf[idx_fr]))
dt = angle(bf[idx_fr]/FFT[idx_fr]) * 180. / pi #radians to degrees
dtau.append(dt)
#--write the data
print station
for idx,[c,period,a,theta] in enumerate(harmonics):
print ' {0:5s} {1:15.7g} {2:15.7g} deg.'.format(c, alpha[idx], dtau[idx])
#plot it
subplot(211)
plot(t/86400., v)
subplot(212)
#plot(freqs[0:nsamples/2],20*log10(abs(FFT[0:nsamples/2]*scale)),linewidth=0.5)
plot(freqs[0:nsamples/2],log10(abs(FFT[0:nsamples/2]*scale)),linewidth=0.5)
return alpha
tide = signals[0]
station_alpha = []
for idx,v in enumerate(signals):
if idx == 0:
station = 'tide'
base = None
else:
station = 'signal {0:02d}'.format(idx+1)
base = tide
alpha = calculate_fft(station,t,v,harmonics,base=base)
station_alpha.append(alpha)
#--reset x axis for subplot 1 to maximum extent
subplot(211)
ylabel('Time, days')
ylabel(r'GMT Z$_0$ (MLLW), meters')
xlim(0,t[-1]/86400.)
subplot(212)
for idx,[c,period,a,theta] in enumerate(harmonics):
period *= 60. * 60.
fr = 1. / period
vlines(fr,-5,0,linewidth=1,color='0.5')
text(fr, 0.1, c,ha='center',va='bottom',size=10)
ylim((-5,0))
ylabel(r'log$_{10}$(amplitude), meters')
xlim((0,0.0001))
xlabel('frequency, 1/sec')
show()
In [8]: