Changeset 1515
- Timestamp:
- 09/29/08 09:12:41 (2 months ago)
- Files:
Legend:
- Unmodified
- Added
- Removed
- Modified
- Copied
- Moved
trunk/talkbox/scikits/talkbox/transforms/dct.py
r1514 r1515 10 10 # - C implementation 11 11 import numpy as np 12 13 # Definition of DCT in 1d (II type) 14 # dct(u) = a(u) \sum_{i=0}^{N-1}{f(x)cos(\pi(x + 0.5)u} 15 16 def direct_dctii(x): 17 """Direct implementation (O(n^2)) of dct II. 18 19 Notes 20 ----- 21 22 Use it as a reference only, it is not suitable for any real computation.""" 23 n = x.size 24 a = np.empty((n, n), dtype = x.dtype) 25 for i in xrange(n): 26 for j in xrange(n): 27 a[i, j] = x[j] * np.cos(np.pi * (0.5 + j) * i / n) 28 29 a[0] *= np.sqrt(1. / n) 30 a[1:] *= np.sqrt(2. / n) 31 32 return a.sum(axis = 1) 33 34 def direct_dctii_2(x): 35 """Direct implementation (O(n^2)) of dct.""" 36 # We are a bit smarter here by computing the coefficient matrix directly, 37 # but still O(N^2) 38 n = x.size 39 40 a = np.cos(np.pi / n * np.linspace(0, n - 1, n)[:, None] 41 * np.linspace(0.5, 0.5 + n - 1, n)[None, :]) 42 a *= x 43 a[0] *= np.sqrt(1. / n) 44 a[1:] *= np.sqrt(2. / n) 45 46 return a.sum(axis = 1) 12 from scipy.fftpack import fft 47 13 48 14 def dctii(x): … … 74 40 y[1:2*n:2] = x 75 41 y[2*n+1::2] = x[-1::-1] 76 y = np.real( np.fft.fft(y))[:n]42 y = np.real(fft(y))[:n] 77 43 y[0] *= np.sqrt(.25 / n) 78 44 y[1:] *= np.sqrt(.5 / n) … … 80 46 81 47 if __name__ == "__main__": 48 from dct_ref import direct_dctii_2 82 49 a = np.linspace(0, 10, 11) 83 50 print direct_dctii_2(a) 84 a = np.linspace(0, 10, 11) 85 print direct_dctii_2(a) 86 b = dctii(a) 87 print b 51 print dctii(a) trunk/talkbox/scikits/talkbox/transforms/tests/test_dct.py
r1334 r1515 2 2 from numpy.testing import TestCase, assert_array_almost_equal 3 3 4 from scikits.talkbox.transforms.dct import direct_dctii, direct_dctii_2, dctii 4 from scikits.talkbox.transforms.dct_ref import direct_dctii, direct_dctii_2 5 from scikits.talkbox.transforms.dct import dctii 5 6 6 7 class TestDCTTypeII(TestCase): … … 27 28 28 29 def test_simple(self): 29 """Test 1d input, 2nd direct implementation."""30 """Test 1d input, fft implementation.""" 30 31 assert_array_almost_equal(self.Y0, dctii(self.X0)) 31 32 assert_array_almost_equal(self.Y1, dctii(self.X1))
