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

[scipy-svn at scipy.org]

Author: cdavid
Date: 2009-01-12 20:01:59 -0600 (Mon, 12 Jan 2009)
New Revision: 5452

Added:
   trunk/scipy/fftpack/src/dct.c
Modified:
   trunk/scipy/fftpack/SConscript
   trunk/scipy/fftpack/fftpack.pyf
   trunk/scipy/fftpack/setup.py
Log:
Start working on dct wrappers of fftpack.

Modified: trunk/scipy/fftpack/SConscript
===================================================================
--- trunk/scipy/fftpack/SConscript	2009-01-13 02:01:42 UTC (rev 5451)
+++ trunk/scipy/fftpack/SConscript	2009-01-13 02:01:59 UTC (rev 5452)
@@ -1,4 +1,4 @@
-# Last Change: Sat Nov 01 10:00 PM 2008 J
+# Last Change: Sat Jan 10 05:00 PM 2009 J
 # vim:syntax=python
 from os.path import join as pjoin
 
@@ -26,7 +26,8 @@
 env.PrependUnique(LIBPATH = ['.'])
 
 # Build _fftpack
-src = ['src/zfft.c','src/drfft.c','src/zrfft.c', 'src/zfftnd.c', 'fftpack.pyf']
+src = ['src/zfft.c','src/drfft.c','src/zrfft.c', 'src/zfftnd.c', 'src/dct.c',
+    'fftpack.pyf']
 env.NumpyPythonExtension('_fftpack', src)
 
 # Build convolve

Modified: trunk/scipy/fftpack/fftpack.pyf
===================================================================
--- trunk/scipy/fftpack/fftpack.pyf	2009-01-13 02:01:42 UTC (rev 5451)
+++ trunk/scipy/fftpack/fftpack.pyf	2009-01-13 02:01:59 UTC (rev 5452)
@@ -163,6 +163,36 @@
          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
+         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
+
+       subroutine dct2(x,n,howmany,normalize)
+         ! y = dct2(x[,n,normalize,overwrite_x])
+         intent(c) dct2
+         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
+
+       subroutine destroy_dct2_cache()
+         intent(c) destroy_dct2_cache
+       end subroutine destroy_dct2_cache
+
+       subroutine destroy_dct1_cache()
+         intent(c) destroy_dct1_cache
+       end subroutine destroy_dct1_cache
+
     end interface 
 end python module _fftpack
 

Modified: trunk/scipy/fftpack/setup.py
===================================================================
--- trunk/scipy/fftpack/setup.py	2009-01-13 02:01:42 UTC (rev 5451)
+++ trunk/scipy/fftpack/setup.py	2009-01-13 02:01:59 UTC (rev 5452)
@@ -18,7 +18,7 @@
                        sources=[join('src/fftpack','*.f')])
 
     sources = ['fftpack.pyf','src/zfft.c','src/drfft.c','src/zrfft.c',
-               'src/zfftnd.c']
+               'src/zfftnd.c', 'src/dct.c']
 
     config.add_extension('_fftpack',
         sources=sources,

Added: trunk/scipy/fftpack/src/dct.c
===================================================================
--- trunk/scipy/fftpack/src/dct.c	2009-01-13 02:01:42 UTC (rev 5451)
+++ trunk/scipy/fftpack/src/dct.c	2009-01-13 02:01:59 UTC (rev 5452)
@@ -0,0 +1,71 @@
+/*
+  Interface to various FFT libraries.
+  Double complex FFT and IFFT.
+  Author: Pearu Peterson, August 2002
+ */
+
+#include "fftpack.h"
+
+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*);
+extern void F_FUNC(dcosqb,DCOSQB)(int*,double*,double*);
+extern void F_FUNC(dcosqf,DCOSQF)(int*,double*,double*);
+
+GEN_CACHE(dct1,(int n)
+	  ,double* wsave;
+	  ,(caches_dct1[i].n==n)
+	  ,caches_dct1[id].wsave = (double*)malloc(sizeof(double)*(3*n+15));
+	   F_FUNC(dcosti,DCOSTI)(&n,caches_dct1[id].wsave);
+	  ,free(caches_dct1[id].wsave);
+	  ,10)
+
+GEN_CACHE(dct2,(int n)
+	  ,double* wsave;
+	  ,(caches_dct2[i].n==n)
+	  ,caches_dct2[id].wsave = (double*)malloc(sizeof(double)*(3*n+15));
+	   F_FUNC(dcosqi,DCOSQI)(&n,caches_dct2[id].wsave);
+	  ,free(caches_dct2[id].wsave);
+	  ,10)
+
+void dct1(double * inout, int n, int howmany, int normalize)
+{
+	int i;
+	double *ptr = inout;
+	double *wsave = NULL;
+
+	wsave = caches_dct1[get_cache_id_dct1(n)].wsave;
+
+        for (i = 0; i < howmany; ++i, ptr += n) {
+                dcost_(&n, (double*)(ptr), wsave);
+        }
+
+	if (normalize) {
+                fprintf(stderr, "dct1: normalize not yet supported=%d\n", 
+                                normalize);
+	} else {
+                ptr = inout;
+                for (i = n * howmany - 1; i >= 0; --i, ++ptr) {
+                        *((double *) (ptr)) *= 0.5;
+                }
+        }
+}
+
+void dct2(double * inout, int n, int howmany, int normalize)
+{
+	int i;
+	double *ptr = inout;
+	double *wsave = NULL;
+
+	wsave = caches_dct2[get_cache_id_dct2(n)].wsave;
+
+        for (i = 0; i < howmany; ++i, ptr += n) {
+                dcosqb_(&n, (double *) (ptr), wsave);
+
+        }
+
+	if (normalize) {
+                fprintf(stderr, "dct2: normalize not yet supported=%d\n", 
+                                normalize);
+	}
+}