In [72]:
%pylab inline
import os
import sys
import scipy.linalg as splin
#--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 [73]:
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 [74]:
# [name, period, amplitude in m, phase in degrees]
harmonics = [['M2',12.4206012,0.298,256.0],['K1',23.93447213,0.031,187.9]]
nfreq = len(harmonics)
ndays = 2
dmin = 6
total_sec = ndays * 24 * 60 * 60
dsec = dmin * 60
nsamples = total_sec / dsec
t = linspace(0.,float(total_sec),nsamples) - 86400. / 4.
# [tidal efficiency, phase shift]
v_parms = array([[1.,0.],
[0.75,45.],
[0.25,180.]])
clrs = ['blue','green','red']
signals = []
for vp in v_parms:
r = random.uniform(low=-0.02, high=0.02, size=nsamples)
v = create_signal(harmonics,t,vp) #+ r
plot(t/86400., v, linewidth=1.0)
signals.append(v)
xlim(t[0]/86400.,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 [75]:
def harmonic_analysis(station,pharmonics,v):
npfreq = len(pharmonics)
npsamples = v.shape[0]
nf2 = 2 * npfreq
nf21, nf22 = nf2+1, nf2+2
aa = zeros((nf22, nf22), float)
ci = zeros(npfreq, float)
si = zeros(npfreq, float)
x = zeros(nf21, float)
trig = zeros(nf22, float)
amplitude = zeros(npfreq, float)
theta = zeros(nfreq, float)
trig[0] = 1.
for i in xrange(npsamples):
for j in xrange(1,nf2,2):
jh = (j-1) / 2
period = pharmonics[jh][1] * 60. * 60.
if period==0.0: print period
sigma = 2. * pi / period
arc = sigma * t[i]
j1 = j + 1
trig[j] = cos(arc)
trig[j1] = sin(arc)
trig[nf22-1] = v[i]
# fill matrix
for i1 in xrange(nf21):
for j1 in xrange(i1,nf22):
aa[i1,j1] = aa[i1,j1] + trig[i1] * trig[j1]
aa[j1,i1] = aa[i1,j1]
# scale diagonal and off diagonals
for i in xrange(nf21):
diagonal = aa[i,i]
aa[i,i] = 1.0
for j in xrange(nf22):
if (j != i):
aa[i,j] = aa[i,j] / diagonal
# solve equations
determinant = splin.det(aa[:nf21,:nf21])
if abs(determinant) < 1.e-20:
print 'singular or ill-conditioned matrix'
print 'determinant = (0}'.format(determinant)
ainv = splin.inv(aa[:nf21,:nf21])
# calculate x
for i in xrange(nf21):
x[i] = 0.
for j in xrange(nf21):
x[i] += ainv[i,j] * aa[j,nf22-1]
# calculate amplitude and phase
for i in xrange(nfreq):
i2 = (i * 2) + 1
ci[i] = x[i2]
si[i] = x[i2+1]
amplitude[i] = sqrt(ci[i] * ci[i] + si[i] * si[i])
theta[i] = 360. - arctan2(si[i], ci[i]) * 180. / pi #degrees
#--write the data
print station
for idx,[c,period,pa,ptheta] in enumerate(pharmonics):
print ' {0:5s} {1:15.7g} {2:15.7g} deg.'.format(c, amplitude[idx], theta[idx])
return amplitude, theta
In [76]:
nf2 = 2 * nfreq
nf21, nf22 = nf2+1, nf2+2
for iv,v in enumerate(signals):
station = 'station {0}'.format(iv)
amplitude, theta = harmonic_analysis(station,harmonics,v)
# calculate simulated value
harmonic_results = []
for idx,[c,period,a,tv] in enumerate(harmonics):
harmonic_results.append([c, period, amplitude[idx], theta[idx]])
vp = [1.,0.]
vr = create_signal(harmonic_results,t,vp)
residual = subtract(vr, v)
subplot(211)
plot(t/86400., vr, linewidth=1.0, color=clrs[iv])
plot(t[::5]/86400., v[::5], linewidth=0.0, marker='x', markersize=3)
subplot(212)
plot(v, residual, linewidth=0.0, marker='o', markersize=3)
subplot(211)
xlim(t[0]/86400.,t[-1]/86400.)
ylabel(r'GMT Z$_0$ (MLLW), meters')
xlabel('Time, days')
subplot(212)
#xlim(0,t[-1]/86400.)
ylabel('residual, meters')
xlabel(r'GMT Z$_0$ (MLLW), meters')
show()
In [76]: