Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/fftpack/pseudo_diffs.py : 15%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Differential and pseudo-differential operators.
3"""
4# Created by Pearu Peterson, September 2002
6__all__ = ['diff',
7 'tilbert','itilbert','hilbert','ihilbert',
8 'cs_diff','cc_diff','sc_diff','ss_diff',
9 'shift']
11from numpy import pi, asarray, sin, cos, sinh, cosh, tanh, iscomplexobj
12from . import convolve
14from scipy.fft._pocketfft.helper import _datacopied
17_cache = {}
20def diff(x,order=1,period=None, _cache=_cache):
21 """
22 Return kth derivative (or integral) of a periodic sequence x.
24 If x_j and y_j are Fourier coefficients of periodic functions x
25 and y, respectively, then::
27 y_j = pow(sqrt(-1)*j*2*pi/period, order) * x_j
28 y_0 = 0 if order is not 0.
30 Parameters
31 ----------
32 x : array_like
33 Input array.
34 order : int, optional
35 The order of differentiation. Default order is 1. If order is
36 negative, then integration is carried out under the assumption
37 that ``x_0 == 0``.
38 period : float, optional
39 The assumed period of the sequence. Default is ``2*pi``.
41 Notes
42 -----
43 If ``sum(x, axis=0) = 0`` then ``diff(diff(x, k), -k) == x`` (within
44 numerical accuracy).
46 For odd order and even ``len(x)``, the Nyquist mode is taken zero.
48 """
49 tmp = asarray(x)
50 if order == 0:
51 return tmp
52 if iscomplexobj(tmp):
53 return diff(tmp.real,order,period)+1j*diff(tmp.imag,order,period)
54 if period is not None:
55 c = 2*pi/period
56 else:
57 c = 1.0
58 n = len(x)
59 omega = _cache.get((n,order,c))
60 if omega is None:
61 if len(_cache) > 20:
62 while _cache:
63 _cache.popitem()
65 def kernel(k,order=order,c=c):
66 if k:
67 return pow(c*k,order)
68 return 0
69 omega = convolve.init_convolution_kernel(n,kernel,d=order,
70 zero_nyquist=1)
71 _cache[(n,order,c)] = omega
72 overwrite_x = _datacopied(tmp, x)
73 return convolve.convolve(tmp,omega,swap_real_imag=order % 2,
74 overwrite_x=overwrite_x)
77del _cache
80_cache = {}
83def tilbert(x, h, period=None, _cache=_cache):
84 """
85 Return h-Tilbert transform of a periodic sequence x.
87 If x_j and y_j are Fourier coefficients of periodic functions x
88 and y, respectively, then::
90 y_j = sqrt(-1)*coth(j*h*2*pi/period) * x_j
91 y_0 = 0
93 Parameters
94 ----------
95 x : array_like
96 The input array to transform.
97 h : float
98 Defines the parameter of the Tilbert transform.
99 period : float, optional
100 The assumed period of the sequence. Default period is ``2*pi``.
102 Returns
103 -------
104 tilbert : ndarray
105 The result of the transform.
107 Notes
108 -----
109 If ``sum(x, axis=0) == 0`` and ``n = len(x)`` is odd, then
110 ``tilbert(itilbert(x)) == x``.
112 If ``2 * pi * h / period`` is approximately 10 or larger, then
113 numerically ``tilbert == hilbert``
114 (theoretically oo-Tilbert == Hilbert).
116 For even ``len(x)``, the Nyquist mode of ``x`` is taken zero.
118 """
119 tmp = asarray(x)
120 if iscomplexobj(tmp):
121 return tilbert(tmp.real, h, period) + \
122 1j * tilbert(tmp.imag, h, period)
124 if period is not None:
125 h = h * 2 * pi / period
127 n = len(x)
128 omega = _cache.get((n, h))
129 if omega is None:
130 if len(_cache) > 20:
131 while _cache:
132 _cache.popitem()
134 def kernel(k, h=h):
135 if k:
136 return 1.0/tanh(h*k)
138 return 0
140 omega = convolve.init_convolution_kernel(n, kernel, d=1)
141 _cache[(n,h)] = omega
143 overwrite_x = _datacopied(tmp, x)
144 return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
147del _cache
150_cache = {}
153def itilbert(x,h,period=None, _cache=_cache):
154 """
155 Return inverse h-Tilbert transform of a periodic sequence x.
157 If ``x_j`` and ``y_j`` are Fourier coefficients of periodic functions x
158 and y, respectively, then::
160 y_j = -sqrt(-1)*tanh(j*h*2*pi/period) * x_j
161 y_0 = 0
163 For more details, see `tilbert`.
165 """
166 tmp = asarray(x)
167 if iscomplexobj(tmp):
168 return itilbert(tmp.real,h,period) + \
169 1j*itilbert(tmp.imag,h,period)
170 if period is not None:
171 h = h*2*pi/period
172 n = len(x)
173 omega = _cache.get((n,h))
174 if omega is None:
175 if len(_cache) > 20:
176 while _cache:
177 _cache.popitem()
179 def kernel(k,h=h):
180 if k:
181 return -tanh(h*k)
182 return 0
183 omega = convolve.init_convolution_kernel(n,kernel,d=1)
184 _cache[(n,h)] = omega
185 overwrite_x = _datacopied(tmp, x)
186 return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
189del _cache
192_cache = {}
195def hilbert(x, _cache=_cache):
196 """
197 Return Hilbert transform of a periodic sequence x.
199 If x_j and y_j are Fourier coefficients of periodic functions x
200 and y, respectively, then::
202 y_j = sqrt(-1)*sign(j) * x_j
203 y_0 = 0
205 Parameters
206 ----------
207 x : array_like
208 The input array, should be periodic.
209 _cache : dict, optional
210 Dictionary that contains the kernel used to do a convolution with.
212 Returns
213 -------
214 y : ndarray
215 The transformed input.
217 See Also
218 --------
219 scipy.signal.hilbert : Compute the analytic signal, using the Hilbert
220 transform.
222 Notes
223 -----
224 If ``sum(x, axis=0) == 0`` then ``hilbert(ihilbert(x)) == x``.
226 For even len(x), the Nyquist mode of x is taken zero.
228 The sign of the returned transform does not have a factor -1 that is more
229 often than not found in the definition of the Hilbert transform. Note also
230 that `scipy.signal.hilbert` does have an extra -1 factor compared to this
231 function.
233 """
234 tmp = asarray(x)
235 if iscomplexobj(tmp):
236 return hilbert(tmp.real)+1j*hilbert(tmp.imag)
237 n = len(x)
238 omega = _cache.get(n)
239 if omega is None:
240 if len(_cache) > 20:
241 while _cache:
242 _cache.popitem()
244 def kernel(k):
245 if k > 0:
246 return 1.0
247 elif k < 0:
248 return -1.0
249 return 0.0
250 omega = convolve.init_convolution_kernel(n,kernel,d=1)
251 _cache[n] = omega
252 overwrite_x = _datacopied(tmp, x)
253 return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
256del _cache
259def ihilbert(x):
260 """
261 Return inverse Hilbert transform of a periodic sequence x.
263 If ``x_j`` and ``y_j`` are Fourier coefficients of periodic functions x
264 and y, respectively, then::
266 y_j = -sqrt(-1)*sign(j) * x_j
267 y_0 = 0
269 """
270 return -hilbert(x)
273_cache = {}
276def cs_diff(x, a, b, period=None, _cache=_cache):
277 """
278 Return (a,b)-cosh/sinh pseudo-derivative of a periodic sequence.
280 If ``x_j`` and ``y_j`` are Fourier coefficients of periodic functions x
281 and y, respectively, then::
283 y_j = -sqrt(-1)*cosh(j*a*2*pi/period)/sinh(j*b*2*pi/period) * x_j
284 y_0 = 0
286 Parameters
287 ----------
288 x : array_like
289 The array to take the pseudo-derivative from.
290 a, b : float
291 Defines the parameters of the cosh/sinh pseudo-differential
292 operator.
293 period : float, optional
294 The period of the sequence. Default period is ``2*pi``.
296 Returns
297 -------
298 cs_diff : ndarray
299 Pseudo-derivative of periodic sequence `x`.
301 Notes
302 -----
303 For even len(`x`), the Nyquist mode of `x` is taken as zero.
305 """
306 tmp = asarray(x)
307 if iscomplexobj(tmp):
308 return cs_diff(tmp.real,a,b,period) + \
309 1j*cs_diff(tmp.imag,a,b,period)
310 if period is not None:
311 a = a*2*pi/period
312 b = b*2*pi/period
313 n = len(x)
314 omega = _cache.get((n,a,b))
315 if omega is None:
316 if len(_cache) > 20:
317 while _cache:
318 _cache.popitem()
320 def kernel(k,a=a,b=b):
321 if k:
322 return -cosh(a*k)/sinh(b*k)
323 return 0
324 omega = convolve.init_convolution_kernel(n,kernel,d=1)
325 _cache[(n,a,b)] = omega
326 overwrite_x = _datacopied(tmp, x)
327 return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
330del _cache
333_cache = {}
336def sc_diff(x, a, b, period=None, _cache=_cache):
337 """
338 Return (a,b)-sinh/cosh pseudo-derivative of a periodic sequence x.
340 If x_j and y_j are Fourier coefficients of periodic functions x
341 and y, respectively, then::
343 y_j = sqrt(-1)*sinh(j*a*2*pi/period)/cosh(j*b*2*pi/period) * x_j
344 y_0 = 0
346 Parameters
347 ----------
348 x : array_like
349 Input array.
350 a,b : float
351 Defines the parameters of the sinh/cosh pseudo-differential
352 operator.
353 period : float, optional
354 The period of the sequence x. Default is 2*pi.
356 Notes
357 -----
358 ``sc_diff(cs_diff(x,a,b),b,a) == x``
359 For even ``len(x)``, the Nyquist mode of x is taken as zero.
361 """
362 tmp = asarray(x)
363 if iscomplexobj(tmp):
364 return sc_diff(tmp.real,a,b,period) + \
365 1j*sc_diff(tmp.imag,a,b,period)
366 if period is not None:
367 a = a*2*pi/period
368 b = b*2*pi/period
369 n = len(x)
370 omega = _cache.get((n,a,b))
371 if omega is None:
372 if len(_cache) > 20:
373 while _cache:
374 _cache.popitem()
376 def kernel(k,a=a,b=b):
377 if k:
378 return sinh(a*k)/cosh(b*k)
379 return 0
380 omega = convolve.init_convolution_kernel(n,kernel,d=1)
381 _cache[(n,a,b)] = omega
382 overwrite_x = _datacopied(tmp, x)
383 return convolve.convolve(tmp,omega,swap_real_imag=1,overwrite_x=overwrite_x)
386del _cache
389_cache = {}
392def ss_diff(x, a, b, period=None, _cache=_cache):
393 """
394 Return (a,b)-sinh/sinh pseudo-derivative of a periodic sequence x.
396 If x_j and y_j are Fourier coefficients of periodic functions x
397 and y, respectively, then::
399 y_j = sinh(j*a*2*pi/period)/sinh(j*b*2*pi/period) * x_j
400 y_0 = a/b * x_0
402 Parameters
403 ----------
404 x : array_like
405 The array to take the pseudo-derivative from.
406 a,b
407 Defines the parameters of the sinh/sinh pseudo-differential
408 operator.
409 period : float, optional
410 The period of the sequence x. Default is ``2*pi``.
412 Notes
413 -----
414 ``ss_diff(ss_diff(x,a,b),b,a) == x``
416 """
417 tmp = asarray(x)
418 if iscomplexobj(tmp):
419 return ss_diff(tmp.real,a,b,period) + \
420 1j*ss_diff(tmp.imag,a,b,period)
421 if period is not None:
422 a = a*2*pi/period
423 b = b*2*pi/period
424 n = len(x)
425 omega = _cache.get((n,a,b))
426 if omega is None:
427 if len(_cache) > 20:
428 while _cache:
429 _cache.popitem()
431 def kernel(k,a=a,b=b):
432 if k:
433 return sinh(a*k)/sinh(b*k)
434 return float(a)/b
435 omega = convolve.init_convolution_kernel(n,kernel)
436 _cache[(n,a,b)] = omega
437 overwrite_x = _datacopied(tmp, x)
438 return convolve.convolve(tmp,omega,overwrite_x=overwrite_x)
441del _cache
444_cache = {}
447def cc_diff(x, a, b, period=None, _cache=_cache):
448 """
449 Return (a,b)-cosh/cosh pseudo-derivative of a periodic sequence.
451 If x_j and y_j are Fourier coefficients of periodic functions x
452 and y, respectively, then::
454 y_j = cosh(j*a*2*pi/period)/cosh(j*b*2*pi/period) * x_j
456 Parameters
457 ----------
458 x : array_like
459 The array to take the pseudo-derivative from.
460 a,b : float
461 Defines the parameters of the sinh/sinh pseudo-differential
462 operator.
463 period : float, optional
464 The period of the sequence x. Default is ``2*pi``.
466 Returns
467 -------
468 cc_diff : ndarray
469 Pseudo-derivative of periodic sequence `x`.
471 Notes
472 -----
473 ``cc_diff(cc_diff(x,a,b),b,a) == x``
475 """
476 tmp = asarray(x)
477 if iscomplexobj(tmp):
478 return cc_diff(tmp.real,a,b,period) + \
479 1j*cc_diff(tmp.imag,a,b,period)
480 if period is not None:
481 a = a*2*pi/period
482 b = b*2*pi/period
483 n = len(x)
484 omega = _cache.get((n,a,b))
485 if omega is None:
486 if len(_cache) > 20:
487 while _cache:
488 _cache.popitem()
490 def kernel(k,a=a,b=b):
491 return cosh(a*k)/cosh(b*k)
492 omega = convolve.init_convolution_kernel(n,kernel)
493 _cache[(n,a,b)] = omega
494 overwrite_x = _datacopied(tmp, x)
495 return convolve.convolve(tmp,omega,overwrite_x=overwrite_x)
498del _cache
501_cache = {}
504def shift(x, a, period=None, _cache=_cache):
505 """
506 Shift periodic sequence x by a: y(u) = x(u+a).
508 If x_j and y_j are Fourier coefficients of periodic functions x
509 and y, respectively, then::
511 y_j = exp(j*a*2*pi/period*sqrt(-1)) * x_f
513 Parameters
514 ----------
515 x : array_like
516 The array to take the pseudo-derivative from.
517 a : float
518 Defines the parameters of the sinh/sinh pseudo-differential
519 period : float, optional
520 The period of the sequences x and y. Default period is ``2*pi``.
521 """
522 tmp = asarray(x)
523 if iscomplexobj(tmp):
524 return shift(tmp.real,a,period)+1j*shift(tmp.imag,a,period)
525 if period is not None:
526 a = a*2*pi/period
527 n = len(x)
528 omega = _cache.get((n,a))
529 if omega is None:
530 if len(_cache) > 20:
531 while _cache:
532 _cache.popitem()
534 def kernel_real(k,a=a):
535 return cos(a*k)
537 def kernel_imag(k,a=a):
538 return sin(a*k)
539 omega_real = convolve.init_convolution_kernel(n,kernel_real,d=0,
540 zero_nyquist=0)
541 omega_imag = convolve.init_convolution_kernel(n,kernel_imag,d=1,
542 zero_nyquist=0)
543 _cache[(n,a)] = omega_real,omega_imag
544 else:
545 omega_real,omega_imag = omega
546 overwrite_x = _datacopied(tmp, x)
547 return convolve.convolve_z(tmp,omega_real,omega_imag,
548 overwrite_x=overwrite_x)
551del _cache