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

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


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

Added:
   trunk/scipy/fftpack/src/dct.c.src
Modified:
   trunk/scipy/fftpack/fftpack.pyf
   trunk/scipy/fftpack/realtransforms.py
   trunk/scipy/fftpack/setup.py
Log:
Use template for dct sources, to support multiple types.

Modified: trunk/scipy/fftpack/fftpack.pyf
===================================================================
--- trunk/scipy/fftpack/fftpack.pyf	2009-01-18 11:23:59 UTC (rev 5481)
+++ trunk/scipy/fftpack/fftpack.pyf	2009-01-18 11:24:18 UTC (rev 5482)
@@ -163,35 +163,35 @@
          intent(c) destroy_rfft_cache
        end subroutine destroy_rfft_cache
 
-       subroutine dct1(x,n,howmany,normalize)
-         ! y = dct1(x[,n,normalize,overwrite_x])
-         intent(c) dct1
+       subroutine ddct1(x,n,howmany,normalize)
+         ! y = ddct1(x[,n,normalize,overwrite_x])
+         intent(c) ddct1
          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 dct1
+       end subroutine ddct1
 
-       subroutine dct2(x,n,howmany,normalize)
-         ! y = dct2(x[,n,normalize,overwrite_x])
-         intent(c) dct2
+       subroutine ddct2(x,n,howmany,normalize)
+         ! y = ddct2(x[,n,normalize,overwrite_x])
+         intent(c) ddct2
          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 dct2
+       end subroutine ddct2
 
-       subroutine destroy_dct2_cache()
-         intent(c) destroy_dct2_cache
-       end subroutine destroy_dct2_cache
+       subroutine destroy_ddct2_cache()
+         intent(c) destroy_ddct2_cache
+       end subroutine destroy_ddct2_cache
 
-       subroutine destroy_dct1_cache()
-         intent(c) destroy_dct1_cache
-       end subroutine destroy_dct1_cache
+       subroutine destroy_ddct1_cache()
+         intent(c) destroy_ddct1_cache
+       end subroutine destroy_ddct1_cache
 
     end interface 
 end python module _fftpack

Modified: trunk/scipy/fftpack/realtransforms.py
===================================================================
--- trunk/scipy/fftpack/realtransforms.py	2009-01-18 11:23:59 UTC (rev 5481)
+++ trunk/scipy/fftpack/realtransforms.py	2009-01-18 11:24:18 UTC (rev 5482)
@@ -8,8 +8,8 @@
 from scipy.fftpack import _fftpack
 
 import atexit
-atexit.register(_fftpack.destroy_dct1_cache)
-atexit.register(_fftpack.destroy_dct2_cache)
+atexit.register(_fftpack.destroy_ddct1_cache)
+atexit.register(_fftpack.destroy_ddct2_cache)
 
 def dct1(x, n=None, axis=-1):
     """
@@ -94,9 +94,9 @@
         raise NotImplemented("Padding/truncating not yet implemented")
 
     if type == 1:
-        f = _fftpack.dct1
+        f = _fftpack.ddct1
     elif type == 2:
-        f = _fftpack.dct2
+        f = _fftpack.ddct2
     else:
         raise ValueError("Type %d not understood" % type)
 

Modified: trunk/scipy/fftpack/setup.py
===================================================================
--- trunk/scipy/fftpack/setup.py	2009-01-18 11:23:59 UTC (rev 5481)
+++ trunk/scipy/fftpack/setup.py	2009-01-18 11:24:18 UTC (rev 5482)
@@ -18,11 +18,12 @@
                        sources=[join('src/fftpack','*.f')])
 
     sources = ['fftpack.pyf','src/zfft.c','src/drfft.c','src/zrfft.c',
-               'src/zfftnd.c', 'src/dct.c']
+               'src/zfftnd.c', 'src/dct.c.src']
 
     config.add_extension('_fftpack',
         sources=sources,
-        libraries=['dfftpack', 'fftpack'])
+        libraries=['dfftpack', 'fftpack'],
+        include_dirs=['src'])
 
     config.add_extension('convolve',
         sources=['convolve.pyf','src/convolve.c'],

Added: trunk/scipy/fftpack/src/dct.c.src
===================================================================
--- trunk/scipy/fftpack/src/dct.c.src	2009-01-18 11:23:59 UTC (rev 5481)
+++ trunk/scipy/fftpack/src/dct.c.src	2009-01-18 11:24:18 UTC (rev 5482)
@@ -0,0 +1,114 @@
+/* vim:syntax=c
+ * vim:sw=4
+ *
+ * Interfaces to the DCT transforms of fftpack
+ */
+#include <math.h>
+
+#include "fftpack.h"
+
+enum normalize {
+    DCT_NORMALIZE_NO = 0,
+    DCT_NORMALIZE_ORTHONORMAL = 1
+};
+
+/**begin repeat
+
+#type=float,double#
+#pref=,d#
+#PREF=,D#
+*/
+extern void F_FUNC(@pref@costi, @PREF@COSTI)(int*, @type@*);
+extern void F_FUNC(@pref@cost, @PREF@COST)(int*, @type@*, @type@*);
+extern void F_FUNC(@pref@cosqi, @PREF@COSQI)(int*, @type@*);
+extern void F_FUNC(@pref@cosqb, @PREF@COSQB)(int*, @type@*, @type@*);
+extern void F_FUNC(@pref@cosqf, @PREF@COSQF)(int*, @type@*, @type@*);
+
+GEN_CACHE(@pref@dct1,(int n)
+      ,@type@* wsave;
+      ,(caches_@pref@dct1[i].n==n)
+      ,caches_@pref@dct1[id].wsave = malloc(sizeof(@type@)*(3*n+15));
+       F_FUNC(@pref@costi, @PREF@COSTI)(&n, caches_@pref@dct1[id].wsave);
+      ,free(caches_@pref@dct1[id].wsave);
+      ,10)
+
+GEN_CACHE(@pref@dct2,(int n)
+      ,@type@* wsave;
+      ,(caches_@pref@dct2[i].n==n)
+      ,caches_@pref@dct2[id].wsave = malloc(sizeof(@type@)*(3*n+15));
+       F_FUNC(@pref@cosqi,@PREF@COSQI)(&n,caches_@pref@dct2[id].wsave);
+      ,free(caches_@pref@dct2[id].wsave);
+      ,10)
+
+void @pref@dct1(@type@ * inout, int n, int howmany, int normalize)
+{
+    int i;
+    @type@ *ptr = inout;
+    @type@ *wsave = NULL;
+
+    wsave = caches_@pref@dct1[get_cache_id_@pref@dct1(n)].wsave;
+
+    for (i = 0; i < howmany; ++i, ptr += n) {
+        @pref@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;
+        }
+    }
+}
+
+void @pref@dct2(@type@ * inout, int n, int howmany, int normalize)
+{
+    int i, j;
+    @type@ *ptr = inout;
+    @type@ *wsave = NULL;
+    @type@ n1, n2;
+
+    wsave = caches_@pref@dct2[get_cache_id_@pref@dct2(n)].wsave;
+
+    for (i = 0; i < howmany; ++i, ptr += n) {
+        @pref@cosqb_(&n, 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 = 0; i < n * howmany; ++i) {
+                ptr[i] *= 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;
+    }
+}
+
+/**end repeat**/



More information about the Scipy-svn mailing list