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

[scipy-svn at scipy.org]

Author: cdavid
Date: 2009-01-17 05:09:41 -0600 (Sat, 17 Jan 2009)
New Revision: 5475

Modified:
   trunk/scipy/fftpack/realtransforms.py
   trunk/scipy/fftpack/src/dct.c
Log:
Add float versions of dct I and II.

Modified: trunk/scipy/fftpack/realtransforms.py
===================================================================
--- trunk/scipy/fftpack/realtransforms.py	2009-01-17 11:09:18 UTC (rev 5474)
+++ trunk/scipy/fftpack/realtransforms.py	2009-01-17 11:09:41 UTC (rev 5475)
@@ -30,7 +30,7 @@
     """
     return _dct(x, 1, n, axis)
 
-def dct2(x, n=None, axis=-1):
+def dct2(x, n=None, axis=-1, norm=None):
     """
     Return Discrete Cosine Transform (type II) of arbitrary type sequence x.
     There are several definitions, we use the following:
@@ -62,9 +62,9 @@
     'A Fast Cosine Transform in One and Two Dimensions', by J. Makhoul, in IEEE
     Transactions on acoustics, speech and signal processing.
     """
-    return _dct(x, 2, n, axis)
+    return _dct(x, 2, n, axis, normalize=norm)
 
-def _dct(x, type, n=None, axis=-1, overwrite_x=0):
+def _dct(x, type, n=None, axis=-1, overwrite_x=0, normalize=None):
     """
     Return Discrete Cosine Transform of arbitrary type sequence x.
 
@@ -100,11 +100,19 @@
     else:
         raise ValueError("Type %d not understood" % type)
 
+    if normalize:
+        if normalize == "ortho":
+            nm = 1
+        else:
+            raise ValueError("Unknown normalize mode %s" % normalize)
+    else:
+        nm = 0
+
     if axis == -1 or axis == len(tmp.shape) - 1:
-        return f(tmp, n, 0, overwrite_x)
+        return f(tmp, n, nm, overwrite_x)
     #else:
     #    raise NotImplementedError("Axis arg not yet implemented")
 
     tmp = np.swapaxes(tmp, axis, -1)
-    tmp = f(tmp, n, 0, overwrite_x)
+    tmp = f(tmp, n, nm, overwrite_x)
     return np.swapaxes(tmp, axis, -1)

Modified: trunk/scipy/fftpack/src/dct.c
===================================================================
--- trunk/scipy/fftpack/src/dct.c	2009-01-17 11:09:18 UTC (rev 5474)
+++ trunk/scipy/fftpack/src/dct.c	2009-01-17 11:09:41 UTC (rev 5475)
@@ -18,6 +18,12 @@
 extern void F_FUNC(dcosqb,DCOSQB)(int*,double*,double*);
 extern void F_FUNC(dcosqf,DCOSQF)(int*,double*,double*);
 
+extern void F_FUNC(costi,DCOSTI)(int*,float*);
+extern void F_FUNC(cost,COST)(int*,float*,float*);
+extern void F_FUNC(cosqi,COSQI)(int*,float*);
+extern void F_FUNC(cosqb,COSQB)(int*,float*,float*);
+extern void F_FUNC(cosqf,COSQF)(int*,float*,float*);
+
 GEN_CACHE(dct1,(int n)
 	  ,double* wsave;
 	  ,(caches_dct1[i].n==n)
@@ -34,6 +40,22 @@
 	  ,free(caches_dct2[id].wsave);
 	  ,10)
 
+GEN_CACHE(fdct1,(int n)
+	  ,float* wsave;
+	  ,(caches_fdct1[i].n==n)
+	  ,caches_fdct1[id].wsave = (float*)malloc(sizeof(float)*(3*n+15));
+	   F_FUNC(costi,COSTI)(&n,caches_fdct1[id].wsave);
+	  ,free(caches_fdct1[id].wsave);
+	  ,10)
+
+GEN_CACHE(fdct2,(int n)
+	  ,float* wsave;
+	  ,(caches_fdct2[i].n==n)
+	  ,caches_fdct2[id].wsave = (float*)malloc(sizeof(float)*(3*n+15));
+	   F_FUNC(cosqi,DCOSQI)(&n,caches_fdct2[id].wsave);
+	  ,free(caches_fdct2[id].wsave);
+	  ,10)
+
 void dct1(double * inout, int n, int howmany, int normalize)
 {
 	int i;
@@ -104,3 +126,74 @@
 		break;
     }
 }
+
+void fdct1(float * inout, int n, int howmany, int normalize)
+{
+	int i;
+	float *ptr = inout;
+	float *wsave = NULL;
+
+	wsave = caches_fdct1[get_cache_id_fdct1(n)].wsave;
+
+        for (i = 0; i < howmany; ++i, ptr += n) {
+                cost_(&n, (float*)(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 = n * howmany - 1; i >= 0; --i, ++ptr) {
+                        *((float *) (ptr)) *= 0.5;
+                }
+        }
+}
+
+void fdct2(float * inout, int n, int howmany, int normalize)
+{
+	int i, j;
+	float *ptr = inout;
+	float *wsave = NULL;
+	float n1, n2;
+
+	wsave = caches_fdct2[get_cache_id_fdct2(n)].wsave;
+
+        for (i = 0; i < howmany; ++i, ptr += n) {
+                cosqb_(&n, (float *) (ptr), wsave);
+
+        }
+
+	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) {
+            *((float *) (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;
+    }
+}