diff --git a/tests/test_wavelets.py b/tests/test_wavelets.py index b3abd6f..48ee968 100644 --- a/tests/test_wavelets.py +++ b/tests/test_wavelets.py @@ -46,7 +46,7 @@ def test_DOG(): npt.assert_almost_equal(wavelets.DOG(m=6).frequency(0), 0, 6) -test_data = np.loadtxt('tests/nino3data.asc', skiprows=3) +test_data = np.loadtxt('nino3data.asc', skiprows=3) nino_time = test_data[:, 0] nino_dt = np.diff(nino_time).mean() @@ -197,7 +197,7 @@ def test_reconstruction_freq(): def test_power_bias(): """See if the global wavelet spectrum is biased or not. - Wavelet transform a signal of 3 distinct fourier frequencies. + Wavelet transform a signal of 3 distinct Fourier frequencies. The power spectrum should contain peaks at the frequencies, all of which should be the same height. @@ -298,11 +298,11 @@ def test_plot_coi(): ax.yaxis.set_ticklabels(ticks.astype(str)) ax.set_ylim(64, 0.5) - # second y scale with equivalent fourier periods to scales + # second y scale with equivalent Fourier periods to scales # except with the ticks at the powers of 2 ax_fourier = ax.twinx() ax_fourier.set_yscale('log') - # match the fourier ticks to the scale ticks + # match the Fourier ticks to the scale ticks ax_fourier.set_yticks(ticks) ax_fourier.set_yticklabels(ticks.astype(str)) ax_fourier.set_ylabel('fourier period (years)') diff --git a/wavelets/transform.py b/wavelets/transform.py index dd8536a..96b5fa4 100644 --- a/wavelets/transform.py +++ b/wavelets/transform.py @@ -12,7 +12,7 @@ def cwt(data, wavelet=None, widths=None, dt=1, frequency=False, axis=-1): - """Continuous wavelet transform using the fourier transform + """Continuous wavelet transform using the Fourier transform convolution as used in Terrence and Compo. (as opposed to the direct convolution method used by @@ -38,7 +38,7 @@ def cwt(data, wavelet=None, widths=None, dt=1, frequency=False, axis=-1): The first parameter is time or frequency. The second is a width parameter, defining the size of the wavelet - (e.g. standard deviation of a gaussian). + (e.g. standard deviation of a Gaussian). The wavelet function, Y, should be such that Int[-inf][inf](|Y|^2) = 1 @@ -101,8 +101,8 @@ def cwt_time(data, wavelet, widths, dt, axis): norm = (dt / width) ** .5 wavelet_data = norm * wavelet(t, width) output[ind, :] = scipy.signal.fftconvolve(data, - wavelet_data[slices], - mode='same') + wavelet_data[slices], + mode='same') return output @@ -164,7 +164,7 @@ class WaveletTransform(object): A *complex* wavelet function will return information about both amplitude and phase and is better adapted for capturing - *osillatory behaviour*. + *oscillatory behaviour*. A *real* wavelet function returns only a single component and can be used to isolate *peaks or discontinuities*. @@ -175,7 +175,7 @@ class WaveletTransform(object): wavelet amplitude. The resolution of the wavelet function is determined by the - balance between the width in real and fourier space. + balance between the width in real and Fourier space. A narrow function in time will have good time resolution but poor frequency resolution and vice versa. @@ -197,10 +197,10 @@ class WaveletTransform(object): The peak wavelet response does not necessarily occur at 1 / s. If we wish to compare wavelet spectra at different scales with - each other and with fourier modes, we need a common set of + each other and with Fourier modes, we need a common set of units. - The equivalent fourier period is defined as where the wavelet + The equivalent Fourier period is defined as where the wavelet power spectrum reaches its maximum and can be found analytically. """ def __init__(self, data=None, time=None, dt=1, @@ -209,7 +209,7 @@ def __init__(self, data=None, time=None, dt=1, """Arguments: data - 1 dimensional input signal time - corresponding times for the input signal - not essential, but the coi will be calculated + not essential, but the COI will be calculated for time starting at zero. dt - sample spacing dj - scale resolution @@ -247,14 +247,14 @@ def __init__(self, data=None, time=None, dt=1, @property def fourier_period(self): - """Return a function that calculates the equivalent fourier + """Return a function that calculates the equivalent Fourier period as a function of scale. """ return getattr(self.wavelet, 'fourier_period') @property def fourier_periods(self): - """Return the equivalent fourier periods for the scales used.""" + """Return the equivalent Fourier periods for the scales used.""" return self.fourier_period(self.scales) @property @@ -270,7 +270,7 @@ def s0(self, value): def find_s0(self): """Find the smallest resolvable scale by finding where the - equivalent fourier period is equal to 2 * dt. For a Morlet + equivalent Fourier period is equal to 2 * dt. For a Morlet wavelet, this is roughly 1. """ dt = self.dt @@ -309,7 +309,7 @@ def compute_optimal_scales(self): choose s0 so that the equivalent Fourier period is 2 * dt. The choice of dj depends on the width in spectral space of - the wavelet function. For the morlet, dj=0.5 is the largest + the wavelet function. For the Morlet, dj=0.5 is the largest that still adequately samples scale. Smaller dj gives finer scale resolution. """ @@ -317,7 +317,7 @@ def compute_optimal_scales(self): # resolution dj = self.dj # smallest resolvable scale, chosen so that the equivalent - # fourier period is approximately 2dt + # Fourier period is approximately 2dt s0 = self.s0 # Largest scale @@ -329,7 +329,7 @@ def compute_optimal_scales(self): # TODO: use np.frompyfunc on this # TODO: can we just replace it with fftfreqs? def w_k(self, k=None): - """Angular frequency as a function of fourier index. + """Angular frequency as a function of Fourier index. If no k, returns an array of all the angular frequencies calculated using the length of the data. @@ -396,7 +396,7 @@ def reconstruction(self, scales=None): x_n = (dj * dt^(1/2)) / (C_d * Y_0(0)) \ * Sum_(j=0)^J { Re(W_n(s_j)) / s_j^(1/2) } - where the factor C_d comes from the recontruction of a delta + where the factor C_d comes from the reconstruction of a delta function from its wavelet transform using the wavelet function Y_0. This C_d is a constant for each wavelet function. @@ -481,7 +481,7 @@ def compute_Cdelta(self): """Compute the parameter C_delta (see self.C_d), used in reconstruction. See section 3.i of TC98. - FIXME: this doesn't work. TC98 gives 0.776 for the morlet + FIXME: this doesn't work. TC98 gives 0.776 for the Morlet wavelet with dj=0.125. """ dj = self.dj @@ -521,7 +521,7 @@ def wavelet_variance(self): Returns the variance of the input data. """ - # TODO: mask coi for calculation of wavelet_variance + # TODO: mask COI for calculation of wavelet_variance # is this possible? how does it change the factors? dj = self.dj dt = self.dt diff --git a/wavelets/wavelets.py b/wavelets/wavelets.py index bcc3348..87a410e 100644 --- a/wavelets/wavelets.py +++ b/wavelets/wavelets.py @@ -14,7 +14,7 @@ class Morlet(object): def __init__(self, w0=6): """w0 is the nondimensional frequency constant. If this is set too low then the wavelet does not sample very well: a - value over 5 should be ok, Terrence and Compo set it to 6. + value over 5 should be ok; Terrence and Compo set it to 6. """ self.w0 = w0 if w0 == 6: @@ -40,7 +40,8 @@ def time(self, t, s=1.0, complete=True): Returns ------- - complex: value of the morlet wavelet at the given time + out : complex + Value of the Morlet wavelet at the given time See Also -------- @@ -54,18 +55,18 @@ def time(self, t, s=1.0, complete=True): This commonly used wavelet is often referred to simply as the Morlet wavelet. Note that this simplified version can cause - admissibility problems at low values of w. + admissibility problems at low values of `w`. The complete version:: pi**-0.25 * (exp(1j*w*x) - exp(-0.5*(w**2))) * exp(-0.5*(x**2)) The complete version of the Morlet wavelet, with a correction - term to improve admissibility. For w greater than 5, the + term to improve admissibility. For `w` greater than 5, the correction term is negligible. Note that the energy of the return wavelet is not normalised - according to s. + according to `s`. The fundamental frequency of this wavelet in Hz is given by ``f = 2*s*w*r / M`` where r is the sampling rate. @@ -86,18 +87,29 @@ def time(self, t, s=1.0, complete=True): # Fourier wavelengths def fourier_period(self, s): - """Equivalent fourier period of morlet""" + """Equivalent Fourier period of Morlet""" return 4 * np.pi * s / (self.w0 + (2 + self.w0 ** 2) ** .5) # Frequency representation def frequency(self, w, s=1.0): - """Frequency representation of morlet. + """Frequency representation of Morlet. - s - scale - w - angular frequency + Parameters + ---------- + w : float + Angular frequency. If `s` is not specified, i.e. set to 1, + this can be used as the non-dimensional angular + frequency w * s. + s : float + Scaling factor. Default is 1. + + Returns + ------- + out : complex + Value of the Morlet wavelet at the given frequency """ x = w * s - # heaviside mock + # Heaviside mock Hw = np.array(w) Hw[w <= 0] = 0 Hw[w > 0] = 1 @@ -117,7 +129,7 @@ def coi(self, s): class Paul(object): def __init__(self, m=4): - """Initialise a Paul wavelet function of order m. + """Initialise a Paul wavelet function of order `m`. """ self.m = m @@ -131,14 +143,15 @@ def time(self, t, s=1.0): Parameters ---------- t : float - Time. If s is not specified, i.e. set to 1, this can be + Time. If `s` is not specified, i.e. set to 1, this can be used as the non-dimensional time t/s. s : float Scaling factor. Default is 1. Returns ------- - complex: value of the paul wavelet at the given time + out : complex + Value of the Paul wavelet at the given time The Paul wavelet is defined (in time) as:: @@ -159,7 +172,7 @@ def time(self, t, s=1.0): # Fourier wavelengths def fourier_period(self, s): - """Equivalent fourier period of Paul""" + """Equivalent Fourier period of Paul""" return 4 * np.pi * s / (2 * self.m + 1) # Frequency representation @@ -169,7 +182,7 @@ def frequency(self, w, s=1.0): Parameters ---------- w : float - Angular frequency. If s is not specified, i.e. set to 1, + Angular frequency. If `s` is not specified, i.e. set to 1, this can be used as the non-dimensional angular frequency w * s. s : float @@ -177,12 +190,13 @@ def frequency(self, w, s=1.0): Returns ------- - complex: value of the paul wavelet at the given time + out : complex + Value of the Paul wavelet at the given frequency """ m = self.m x = w * s - # heaviside mock + # Heaviside mock Hw = 0.5 * (np.sign(x) + 1) # prefactor @@ -208,7 +222,7 @@ def coi(self, s): class DOG(object): def __init__(self, m=2): - """Initialise a Derivative of Gaussian wavelet of order m.""" + """Initialise a Derivative of Gaussian wavelet of order `m`.""" if m == 2: # value of C_d from TC98 self.C_d = 3.541 @@ -223,7 +237,7 @@ def __call__(self, *args, **kwargs): def time(self, t, s=1.0): """ - Return a DOG wavelet, + Return a Derivative of Gaussian wavelet, When m = 2, this is also known as the "Mexican hat", "Marr" or "Ricker" wavelet. @@ -236,24 +250,24 @@ def time(self, t, s=1.0): and ``x = t / s``. Note that the energy of the return wavelet is not normalised - according to s. + according to `s`. Parameters ---------- t : float - Time. If s is not specified, this can be used as the + Time. If `s` is not specified, this can be used as the non-dimensional time t/s. s : scalar Width parameter of the wavelet. Returns ------- - float : value of the ricker wavelet at the given time - + out : float + Value of the DOG wavelet at the given time Notes ----- - The derivative of the gaussian has a polynomial representation: + The derivative of the Gaussian has a polynomial representation: from http://en.wikipedia.org/wiki/Gaussian_function: @@ -271,8 +285,8 @@ def time(self, t, s=1.0): x = t / s m = self.m - # compute the hermite polynomial (used to evaluate the - # derivative of a gaussian) + # compute the Hermite polynomial (used to evaluate the + # derivative of a Gaussian) He_n = scipy.special.hermitenorm(m) gamma = scipy.special.gamma @@ -282,16 +296,16 @@ def time(self, t, s=1.0): return const * function def fourier_period(self, s): - """Equivalent fourier period of derivative of gaussian""" + """Equivalent Fourier period of derivative of Gaussian""" return 2 * np.pi * s / (self.m + 0.5) ** .5 def frequency(self, w, s=1.0): - """Frequency representation of derivative of gaussian. + """Frequency representation of derivative of Gaussian. Parameters ---------- w : float - Angular frequency. If s is not specified, i.e. set to 1, + Angular frequency. If `s` is not specified, i.e. set to 1, this can be used as the non-dimensional angular frequency w * s. s : float @@ -299,8 +313,9 @@ def frequency(self, w, s=1.0): Returns ------- - complex: value of the derivative of gaussian wavelet at the - given time + out : complex + Value of the derivative of Gaussian wavelet at the + given time """ m = self.m x = s * w @@ -324,7 +339,7 @@ def coi(self, s): class Ricker(DOG): def __init__(self): """The Ricker, aka Marr / Mexican Hat, wavelet is a - derivative of gaussian order 2. + derivative of Gaussian order 2. """ DOG.__init__(self, m=2) # value of C_d from TC98