diff --git a/okean/filters.py b/okean/filters.py index b8c837b..7a1e310 100644 --- a/okean/filters.py +++ b/okean/filters.py @@ -5,6 +5,7 @@ ''' import numpy as np +from scipy.signal import detrend as dtrend class plfilt(object): """docstring for plfilt""" @@ -39,19 +40,66 @@ def __init__(self, dt=1.0): self.pl33 = np.interp(self.filter_time, self._dt, self._pl33) self.pl33 /= self.pl33.sum() - def __call__(self, u, t=None, mode='valid'): - uf = np.convolve(u, self.pl33, mode=mode) + def __call__(self, u, t=None, mode='valid', detrend = False): + if detrend: + u_dt = dtrend(u, type = 'linear') + du = u - u_dt + u = u_dt + else: + du = 0 + uf = np.convolve(u, self.pl33, mode=mode)+du if t is None: return uf else: - tf = t[self.Nt-1:-self.Nt+1] - return uf, tf + tf = {'valid' : t[self.Nt-1:-self.Nt+1], 'same' : t} + return uf, tf[mode] -def pl33(u,dt=1,t=None): - a=plfilt(dt) - return a(u,t) +def pl33(u,dt=1,t=None, extend_mode = 'same', detrend = False): + """Applies a tide-killer type convolution filter. + + The pl33 filter is described on p. 21, Rosenfeld (1983), WHOI + Technical Report 85-35. Filter half amplitude period = 33 h., half + power period = 38 h. + + Parameters + ---------- + u : (M,) array_like + Time series signal to filter. + dt : float (optional) + Sampling rate in multiples of the unit used in `t`. Default + is 1. This usually is represented in multiples of hrs. If the + series is sampled every hour and is set to 1 h, the filter will apply + a 33 h half amplitude period (`T`). `dt` can be set different from the + actual sampling rate to change the cutoff period, e.g. if sampling + rate is equal to 1h and dt is set to 0.5, `T` will be equal to 66 h. + t : (M,) array_like (optional) + Time stamp series with the same length of `u`. + extend_mode : {`valid`, `same`}, (optional) + Keyword for the type of padding made by the convolution function. + For a window W = 1/dt*33, N = M-2*W. + `same`: + Mode `same` returns output of length max(M, N). Boundary effects + are still visible. + `valid`: + Mode `valid` returns output of length max(M, N) - min(M, N) + 1. + The convolution product is only given for points where the signals + overlap completely. Values outside the signal boundary have no + effect. + detrend : bool + Apply linear detrend prior to applying convolution. + + Returns + ------- + uf : numpy.ndarray + Filtered signal. + tf : numpy.ndarray + Time stamps cut to the same length as the signal. + """ + + a=plfilt(dt) + return a(u,t, mode = extend_mode, detrend = detrend) #if __name__ == '__main__': # u = np.zeros((1000,),'d')