[Scipy-svn] r5449 - trunk/scipy/fftpack/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Jan 12 21:01:12 EST 2009
Author: cdavid
Date: 2009-01-12 20:01:03 -0600 (Mon, 12 Jan 2009)
New Revision: 5449
Modified:
trunk/scipy/fftpack/tests/test_real_transforms.py
Log:
Few more tests.
Modified: trunk/scipy/fftpack/tests/test_real_transforms.py
===================================================================
--- trunk/scipy/fftpack/tests/test_real_transforms.py 2009-01-13 02:00:41 UTC (rev 5448)
+++ trunk/scipy/fftpack/tests/test_real_transforms.py 2009-01-13 02:01:03 UTC (rev 5449)
@@ -8,19 +8,26 @@
from scipy.io import loadmat
TDATA = loadmat(join(dirname(__file__), 'test.mat'),
- squeeze_me=True, struct_as_record=True)
+ squeeze_me=True, struct_as_record=True, mat_dtype=True)
X = [TDATA['x%d' % i] for i in range(8)]
Y = [TDATA['y%d' % i] for i in range(8)]
-def direct_dct(x):
+def direct_fft_dct(x, matlab=False):
"""Compute a Discrete Cosine Transform, type II.
- The DCT type II is defined as:
+ The DCT type II is defined as (matlab=False):
\forall u \in 0...N-1,
- dct(u) = a(u) sum_{i=0}^{N-1}{f(i)cos((i + 0.5)\pi u}
+ dct(u) = 2 * sum_{i=0}^{N-1}{f(i)cos((i + 0.5)\pi u/N}
- Where a(0) = sqrt(1/(4N)), a(u) = sqrt(1/(2N)) for u > 0
+ Or (matlab=True)
+
+ \forall u \in 0...N-1,
+ dct(u) = a(u) sum_{i=0}^{N-1}{f(i)cos((i + 0.5)\pi u/N}
+
+ Where a(0) = sqrt(1/N), a(u) = sqrt(2/N) for u > 0
+
+ This is the exact same formula as Matlab.
"""
x = np.asarray(x)
if not np.isrealobj(x):
@@ -30,13 +37,74 @@
y[1:2*n:2] = x
y[2*n+1::2] = x[-1::-1]
y = np.real(numfft(y))[:n]
- y[0] *= np.sqrt(.25 / n)
- y[1:] *= np.sqrt(.5 / n)
+ if matlab:
+ y[0] *= np.sqrt(.25 / n)
+ y[1:] *= np.sqrt(.5 / n)
return y
-def test_ref():
+def direct_dct(x):
+ """Direct implementation (O(n^2)) of dct II.
+
+ dct(u) = 2 * sum_{i=0}^{N-1}{f(i)cos((i + 0.5)\pi u/N}
+
+ Note that it is not 'normalized'
+ """
+ n = x.size
+ a = np.empty((n, n), dtype = x.dtype)
+ for i in xrange(n):
+ for j in xrange(n):
+ a[i, j] = x[j] * np.cos(np.pi * (0.5 + j) * i / n)
+
+ return 2 * a.sum(axis = 1)
+
+def fdct(x):
+ """Compute a 'Fast' Discrete Cosine Transform, type II, using a N point fft
+ instead of a direct 4n point DFT
+
+ \forall u \in 0...N-1,
+ dct(u) = sum_{i=0}^{N-1}{f(i)cos((i + 0.5)\pi u/N}
+
+ See 'A Fast Cosine Transform in One and Two Dimensions', by J. Makhoul, in
+ IEEE Transactions on acoustics, speech and signal processing.
+
+ Note that it is not 'normalized'
+ """
+ x = np.asarray(x)
+ n = x.size
+ v = np.empty(x.size, x.dtype)
+ if (n/2) * 2 == n:
+ iseven = True
+ else:
+ iseven = False
+ cut = (n-1)/2 + 1
+ v[:cut] = x[::2]
+ if iseven:
+ v[cut:] = x[-1:0:-2]
+ else:
+ v[cut:] = x[-2::-2]
+ t = 2 * numfft(v) * np.exp(-1j * np.pi * 0.5 / n * np.linspace(0, n-1, n))
+ v[:n/2+1] = np.real(t)[:n/2+1]
+ if iseven:
+ v[n/2+1:] = -np.imag(t)[n/2-1:0:-1]
+ else:
+ v[n/2+1:] = -np.imag(t)[n/2:0:-1]
+ return v
+
+def test_refs():
for i in range(len(X)):
- assert_array_almost_equal(direct_dct(X[i]), Y[i])
+ assert_array_almost_equal(direct_fft_dct(X[i], matlab=True), Y[i])
+ assert_array_almost_equal(direct_fft_dct(X[i], matlab=False), direct_dct(X[i]))
+ for i in range(len(X)):
+ x = X[i]
+ y = direct_fft_dct(x, matlab=True)
+ y[0] *= np.sqrt(x.size*4)
+ y[1:] *= np.sqrt(x.size*2)
+ assert_array_almost_equal(y, direct_dct(x))
+
+def test_fdct():
+ for i in range(len(X)):
+ assert_array_almost_equal(direct_dct(X[i]), fdct(X[i]))
+
if __name__ == "__main__":
np.testing.run_module_suite()
More information about the Scipy-svn
mailing list