Aluode/PerceptionLabPortable
0
1# Copyright (c) 2006-2012 Filip Wasilewski <http://en.ig.ma/>
2# Copyright (c) 2012-2016 The PyWavelets Developers
3# <https://github.com/PyWavelets/pywt>
4# See COPYING for license details.
5
6"""
7Other wavelet related functions.
8"""
9
10
11import warnings
12
13import numpy as np
14from numpy.fft import fft
15
16from ._extensions._pywt import ContinuousWavelet, DiscreteContinuousWavelet, Wavelet
17
18__all__ = ["integrate_wavelet", "central_frequency",
19 "scale2frequency", "frequency2scale", "qmf",
20 "orthogonal_filter_bank",
21 "intwave", "centrfrq", "scal2frq", "orthfilt"]
22
23
24_DEPRECATION_MSG = ("`{old}` has been renamed to `{new}` and will "
25 "be removed in a future version of pywt.")
26
27
28def _integrate(arr, step):
29 integral = np.cumsum(arr)
30 integral *= step
31 return integral
32
33
34def intwave(*args, **kwargs):
35 msg = _DEPRECATION_MSG.format(old='intwave', new='integrate_wavelet')
36 warnings.warn(msg, DeprecationWarning)
37 return integrate_wavelet(*args, **kwargs)
38
39
40def centrfrq(*args, **kwargs):
41 msg = _DEPRECATION_MSG.format(old='centrfrq', new='central_frequency')
42 warnings.warn(msg, DeprecationWarning)
43 return central_frequency(*args, **kwargs)
44
45
46def scal2frq(*args, **kwargs):
47 msg = _DEPRECATION_MSG.format(old='scal2frq', new='scale2frequency')
48 warnings.warn(msg, DeprecationWarning)
49 return scale2frequency(*args, **kwargs)
50
51
52def orthfilt(*args, **kwargs):
53 msg = _DEPRECATION_MSG.format(old='orthfilt', new='orthogonal_filter_bank')
54 warnings.warn(msg, DeprecationWarning)
55 return orthogonal_filter_bank(*args, **kwargs)
56
57
58def integrate_wavelet(wavelet, precision=8):
59 """
60 Integrate `psi` wavelet function from -Inf to x using the rectangle
61 integration method.
62
63 Parameters
64 ----------
65 wavelet : Wavelet instance or str
66 Wavelet to integrate. If a string, should be the name of a wavelet.
67 precision : int, optional
68 Precision that will be used for wavelet function
69 approximation computed with the wavefun(level=precision)
70 Wavelet's method (default: 8).
71
72 Returns
73 -------
74 [int_psi, x] :
75 for orthogonal wavelets
76 [int_psi_d, int_psi_r, x] :
77 for other wavelets
78
79
80 Examples
81 --------
82 >>> from pywt import Wavelet, integrate_wavelet
83 >>> wavelet1 = Wavelet('db2')
84 >>> [int_psi, x] = integrate_wavelet(wavelet1, precision=5)
85 >>> wavelet2 = Wavelet('bior1.3')
86 >>> [int_psi_d, int_psi_r, x] = integrate_wavelet(wavelet2, precision=5)
87
88 """
89 # FIXME: this function should really use scipy.integrate.quad
90
91 if type(wavelet) in (tuple, list):
92 msg = ("Integration of a general signal is deprecated "
93 "and will be removed in a future version of pywt.")
94 warnings.warn(msg, DeprecationWarning)
95 elif not isinstance(wavelet, (Wavelet, ContinuousWavelet)):
96 wavelet = DiscreteContinuousWavelet(wavelet)
97
98 if type(wavelet) in (tuple, list):
99 psi, x = np.asarray(wavelet[0]), np.asarray(wavelet[1])
100 step = x[1] - x[0]
101 return _integrate(psi, step), x
102
103 functions_approximations = wavelet.wavefun(precision)
104
105 if len(functions_approximations) == 2: # continuous wavelet
106 psi, x = functions_approximations
107 step = x[1] - x[0]
108 return _integrate(psi, step), x
109
110 elif len(functions_approximations) == 3: # orthogonal wavelet
111 phi, psi, x = functions_approximations
112 step = x[1] - x[0]
113 return _integrate(psi, step), x
114
115 else: # biorthogonal wavelet
116 phi_d, psi_d, phi_r, psi_r, x = functions_approximations
117 step = x[1] - x[0]
118 return _integrate(psi_d, step), _integrate(psi_r, step), x
119
120
121def central_frequency(wavelet, precision=8):
122 """
123 Computes the central frequency of the `psi` wavelet function.
124
125 Parameters
126 ----------
127 wavelet : Wavelet instance, str or tuple
128 Wavelet to integrate. If a string, should be the name of a wavelet.
129 precision : int, optional
130 Precision that will be used for wavelet function
131 approximation computed with the wavefun(level=precision)
132 Wavelet's method (default: 8).
133
134 Returns
135 -------
136 scalar
137
138 """
139
140 if not isinstance(wavelet, (Wavelet, ContinuousWavelet)):
141 wavelet = DiscreteContinuousWavelet(wavelet)
142
143 functions_approximations = wavelet.wavefun(precision)
144
145 if len(functions_approximations) == 2:
146 psi, x = functions_approximations
147 else:
148 # (psi, x) for (phi, psi, x)
149 # (psi_d, x) for (phi_d, psi_d, phi_r, psi_r, x)
150 psi, x = functions_approximations[1], functions_approximations[-1]
151
152 domain = float(x[-1] - x[0])
153 assert domain > 0
154
155 index = np.argmax(abs(fft(psi)[1:])) + 2
156 if index > len(psi) / 2:
157 index = len(psi) - index + 2
158
159 return 1.0 / (domain / (index - 1))
160
161
162def scale2frequency(wavelet, scale, precision=8):
163 """Convert from CWT "scale" to normalized frequency.
164
165 Parameters
166 ----------
167 wavelet : Wavelet instance or str
168 Wavelet to integrate. If a string, should be the name of a wavelet.
169 scale : scalar
170 The scale of the CWT.
171 precision : int, optional
172 Precision that will be used for wavelet function approximation computed
173 with ``wavelet.wavefun(level=precision)``. Default is 8.
174
175 Returns
176 -------
177 freq : scalar
178 Frequency normalized to the sampling frequency. In other words, for a
179 sampling interval of `dt` seconds, the normalized frequency of 1.0
180 corresponds to (`1/dt` Hz).
181
182 """
183 return central_frequency(wavelet, precision=precision) / scale
184
185def frequency2scale(wavelet, freq, precision=8):
186 """Convert from to normalized frequency to CWT "scale".
187
188 Parameters
189 ----------
190 wavelet : Wavelet instance or str
191 Wavelet to integrate. If a string, should be the name of a wavelet.
192 freq : scalar
193 Frequency, normalized so that the sampling frequency corresponds to a
194 value of 1.0.
195 precision : int, optional
196 Precision that will be used for wavelet function approximation computed
197 with ``wavelet.wavefun(level=precision)``. Default is 8.
198
199 Returns
200 -------
201 scale : scalar
202
203 """
204 return central_frequency(wavelet, precision=precision) / freq
205
206def qmf(filt):
207 """
208 Returns the Quadrature Mirror Filter(QMF).
209
210 The magnitude response of QMF is mirror image about `pi/2` of that of the
211 input filter.
212
213 Parameters
214 ----------
215 filt : array_like
216 Input filter for which QMF needs to be computed.
217
218 Returns
219 -------
220 qm_filter : ndarray
221 Quadrature mirror of the input filter.
222
223 """
224 qm_filter = np.array(filt)[::-1]
225 qm_filter[1::2] = -qm_filter[1::2]
226 return qm_filter
227
228
229def orthogonal_filter_bank(scaling_filter):
230 """
231 Returns the orthogonal filter bank.
232
233 The orthogonal filter bank consists of the HPFs and LPFs at
234 decomposition and reconstruction stage for the input scaling filter.
235
236 Parameters
237 ----------
238 scaling_filter : array_like
239 Input scaling filter (father wavelet).
240
241 Returns
242 -------
243 orth_filt_bank : tuple of 4 ndarrays
244 The orthogonal filter bank of the input scaling filter in the order :
245 1] Decomposition LPF
246 2] Decomposition HPF
247 3] Reconstruction LPF
248 4] Reconstruction HPF
249
250 """
251 if not (len(scaling_filter) % 2 == 0):
252 raise ValueError("`scaling_filter` length has to be even.")
253
254 scaling_filter = np.asarray(scaling_filter, dtype=np.float64)
255
256 rec_lo = np.sqrt(2) * scaling_filter / np.sum(scaling_filter)
257 dec_lo = rec_lo[::-1]
258
259 rec_hi = qmf(rec_lo)
260 dec_hi = rec_hi[::-1]
261
262 orth_filt_bank = (dec_lo, dec_hi, rec_lo, rec_hi)
263 return orth_filt_bank
264