[Scipy-svn] r6205 - in trunk/scipy/signal: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Jan 17 23:56:47 EST 2010
Author: josef
Date: 2010-01-17 22:56:46 -0600 (Sun, 17 Jan 2010)
New Revision: 6205
Modified:
trunk/scipy/signal/signaltools.py
trunk/scipy/signal/tests/test_signaltools.py
Log:
signal hilbert: correction to axis handling, new axis argument, with tests, ticket:1093
Modified: trunk/scipy/signal/signaltools.py
===================================================================
--- trunk/scipy/signal/signaltools.py 2010-01-16 10:29:43 UTC (rev 6204)
+++ trunk/scipy/signal/signaltools.py 2010-01-18 04:56:46 UTC (rev 6205)
@@ -1028,22 +1028,24 @@
return w
-def hilbert(x, N=None):
+def hilbert(x, N=None, axis=-1):
"""Compute the analytic signal.
- The transformation is done along the first axis.
+ The transformation is done along the last axis by default.
Parameters
----------
x : array-like
Signal data
N : int, optional
- Number of Fourier components. Default: ``x.shape[0]``
+ Number of Fourier components. Default: ``x.shape[axis]``
+ axis : int, optional
+
Returns
-------
- xa : ndarray, shape (N,) + x.shape[1:]
- Analytic signal of `x`
+ xa : ndarray
+ Analytic signal of `x`, of each 1d array along axis
Notes
-----
@@ -1054,6 +1056,8 @@
where ``F`` is the Fourier transform, ``U`` the unit step function,
and ``y`` the Hilbert transform of ``x``. [1]
+ changes in scipy 0.8.0: new axis argument, new default axis=-1
+
References
----------
.. [1] Wikipedia, "Analytic signal".
@@ -1062,13 +1066,13 @@
"""
x = asarray(x)
if N is None:
- N = len(x)
+ N = x.shape[axis]
if N <=0:
raise ValueError, "N must be positive."
if iscomplexobj(x):
print "Warning: imaginary part of x ignored."
x = real(x)
- Xf = fft(x,N,axis=0)
+ Xf = fft(x, N, axis=axis)
h = zeros(N)
if N % 2 == 0:
h[0] = h[N/2] = 1
@@ -1078,8 +1082,10 @@
h[1:(N+1)/2] = 2
if len(x.shape) > 1:
- h = h[:, newaxis]
- x = ifft(Xf*h)
+ ind = [newaxis]*x.ndim
+ ind[axis] = slice(None)
+ h = h[ind]
+ x = ifft(Xf*h, axis=axis)
return x
def hilbert2(x,N=None):
Modified: trunk/scipy/signal/tests/test_signaltools.py
===================================================================
--- trunk/scipy/signal/tests/test_signaltools.py 2010-01-16 10:29:43 UTC (rev 6204)
+++ trunk/scipy/signal/tests/test_signaltools.py 2010-01-18 04:56:46 UTC (rev 6205)
@@ -5,7 +5,7 @@
from numpy.testing import *
import scipy.signal as signal
-from scipy.signal import lfilter, correlate, convolve, convolve2d
+from scipy.signal import lfilter, correlate, convolve, convolve2d, hilbert
from numpy import array, arange
@@ -704,5 +704,85 @@
x = np.arange(6)
assert_array_equal(signal.decimate(x, 2, n=1).round(), x[::2])
+
+class TestHilbert:
+ def test_hilbert_theoretical(self):
+ #test cases by Ariel Rokem
+ decimal = 14
+
+ pi = np.pi
+ t = np.arange(0, 2*pi, pi/256)
+ a0 = np.sin(t)
+ a1 = np.cos(t)
+ a2 = np.sin(2*t)
+ a3 = np.cos(2*t)
+ a = np.vstack([a0,a1,a2,a3])
+
+ h = hilbert(a)
+ h_abs = np.abs(h)
+ h_angle = np.angle(h)
+ h_real = np.real(h)
+
+ #The real part should be equal to the original signals:
+ assert_almost_equal(h_real, a, decimal)
+ #The absolute value should be one everywhere, for this input:
+ assert_almost_equal(h_abs, np.ones(a.shape), decimal)
+ #For the 'slow' sine - the phase should go from -pi/2 to pi/2 in
+ #the first 256 bins:
+ assert_almost_equal(h_angle[0,:256], np.arange(-pi/2,pi/2,pi/256),
+ decimal)
+ #For the 'slow' cosine - the phase should go from 0 to pi in the
+ #same interval:
+ assert_almost_equal(h_angle[1,:256], np.arange(0,pi,pi/256), decimal)
+ #The 'fast' sine should make this phase transition in half the time:
+ assert_almost_equal(h_angle[2,:128], np.arange(-pi/2,pi/2,pi/128),
+ decimal)
+ #Ditto for the 'fast' cosine:
+ assert_almost_equal(h_angle[3,:128], np.arange(0,pi,pi/128), decimal)
+
+ #The imaginary part of hilbert(cos(t)) = sin(t) Wikipedia
+ assert_almost_equal(h[1].imag, a0, decimal)
+
+ def test_hilbert_axisN(self):
+ # tests for axis and N arguments
+ a = np.arange(18).reshape(3,6)
+ # test axis
+ aa = hilbert(a, axis=-1)
+ yield assert_equal, hilbert(a.T, axis=0), aa.T
+ # test 1d
+ yield assert_equal, hilbert(a[0]), aa[0]
+
+ # test N
+ aan = hilbert(a, N=20, axis=-1)
+ yield assert_equal, aan.shape, [3,20]
+ yield assert_equal, hilbert(a.T, N=20, axis=0).shape, [20,3]
+ #the next test is just a regression test,
+ #no idea whether numbers make sense
+ a0hilb = np.array(
+ [ 0.000000000000000e+00-1.72015830311905j ,
+ 1.000000000000000e+00-2.047794505137069j,
+ 1.999999999999999e+00-2.244055555687583j,
+ 3.000000000000000e+00-1.262750302935009j,
+ 4.000000000000000e+00-1.066489252384493j,
+ 5.000000000000000e+00+2.918022706971047j,
+ 8.881784197001253e-17+3.845658908989067j,
+ -9.444121133484362e-17+0.985044202202061j,
+ -1.776356839400251e-16+1.332257797702019j,
+ -3.996802888650564e-16+0.501905089898885j,
+ 1.332267629550188e-16+0.668696078880782j,
+ -1.192678053963799e-16+0.235487067862679j,
+ -1.776356839400251e-16+0.286439612812121j,
+ 3.108624468950438e-16+0.031676888064907j,
+ 1.332267629550188e-16-0.019275656884536j,
+ -2.360035624836702e-16-0.1652588660287j ,
+ 0.000000000000000e+00-0.332049855010597j,
+ 3.552713678800501e-16-0.403810179797771j,
+ 8.881784197001253e-17-0.751023775297729j,
+ 9.444121133484362e-17-0.79252210110103j ])
+ yield assert_almost_equal, aan[0], a0hilb, 14, 'N regression'
+
+
+
+
if __name__ == "__main__":
run_module_suite()
More information about the Scipy-svn
mailing list