[Scipy-svn] r5473 - trunk/scipy/fftpack/src

[scipy-svn at scipy.org]

Author: cdavid
Date: 2009-01-17 05:08:51 -0600 (Sat, 17 Jan 2009)
New Revision: 5473

Modified:
   trunk/scipy/fftpack/src/dct.c
Log:
Handle orthonormal mode for normalization (DCT matrix orthonormal).

Modified: trunk/scipy/fftpack/src/dct.c
===================================================================
--- trunk/scipy/fftpack/src/dct.c	2009-01-17 05:16:25 UTC (rev 5472)
+++ trunk/scipy/fftpack/src/dct.c	2009-01-17 11:08:51 UTC (rev 5473)
@@ -3,9 +3,15 @@
   Double complex FFT and IFFT.
   Author: Pearu Peterson, August 2002
  */
+#include <math.h>
 
 #include "fftpack.h"
 
+enum normalize {
+	DCT_NORMALIZE_NO = 0,
+	DCT_NORMALIZE_ORTHONORMAL = 1
+};
+
 extern void F_FUNC(dcosti,DCOSTI)(int*,double*);
 extern void F_FUNC(dcost,DCOST)(int*,double*,double*);
 extern void F_FUNC(dcosqi,DCOSQI)(int*,double*);
@@ -45,6 +51,9 @@
                                 normalize);
 	} else {
                 ptr = inout;
+		/* 0.5 coeff comes from fftpack defining DCT as
+		 * 4 * sum(cos(something)), whereas most definition 
+		 * use 2 */
                 for (i = n * howmany - 1; i >= 0; --i, ++ptr) {
                         *((double *) (ptr)) *= 0.5;
                 }
@@ -53,9 +62,10 @@
 
 void dct2(double * inout, int n, int howmany, int normalize)
 {
-	int i;
+	int i, j;
 	double *ptr = inout;
 	double *wsave = NULL;
+	double n1, n2;
 
 	wsave = caches_dct2[get_cache_id_dct2(n)].wsave;
 
@@ -64,13 +74,33 @@
 
         }
 
-	if (normalize) {
-                fprintf(stderr, "dct2: normalize not yet supported=%d\n",
-                                normalize);
-	} else {
-                ptr = inout;
-                for (i = n * howmany - 1; i >= 0; --i, ++ptr) {
-                        *((double *) (ptr)) *= 0.5;
-                }
+	switch (normalize) {
+	case DCT_NORMALIZE_NO:
+        ptr = inout;
+		/* 0.5 coeff comes from fftpack defining DCT as
+		 * 4 * sum(cos(something)), whereas most definition 
+		 * use 2 */
+        for (i = n * howmany - 1; i >= 0; --i, ++ptr) {
+            *((double *) (ptr)) *= 0.5;
         }
+		break;
+	case DCT_NORMALIZE_ORTHONORMAL:
+        ptr = inout;
+		/* 0.5 coeff comes from fftpack defining DCT as
+		 * 4 * sum(cos(something)), whereas most definition 
+		 * use 2 */
+		n1 = 0.25 * sqrt(1./n);
+		n2 = 0.25 * sqrt(2./n);
+        for (i = 0; i < howmany; ++i, ptr+=n) {
+            ptr[0] *= n1;
+            for (j = 1; j < n; ++j) {
+                ptr[j] *= n2;
+            }
+        } 
+		break;
+	default:
+        fprintf(stderr, "dct2: normalize not yet supported=%d\n",
+                        normalize);
+		break;
+    }
 }