[Scipy-svn] r5482 - in trunk/scipy/fftpack: . src
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Jan 18 06:24:26 EST 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 at costi, @PREF at COSTI)(int*, @type@*);
+extern void F_FUNC(@pref at cost, @PREF at COST)(int*, @type@*, @type@*);
+extern void F_FUNC(@pref at cosqi, @PREF at COSQI)(int*, @type@*);
+extern void F_FUNC(@pref at cosqb, @PREF at COSQB)(int*, @type@*, @type@*);
+extern void F_FUNC(@pref at cosqf, @PREF at COSQF)(int*, @type@*, @type@*);
+
+GEN_CACHE(@pref at dct1,(int n)
+ , at type@* wsave;
+ ,(caches_ at pref@dct1[i].n==n)
+ ,caches_ at pref@dct1[id].wsave = malloc(sizeof(@type@)*(3*n+15));
+ F_FUNC(@pref at costi, @PREF at COSTI)(&n, caches_ at pref@dct1[id].wsave);
+ ,free(caches_ at pref@dct1[id].wsave);
+ ,10)
+
+GEN_CACHE(@pref at dct2,(int n)
+ , at type@* wsave;
+ ,(caches_ at pref@dct2[i].n==n)
+ ,caches_ at pref@dct2[id].wsave = malloc(sizeof(@type@)*(3*n+15));
+ F_FUNC(@pref at cosqi, at PREF@COSQI)(&n,caches_ at pref@dct2[id].wsave);
+ ,free(caches_ at pref@dct2[id].wsave);
+ ,10)
+
+void @pref at dct1(@type@ * inout, int n, int howmany, int normalize)
+{
+ int i;
+ @type@ *ptr = inout;
+ @type@ *wsave = NULL;
+
+ wsave = caches_ at pref@dct1[get_cache_id_ at pref@dct1(n)].wsave;
+
+ for (i = 0; i < howmany; ++i, ptr += n) {
+ @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;
+ }
+ }
+}
+
+void @pref at dct2(@type@ * inout, int n, int howmany, int normalize)
+{
+ int i, j;
+ @type@ *ptr = inout;
+ @type@ *wsave = NULL;
+ @type@ n1, n2;
+
+ wsave = caches_ at pref@dct2[get_cache_id_ at pref@dct2(n)].wsave;
+
+ for (i = 0; i < howmany; ++i, ptr += n) {
+ @pref at 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