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

[scipy-svn at scipy.org]

Author: cdavid
Date: 2009-01-19 02:31:02 -0600 (Mon, 19 Jan 2009)
New Revision: 5492

Modified:
   trunk/scipy/fftpack/realtransforms.py
   trunk/scipy/fftpack/src/dct.c.src
Log:
Implement dct 1 normalization mode similarly to other types.

Modified: trunk/scipy/fftpack/realtransforms.py
===================================================================
--- trunk/scipy/fftpack/realtransforms.py	2009-01-18 11:27:04 UTC (rev 5491)
+++ trunk/scipy/fftpack/realtransforms.py	2009-01-19 08:31:02 UTC (rev 5492)
@@ -11,7 +11,7 @@
 atexit.register(_fftpack.destroy_ddct1_cache)
 atexit.register(_fftpack.destroy_ddct2_cache)
 
-def dct1(x, n=None, axis=-1):
+def dct1(x, n=None, axis=-1, norm=None):
     """
     Return the Discrete Cosine Transform (type I) of arbitrary type sequence x.
 
@@ -28,7 +28,7 @@
     -------
     y : real ndarray
     """
-    return _dct(x, 1, n, axis)
+    return _dct(x, 1, n, axis, normalize=norm)
 
 def dct2(x, n=None, axis=-1, norm=None):
     """

Modified: trunk/scipy/fftpack/src/dct.c.src
===================================================================
--- trunk/scipy/fftpack/src/dct.c.src	2009-01-18 11:27:04 UTC (rev 5491)
+++ trunk/scipy/fftpack/src/dct.c.src	2009-01-19 08:31:02 UTC (rev 5492)
@@ -42,8 +42,8 @@
 
 void @pref at dct1(@type@ * inout, int n, int howmany, int normalize)
 {
-    int i;
-    @type@ *ptr = inout;
+    int i, j;
+    @type@ *ptr = inout, n1, n2;
     @type@ *wsave = NULL;
 
     wsave = caches_ at pref@dct1[get_cache_id_ at pref@dct1(n)].wsave;
@@ -52,17 +52,26 @@
         @pref at cost_(&n, ptr, wsave);
     }
 
-    if (normalize) {
-        fprintf(stderr, "dct1: normalize not yet supported=%d\n",
-                normalize);
-    } else {
-        ptr = inout;
-        /* 0.5 coeff comes from fftpack defining DCT as
-         * 4 * sum(cos(something)), whereas most definition
-         * use 2 */
-        for (i = 0; i < n * howmany; ++i) {
-            ptr[i] *= 0.5;
-        }
+    switch (normalize) {
+        case DCT_NORMALIZE_NO:
+            break;
+#if 0
+        case DCT_NORMALIZE_ORTHONORMAL:
+            ptr = inout;
+            n1 = sqrt(0.5 / (n-1));
+            n2 = sqrt(1. / (n-1));
+            for (i = 0; i < howmany; ++i, ptr+=n) {
+                ptr[0] *= n1;
+                for (j = 1; j < n-1; ++j) {
+                    ptr[j] *= n2;
+                }
+            }
+            break;
+#endif
+        default:
+            fprintf(stderr, "dct1: normalize not yet supported=%d\n",
+                    normalize);
+            break;
     }
 }