Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(667)

Unified Diff: scipy/signal/cwt.py

Issue 154092: CWT support for SciPy (update)
Patch Set: Created 14 years, 4 months ago
Use n/p to move between diff chunks; N/P to move between comments. Please Sign in to add in-line comments.
Jump to:
View side-by-side diff with in-line comments
Download patch
« no previous file with comments | « scipy/signal/__init__.py ('k') | scipy/version.py » ('j') | no next file with comments »
Expand Comments ('e') | Collapse Comments ('c') | Show Comments Hide Comments ('s')
Index: scipy/signal/cwt.py
diff --git a/scipy/signal/cwt.py b/scipy/signal/cwt.py
new file mode 100644
index 0000000000000000000000000000000000000000..9f954c8040194180834347ea5d79e205b532b58b
--- /dev/null
+++ b/scipy/signal/cwt.py
@@ -0,0 +1,741 @@
+__all__ = ['cwt', 'ccwt', 'icwt', 'SDG', 'Morlet']
+
+import numpy as np
+from scipy.fftpack import fft, ifft, fftshift
+
+class Wavelet(object):
+ '''
+ Class for Wavelet object
+
+ The Wavelet object holds the wavelet coefficients as well as information on
+ how they were obtained.
+ '''
+
+ def __init__(self,wt,wavelet,weighting_function,signal_dtype,deep_copy=True):
+ '''
+ Initialization of Wavelet object.
+
+ Parameters
+ ----------
+ wt : array
+ array of wavelet coefficients
+ wavelet : object
+ mother wavelet object used in the creation of `wt`
+ weighting_function : function
+ function used in the creation of `wt`
+ signal_dtype : dtype
+ dtype of signal used in the creation of `wt`
+ deep_copy : bool
+ if true (default), the mother wavelet object used in the creation of
+ the wavelet object will be fully copied and accessible through
+ wavelet.motherwavelet; if false, wavelet.motherwavelet will be a
+ reference to the motherwavelet object.
+
+ Returns
+ -------
+
+ Returns an instance of the Wavelet class.
+ '''
+
+ from copy import deepcopy
+ self.coefs = wt[:,0:wavelet.len_signal]
+
+ if wavelet.len_signal != wavelet.len_wavelet:
+ self._pad_coefs = wt[:,wavelet.len_signal:]
+ else:
+ self._pad_coefs = None
+ if deep_copy:
+ self.motherwavelet = deepcopy(wavelet)
+ else:
+ self.motherwavelet = wavelet
+
+ self.weighting_function = weighting_function
+ self._signal_dtype = signal_dtype
+
+ def get_gws(self):
+ """
+ Calculate Global Wavelet Spectrum as defined in Torrence and Compo (1998)
+ """
+
+ gws = self.get_wavelet_var()
+
+ return gws
+
+
+ def get_wes(self):
+ """
+ Calculate Wavelet Energy Spectrum
+ """
+
+ from scipy.integrate import trapz
+
+ coef = 1. / (self.motherwavelet.fc * self.motherwavelet.cg)
+
+ wes = coef * trapz(np.power(np.abs(self.coefs),2),axis=1);
+
+ return wes
+
+ def get_wps(self):
+ """
+ Calculate Wavelet Power Spectrum
+ """
+
+ wps = 1./ (self.motherwavelet.len_signal) * self.get_wes()
+
+ return wps
+
+ def get_wavelet_var(self):
+ """
+ Calculate Wavelet Variance (a.k.a. the Global Wavelet Spectrum of
+ Torrence and Compo (1998))
+
+ References
+
+ Torrence, C., and G. P. Compo, 1998: A Practical Guide to Wavelet
+ Analysis. Bulletin of the American Meteorological Society, 79, 1,
+ pp. 61-78.
+ """
+
+ coef = self.motherwavelet.cg * self.motherwavelet.fc
+
+ wvar = (coef / self.motherwavelet.len_signal) * self.get_wes()
+
+ return wvar
+
+ def scalogram(self,show_coi=False,show_wps=False,ts = None,time = None, use_period = True, ylog_base = None, xlog_base = None,origin='top',figname = None):
+ """
+ Creates a simple plot of scalogram, with optional wavelet power specturm and
+ time series of the transformed signal.
+
+ Parameters
+ ----------
+ show_coi : bool
+ set to true to see Cone of Influence
+ show_wps : bool
+ set to true to see the Wavelet Power Spectrum
+ ts : array
+ 1D array containing time series data used in wavelet transform. If set,
+ time series will be plotted.
+ time : array
+ 1D array containing time information
+ use_period : bool
+ set to true to see figures use period instead of scale
+ ylog_base : float
+ if a log scale is desired, set `ylog_base` as float. (for log 10, set
+ ylog_base = 10)
+ xlog_base : float
+ if a log scale is desired, set `xlog_base` as float. (for log 10, set
+ xlog_base = 10) *note that this option is only valid for the wavelet power
+ spectrum figure.
+ origin : 'top' or 'bottom'
+ set origin of scale axis to top or bottom of figure
+
+ Returns
+ -------
+
+ None
+
+ Examples
+ --------
+
+ Create instance of SDG mother wavelet, normalized, using 10 scales and the
+ center frequency of the Fourier transform as the characteristic frequency.
+ Then, perform the continuous wavelet transform and plot the scalogram.
+
+ x = numpy.arange(0,2*numpy.pi,numpy.pi/8.)
+ data = numpy.sin(x**2)
+ scales = numpy.arange(10)
+
+ mother_wavelet = SDG(len_signal = len(data), scales = np.arange(10), normalize = True, fc = 'center')
+ wavelet = cwt(data, mother_wavelet)
+ wave_coefs.scalogram(origin = 'bottom')
+ """
+
+ import matplotlib.pyplot as plt
+ import matplotlib.cm as cm
+ from pylab import poly_between
+
+ if ts is not None:
+ show_ts = True
+ else:
+ show_ts = False
+
+ if not show_wps and not show_ts:
+ #only show scalogram
+ figrow = 1
+ figcol = 1
+ elif show_wps and not show_ts:
+ #show scalogram and wps
+ figrow = 1
+ figcol = 4
+ elif not show_wps and show_ts:
+ #show scalogram and ts
+ figrow = 2
+ figcol = 1
+ else:
+ #show scalogram, wps, and ts
+ figrow = 2
+ figcol = 4
+
+ if time is None:
+ x = np.arange(self.motherwavelet.len_signal)
+ else:
+ x = time
+
+ if use_period:
+ y = self.motherwavelet.scales / self.motherwavelet.fc
+ else:
+ y = self.motherwavelet.scales
+
+ fig = plt.figure(figsize=(16, 12),dpi=160)
+ ax1 = fig.add_subplot(figrow,figcol,1)
+
+ # if show wps, give 3/4 space to scalogram, 1/4 to wps
+ if show_wps:
+ # create temp axis at 3 or 4 col of row 1
+ axt = fig.add_subplot(figrow,figcol,3)
+ # get location of axtmp and ax1
+ axt_pos = axt.get_position()
+ ax1_pos = ax1.get_position()
+ axt_points = axt_pos.get_points()
+ ax1_points = ax1_pos.get_points()
+ # set axt_pos left bound to that of ax1
+ axt_points[0][0] = ax1_points[0][0]
+ ax1.set_position(axt_pos)
+ fig.delaxes(axt)
+
+ if show_coi:
+ # coi_coef is defined using the assumption that you are using
+ # period, not scale, in plotting - this handles that behavior
+ if use_period:
+ coi = self.motherwavelet.get_coi() / self.motherwavelet.fc / self.motherwavelet.sampf
+ else:
+ coi = self.motherwavelet.get_coi()
+
+ coi[coi==0]=y.min() - 0.1*y.min()
+ xs,ys = poly_between(np.arange(0,len(coi)),np.max(y),coi)
+ ax1.fill(xs,ys,'k',alpha=0.4,zorder = 2)
+
+ contf=ax1.contourf(x,y,np.abs(self.coefs)**2)
+ fig.colorbar(contf, ax=ax1, orientation = 'vertical',format='%2.1f')
+
+ if ylog_base is not None:
+ ax1.axes.set_yscale('log', basey=ylog_base)
+
+ if origin is 'top':
+ ax1.set_ylim((y[-1],y[0]))
+ elif origin is 'bottom':
+ ax1.set_ylim((y[0],y[-1]))
+ else:
+ raise OriginError('`origin` must be set to "top" or "bottom"')
+
+ ax1.set_xlim((x[0],x[-1]))
+ ax1.set_title('scalogram')
+ ax1.set_ylabel('time')
+ if use_period:
+ ax1.set_ylabel('period')
+ ax1.set_xlabel('time')
+ else:
+ ax1.set_ylabel('scales')
+ if time is not None:
+ ax1.set_xlabel('time')
+ else:
+ ax1.set_xlabel('sample')
+
+ if show_wps:
+ ax2 = fig.add_subplot(figrow,figcol,4,sharey=ax1)
+ if use_period:
+ ax2.plot(self.get_wps(),y,'k')
+ else:
+ ax2.plot(self.motherwavelet.fc * self.get_wps(),y,'k')
+
+ if ylog_base is not None:
+ ax2.axes.set_yscale('log', basey=ylog_base)
+ if xlog_base is not None:
+ ax2.axes.set_xscale('log', basey=xlog_base)
+ if origin is 'top':
+ ax2.set_ylim((y[-1],y[0]))
+ else:
+ ax2.set_ylim((y[0],y[-1]))
+ if use_period:
+ ax2.set_ylabel('period')
+ else:
+ ax2.set_ylabel('scales')
+ ax2.grid()
+ ax2.set_title('wavelet power spectrum')
+
+ if show_ts:
+ ax3 = fig.add_subplot(figrow,2,3,sharex=ax1)
+ ax3.plot(x,ts)
+ ax3.set_xlim((x[0],x[-1]))
+ ax3.legend(['time series'])
+ ax3.grid()
+ # align time series fig with scalogram fig
+ t = ax3.get_position()
+ ax3pos=t.get_points()
+ ax3pos[1][0]=ax1.get_position().get_points()[1][0]
+ t.set_points(ax3pos)
+ ax3.set_position(t)
+ if (time is not None) or use_period:
+ ax3.set_xlabel('time')
+ else:
+ ax3.set_xlabel('sample')
+
+ if figname is None:
+ plt.show()
+ else:
+ plt.savefig(figname)
+ plt.close('all')
+
+def cwt(x,wavelet,weighting_function = lambda x: x**(-0.5), deep_copy = True):
+ """
+ Computes the continuous wavelet transform of x using the mother wavelet
+ `wavelet`.
+
+ This function computes the continuous wavelet transform of x using an
+ instance a mother wavelet object.
+
+ The cwt is defined as:
+
+ T(a,b) = w(a) integral(-inf,inf)(x(t) * psi*{(t-b)/a} dt
+
+ which is a convolution. In this algorithm, the convolution in the time
+ domain is implemented as a multiplication in the Fourier domain.
+
+ Parameters
+ ----------
+ x : 1D array
+ time series to be transformed by the cwt
+ wavelet : Instance of the MotherWavelet class
+ instance of the MotherWavelet class for a particular wavelet family
+ weighting_function: Function used to weight
+ Typically w(a) = a^(-0.5) is chosen as it ensures that the
+ wavelets at every scale have the same energy.
+ deep_copy : bool
+ if true (default), the mother wavelet object used in the creation of
+ the wavelet object will be fully copied and accessible through
+ wavelet.motherwavelet; if false, wavelet.motherwavelet will be a
+ reference to the motherwavelet object.
+
+ Returns
+ -------
+
+ Returns an instance of the Wavelet class. The coefficients of the transform
+ can be obtain by the coefs() method (i.e. wavelet.coefs() )
+
+ Examples
+ --------
+
+ Create instance of SDG mother wavelet, normalized, using 10 scales and the
+ center frequency of the Fourier transform as the characteristic frequency.
+ Then, perform the continuous wavelet transform and plot the scalogram.
+
+ x = numpy.arange(0,2*numpy.pi,numpy.pi/8.)
+ data = numpy.sin(x**2)
+ scales = numpy.arange(10)
+
+ mother_wavelet = SDG(len_signal = len(data), scales = np.arange(10), normalize = True, fc = 'center')
+ wavelet = cwt(data, mother_wavelet)
+ wave_coefs.scalogram()
+
+ References
+ ----------
+
+ Addison, P. S., 2002: The Illustrated Wavelet Transform Handbook. Taylor
+ and Francis Group, New York/London. 353 pp.
+ """
+
+ signal_dtype = x.dtype
+
+ if len(x) < wavelet.len_wavelet:
+ n = len(x)
+ x = np.resize(x, (wavelet.len_wavelet,))
+ x[n:] = 0
+
+ # Transform the signal and motherwavelet into the Fourier domain
+
+ xf=fft(x)
+ mwf=fft(wavelet.coefs.conj(),axis=1)
+
+ # Convolve (mult. in Fourier space)
+ wt_tmp=ifft(mwf*xf[np.newaxis,:],axis=1)
+
+ # shift output from ifft and multiply by weighting function
+ wt = fftshift(wt_tmp,axes=[1])*weighting_function(wavelet.scales[:,np.newaxis])
+
+ # if motherwavelet and signal are real, only keep real part of transform
+ wt=wt.astype(np.lib.common_type(wavelet.coefs,x))
+
+ return Wavelet(wt,wavelet,weighting_function,signal_dtype,deep_copy)
+
+def ccwt(x1,x2,wavelet):
+ '''
+ Compute the continuous cross-wavelet transform of 'x1' and 'x2' using the
+ mother wavelet 'wavelet', which is an instance of the MotherWavelet class.
+
+ Parameters
+ ----------
+
+ x1,x2 : 1D array
+ time series used to compute cross-wavelet transform
+ wavelet : Instance of the MotherWavelet class
+ instance of the MotherWavelet class for a particular wavelet family
+
+ Returns
+ -------
+
+ Returns an instance of the Wavelet class.
+ '''
+
+ xwt=cwt(x1,wavelet)*np.conjugate(cwt(x2,wavelet))
+
+ return xwt
+
+def icwt(wavelet):
+ """
+ Compute the inverse continuous wavelet transform.
+
+ Parameters
+ ----------
+ wavelet : Instance of the MotherWavelet class
+ instance of the MotherWavelet class for a particular wavelet family
+
+ Examples
+ --------
+
+ Use the Morlet mother wavelet to perform wavelet transform on 'data', then
+ use icwt to compute the inverse wavelet transform to come up with an esitmate
+ of data ('data2'). Note that data2 is not exactly equal data.
+
+ import matplotlib.pyplot as plt
+ from scipy.signal import SDG, Morlet, cwt, icwt, fft, ifft
+ import numpy as np
+
+ x = np.arange(0,2*np.pi,np.pi/64)
+ data = np.sin(8*x)
+ scales=np.arange(0.5,17)
+
+ mother_wavelet = Morlet(len_signal = len(data), scales = scales)
+ wave_coefs=cwt(data, mother_wavelet)
+ data2 = icwt(wave_coefs)
+
+ plt.plot(data)
+ plt.plot(data2)
+ plt.show()
+
+ References
+ ----------
+
+ Addison, P. S., 2002: The Illustrated Wavelet Transform Handbook. Taylor
+ and Francis Group, New York/London. 353 pp.
+ """
+ from scipy.integrate import trapz
+
+ # if origional wavelet was created using padding, make sure to include
+ # information that is missing after truncation (see self.coefs under __init__
+ # in class Wavelet.
+
+ if wavelet.motherwavelet.len_signal != wavelet.motherwavelet.len_wavelet:
+ full_wc = np.c_[wavelet.coefs,wavelet._pad_coefs]
+ else:
+ full_wc = wavelet.coefs
+ # get wavelet coefficents and take fft
+ wcf = fft(full_wc,axis=1)
+ # get motherwavelet coeffifientts and take fft
+ mwf = fft(wavelet.motherwavelet.coefs,axis=1)
+ # perform inverse continuous wavelet transform and make sure the result is the same type
+ # (real or complex) as the origional data used in the transform
+ x = (1. / wavelet.motherwavelet.cg) * trapz(fftshift(ifft(
+ wcf * mwf,axis=1),axes=[1])/(wavelet.motherwavelet.scales[:,np.newaxis]**2),
+ dx = 1 / wavelet.motherwavelet.sampf,axis=0)
+
+
+ return x[0:wavelet.motherwavelet.len_signal].astype(wavelet._signal_dtype)
+
+class MotherWavelet(object):
+ """
+ Class for MotherWavelets
+
+ Contains methods related to mother wavelets. Also used to ensure that new
+ mother wavelet objects contain the minimum requirements to be used in the
+ cwt related functions.
+ """
+
+ @staticmethod
+ def get_coefs(self):
+ """
+ raise error method for calculating mother wavelet coefficients is
+ missing! To follow the convention in the literature, please define your
+ COI coef as a function of period, not scale - this will ensure
+ compatibility with the scalogram method.
+ """
+
+ raise NotImplementedError('get_coefs needs to be implemented for the mother wavelet')
+
+ @staticmethod
+ def get_coi_coef(sampf):
+ """
+ raise error if Cone of Influence coefficient is not set in subclass wavelet
+ """
+
+ raise NotImplementedError('coi_coef needs to be implemented in subclass wavelet')
+
+ #add methods for computing cone of influence and mask
+ def get_coi(self):
+ """
+ Compute cone of influence
+ """
+
+ y1 = self.coi_coef * np.arange(0,self.len_signal/2)
+ y2 = -self.coi_coef * np.arange(0,self.len_signal/2)+y1[-1]
+ coi = np.r_[y1,y2]
+ self.coi = coi
+ return coi
+
+ def get_mask(self):
+ """
+ get mask for cone of influence.
+
+ Sets self.mask as an array of bools for use in np.ma.array('',mask=mask)
+ """
+
+ mask = np.ones(self.coefs.shape)
+ masks = self.coi_coef*self.scales
+ for s in range(0,len(self.scales)):
+ if (s != 0) and (int(np.ceil(masks[s])) < mask.shape[1]):
+ mask[s,np.ceil(int(masks[s])):-np.ceil(int(masks[s]))]=0
+ self.mask = mask.astype(bool)
+ return self.mask
+
+
+class SDG(MotherWavelet):
+ """
+ Class for the SDG MotherWavelet (a subclass of MotherWavelet).
+
+ SDG(self, len_signal = None, pad_to = None, scales = None, sampf = 1,
+ normalize = True, fc = 'bandpass')
+
+ Parameters
+ ----------
+
+ len_signal : int
+ length of time series to be decomposed
+ pad_to : int
+ pad time series to a total length `pad_to` using zero padding (note,
+ the signal will be zero padded automatically during continuous wavelet
+ transform if pad_to is set).
+ scales : array
+ array of scales used to initialize the mother wavelet
+ sampf : float
+ sample frequency of the time series to be decomposed
+ normalize : bool
+ If True, the normalized version of the mother wavelet will be used (i.e.
+ the mother wavelet will have unit energy)
+ fc : string
+ Characteristic frequency - use the 'bandpass' or 'center' frequency of
+ the Fourier spectrum of the mother wavelet to relate scale to period
+ (default is 'bandpass')
+
+ Returns
+ -------
+
+ Returns an instance of the MotherWavelet class which is used in the cwt and
+ icwt functions.
+
+ Examples
+ --------
+
+ Create instance of SDG mother wavelet, normalized, using 10 scales and the
+ center frequency of the Fourier transform as the characteristic frequency.
+ Then, perform the continuous wavelet transform and plot the scalogram.
+
+ x = numpy.arange(0,2*numpy.pi,numpy.pi/8.)
+ data = numpy.sin(x**2)
+ scales = numpy.arange(10)
+
+ mother_wavelet = SDG(len_signal = len(data), scales = np.arange(10),normalize = True, fc = 'center')
+ wavelet = cwt(data, mother_wavelet)
+ wave_coefs.scalogram()
+
+ Notes
+ -----
+
+ None
+
+ References
+ ----------
+ Addison, P. S., 2002: The Illustrated Wavelet Transform Handbook. Taylor
+ and Francis Group, New York/London. 353 pp.
+ """
+
+ def __init__(self,len_signal=None,pad_to=None,scales=None,sampf=1,normalize=True, fc = 'bandpass'):
+ self.name='second degree of a Gaussian (mexican hat)'
+ self.sampf = sampf
+ self.scales = scales
+ self.len_signal = len_signal
+ self.normalize = normalize
+
+ #set total length of wavelet to account for zero padding
+ if pad_to is None:
+ self.len_wavelet = len_signal
+ else:
+ self.len_wavelet = pad_to
+
+ #set admissibility constant
+ if normalize:
+ self.cg = 4 * np.sqrt(np.pi) / 3.
+ else:
+ self.cg = np.pi
+
+ #define characteristic frequency
+ if fc is 'bandpass':
+ self.fc = np.sqrt(5./2.)*self.sampf/(2*np.pi)
+ elif fc is 'center':
+ self.fc = np.sqrt(2.)*self.sampf/(2*np.pi)
+ else:
+ raise CharacteristicFrequencyError("fc = %s not defined"%(fc,))
+
+ # coi_coef defined under the assumption that period is used, not scale
+ self.coi_coef = 2*np.pi*np.sqrt(2./5.)*self.fc ;#Torrence and Compo 1998
+
+ #compute coefficients for the dilated mother wavelet
+
+ self.coefs = self.get_coefs()
+
+ def get_coefs(self):
+ """
+ Calculate the coefficients for the mother wavelet SDG
+ """
+
+ #Create array containing values used to evaluate the wavelet function
+ xi=np.arange(-self.len_wavelet/2.,self.len_wavelet/2.)
+
+ #find motherwavelet coefficients at each scale
+ xsd = -xi * xi / (self.scales[:,np.newaxis] * self.scales[:,np.newaxis])
+
+ if self.normalize is True:
+ c=2./(np.sqrt(3)*np.power(np.pi,0.25))
+ else:
+ c=1.
+
+ mw = c * (1. + xsd) * np.exp(xsd / 2.)
+
+ self.coefs = mw
+
+ return mw
+
+class Morlet(MotherWavelet):
+ """
+ Class for the SDG MotherWavelet (a subclass of MotherWavelet).
+
+ Morlet(self, len_signal = None, pad_to = None, scales = None,
+ sampf = 1, f0 = 0.849)
+
+ Parameters
+ ----------
+
+ len_signal : int
+ length of time series to be decomposed
+ pad_to : int
+ pad time series to a total length `pad_to` using zero padding (note,
+ the signal will be zero padded automatically during continuous wavelet
+ transform if pad_to is set).
+ scales : array
+ array of scales used to initialize the mother wavelet
+ sampf : float
+ sample frequency of the time series to be decomposed
+ f0 : float
+ central frequency of the Morlet mother wavelet. The Fourier spectrum of
+ the Morlet wavelet appears as a Gaussian centered on f0. f0 defaults
+ to a value of 0.849 (the angular frequency would be ~5.336)
+
+ Returns
+ -------
+ Returns an instance of the MotherWavelet class which is used in the cwt
+ and icwt functions.
+
+ Examples
+ --------
+
+ Create instance of Morlet mother wavelet using 10 scales, perform the
+ continuous wavelet transform, and plot the resulting scalogram.
+
+ x = numpy.arange(0,2*numpy.pi,numpy.pi/8.)
+ data = numpy.sin(x**2)
+ scales = numpy.arange(10)
+
+ mother_wavelet = Morlet(len_signal=len(data), scales = np.arange(10))
+ wavelet = cwt(data, mother_wavelet)
+ wave_coefs.scalogram()
+
+ Notes
+ -----
+
+ * Morlet wavelet is defined as having unit energy, so the `normalize` flag
+ will always be set to True
+
+ * The Morlet wavelet will always use f0 as it's characteristic frequency, so
+ fc is set as f0
+
+ References
+ ----------
+
+ Addison, P. S., 2002: The Illustrated Wavelet Transform Handbook. Taylor
+ and Francis Group, New York/London. 353 pp.
+
+ """
+
+ def __init__(self,len_signal=None,pad_to=None,scales=None,sampf=1,normalize=True, f0 = 0.849):
+
+ from scipy.integrate import trapz
+
+ self.sampf = sampf
+ self.scales = scales
+ self.len_signal = len_signal
+ self.normalize = True
+ self.name='Morlet'
+
+ #set total length of wavelet to account for zero padding
+ if pad_to is None:
+ self.len_wavelet = len_signal
+ else:
+ self.len_wavelet = pad_to
+
+ #define characteristic frequency
+ self.fc = f0
+
+ # Cone of influence coefficient
+ self.coi_coef = 2* self.sampf/(self.fc + np.sqrt(2. + self.fc**2) \
+ *np.sqrt(2)); #Torrence and Compo 1998 (in code)
+
+ #set admissibility constant
+ #based on the simpilified Morlet wavelet energy spectrum
+ #in Addison (2002), eqn (2.39) - sould be ok for f0 >0.84
+ f = np.arange(0.001,50,0.001)
+ y = 2.*np.sqrt(np.pi)*np.exp(-np.power((2.*np.pi*f-2.*np.pi*self.fc),2))
+ self.cg = trapz(y[1:]/f[1:])*(f[1]-f[0])
+
+ #compute coefficients for the dilated mother wavelet
+ self.coefs = self.get_coefs()
+
+ def get_coefs(self):
+ """
+ Calculate the coefficients for the mother wavelet SDG
+ """
+
+ #Create array containing values used to evaluate the wavelet function
+
+ xi=np.arange(-self.len_wavelet/2.,self.len_wavelet/2.)
+
+ #find motherwavelet coefficients at each scale
+
+ xsd = xi / (self.scales[:,np.newaxis])
+
+ mw = np.power(np.pi,-0.25) * \
+ (np.exp(np.complex(1j) * 2. * np.pi * self.fc * xsd) - \
+ np.exp(-np.power((2. * np.pi * self.fc),2) / 2.)) * \
+ np.exp(-np.power(xsd,2) / 2.)
+
+ self.coefs = mw
+
+ return mw
« no previous file with comments | « scipy/signal/__init__.py ('k') | scipy/version.py » ('j') | no next file with comments »

Powered by Google App Engine
RSS Feeds Recent Issues | This issue
This is Rietveld f62528b