[Scipy-svn] r5483 - in trunk/scipy/fftpack: . src tests

scipy-svn@scip... scipy-svn@scip...
Sun Jan 18 05:24:48 CST 2009


Author: cdavid
Date: 2009-01-18 05:24:37 -0600 (Sun, 18 Jan 2009)
New Revision: 5483

Modified:
   trunk/scipy/fftpack/fftpack.pyf
   trunk/scipy/fftpack/realtransforms.py
   trunk/scipy/fftpack/src/dct.c.src
   trunk/scipy/fftpack/tests/test_real_transforms.py
Log:
Add dct III.

Modified: trunk/scipy/fftpack/fftpack.pyf
===================================================================
--- trunk/scipy/fftpack/fftpack.pyf	2009-01-18 11:24:18 UTC (rev 5482)
+++ trunk/scipy/fftpack/fftpack.pyf	2009-01-18 11:24:37 UTC (rev 5483)
@@ -185,6 +185,17 @@
          integer optional,intent(c,in) :: normalize = 0
        end subroutine ddct2
 
+       subroutine ddct3(x,n,howmany,normalize)
+         ! y = ddct3(x[,n,normalize,overwrite_x])
+         intent(c) ddct3
+         real*8 intent(c,in,out,copy,out=y) :: x(*)
+         integer optional,depend(x),intent(c,in) :: n=size(x)
+         check(n>0&&n<=size(x)) n
+         integer depend(x,n),intent(c,hide) :: howmany = size(x)/n
+         check(n*howmany==size(x)) howmany
+         integer optional,intent(c,in) :: normalize = 0
+       end subroutine ddct3
+
        subroutine destroy_ddct2_cache()
          intent(c) destroy_ddct2_cache
        end subroutine destroy_ddct2_cache

Modified: trunk/scipy/fftpack/realtransforms.py
===================================================================
--- trunk/scipy/fftpack/realtransforms.py	2009-01-18 11:24:18 UTC (rev 5482)
+++ trunk/scipy/fftpack/realtransforms.py	2009-01-18 11:24:37 UTC (rev 5483)
@@ -64,6 +64,43 @@
     """
     return _dct(x, 2, n, axis, normalize=norm)
 
+def dct3(x, n=None, axis=-1, norm=None):
+    """
+    Return Discrete Cosine Transform (type III) of arbitrary type sequence x.
+
+    There are several definitions, we use the following:
+
+                          N-1
+        y[k] = x[0] + 2 * sum x[n]*cos(pi*(k+0.5)*n/N), 0 <= k < N.
+                          n=0
+
+    The DCT-III is the inverse of DCT-II up to a scaling factor.
+
+    Parameters
+    ----------
+    x : array-like
+        input array.
+    n : int, optional
+        Length of the transform.
+    axis : int, optional
+        axis over which to compute the transform.
+
+    Returns
+    -------
+    y : real ndarray
+
+    Notes
+    -----
+    The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up
+    to a factor 2*N.
+    
+    Examples
+    --------
+    >>> x = np.linspace(0, 9, 10)
+    >>> np.testing.assert_array_almost_equal(dct3(dct2(x)) / 20, x)
+    """
+    return _dct(x, 3, n, axis, normalize=norm)
+
 def _dct(x, type, n=None, axis=-1, overwrite_x=0, normalize=None):
     """
     Return Discrete Cosine Transform of arbitrary type sequence x.
@@ -97,6 +134,8 @@
         f = _fftpack.ddct1
     elif type == 2:
         f = _fftpack.ddct2
+    elif type == 3:
+        f = _fftpack.ddct3
     else:
         raise ValueError("Type %d not understood" % type)
 

Modified: trunk/scipy/fftpack/src/dct.c.src
===================================================================
--- trunk/scipy/fftpack/src/dct.c.src	2009-01-18 11:24:18 UTC (rev 5482)
+++ trunk/scipy/fftpack/src/dct.c.src	2009-01-18 11:24:37 UTC (rev 5483)
@@ -111,4 +111,26 @@
     }
 }
 
+void @pref@dct3(@type@ * inout, int n, int howmany, int normalize)
+{
+    int i;
+    @type@ *ptr = inout;
+    @type@ *wsave = NULL;
+
+    wsave = caches_@pref@dct2[get_cache_id_@pref@dct2(n)].wsave;
+
+    for (i = 0; i < howmany; ++i, ptr += n) {
+        @pref@cosqf_(&n, ptr, wsave);
+
+    }
+
+    switch (normalize) {
+        case DCT_NORMALIZE_NO:
+            break;
+        default:
+            fprintf(stderr, "dct3: normalize not yet supported=%d\n",
+                    normalize);
+            break;
+    }
+}
 /**end repeat**/

Modified: trunk/scipy/fftpack/tests/test_real_transforms.py
===================================================================
--- trunk/scipy/fftpack/tests/test_real_transforms.py	2009-01-18 11:24:18 UTC (rev 5482)
+++ trunk/scipy/fftpack/tests/test_real_transforms.py	2009-01-18 11:24:37 UTC (rev 5483)
@@ -6,7 +6,7 @@
 from numpy.testing import assert_array_almost_equal, TestCase
 
 from scipy.io import loadmat
-from scipy.fftpack.realtransforms import dct1, dct2
+from scipy.fftpack.realtransforms import dct1, dct2, dct3
 
 TDATA = loadmat(join(dirname(__file__), 'test.mat'),
                 squeeze_me=True,  struct_as_record=True, mat_dtype=True)
@@ -153,5 +153,34 @@
     def setUp(self):
         self.rdt = np.double
 
+class _TestDCTIIIBase(TestCase):
+    def setUp(self):
+        self.rdt = None
+
+    def test_definition(self):
+        for i in range(len(X)):
+            x = np.array(X[i], dtype=self.rdt)
+            y = dct3(x)
+            self.failUnless(y.dtype == self.rdt,
+                    "Output dtype is %s, expected %s" % (y.dtype, self.rdt))
+            assert_array_almost_equal(dct2(y) / (2*x.size), x)
+
+    def test_axis(self):
+        nt = 2
+        for i in [7, 8, 9, 16, 32, 64]:
+            x = np.random.randn(nt, i)
+            y = dct3(x)
+            for j in range(nt):
+                assert_array_almost_equal(y[j], dct3(x[j]))
+
+            x = x.T
+            y = dct3(x, axis=0)
+            for j in range(nt):
+                assert_array_almost_equal(y[:,j], dct3(x[:,j]))
+
+class TestDCTIIIDouble(_TestDCTIIIBase):
+    def setUp(self):
+        self.rdt = np.double
+
 if __name__ == "__main__":
     np.testing.run_module_suite()



More information about the Scipy-svn mailing list