From b535b43980a30b02ae367dc1fb58b9ac978e0025 Mon Sep 17 00:00:00 2001 From: dwarning Date: Mon, 2 Dec 2013 22:16:22 +0100 Subject: [PATCH] 1-f-code.c, use fftw3 --- src/frontend/trannoise/1-f-code.c | 73 +++++++++++++++++++++++++--- src/spicelib/devices/isrc/isrcacct.c | 6 ++- src/spicelib/devices/vsrc/vsrcacct.c | 6 ++- 3 files changed, 73 insertions(+), 12 deletions(-) diff --git a/src/frontend/trannoise/1-f-code.c b/src/frontend/trannoise/1-f-code.c index 9db935714..f00f44c90 100644 --- a/src/frontend/trannoise/1-f-code.c +++ b/src/frontend/trannoise/1-f-code.c @@ -15,6 +15,9 @@ #include "ngspice/fftext.h" #include "ngspice/wallace.h" +#ifdef HAVE_LIBFFTW3 +#include "fftw3.h" +#endif extern void controlled_exit(int status); @@ -22,15 +25,27 @@ extern void controlled_exit(int status); void f_alpha(int n_pts, int n_exp, double X[], double Q_d, double alpha) { - int i; - double *hfa, *wfa; + int i, length; double ha; + double *hfa, *wfa; - ha = alpha/2.0f; +#ifdef HAVE_LIBFFTW3 + fftw_complex *out = NULL; + fftw_plan plan_forward = NULL; + fftw_plan plan_backward = NULL; +#endif + + ha = alpha/2.0; // Q_d = sqrt(Q_d); /* find the deviation of the noise */ - hfa = TMALLOC(double,n_pts); - wfa = TMALLOC(double,n_pts); - hfa[0] = 1.0f; +#ifdef HAVE_LIBFFTW3 + length = 2 * (n_pts/2 + 1); +#else + length = n_pts; +#endif + hfa = TMALLOC(double, length); + wfa = TMALLOC(double, length); + + hfa[0] = 1.0; wfa[0] = Q_d * GaussWa; /* generate the coefficients hk */ for (i = 1; i < n_pts; i++) { @@ -40,8 +55,42 @@ f_alpha(int n_pts, int n_exp, double X[], double Q_d, double alpha) wfa[i] = Q_d * GaussWa; } - // for (i = 0; i < n_pts; i++) - // printf("rnd %e, hk %e\n", wfa[i], hfa[i]); +#ifdef HAVE_LIBFFTW3 + + /* in-place transformation needs zero padding on the end */ + hfa[n_pts] = 0.0; + wfa[n_pts] = 0.0; + hfa[n_pts+1] = 0.0; + wfa[n_pts+1] = 0.0; + + /* perform the discrete Fourier transform */ + plan_forward = fftw_plan_dft_r2c_1d(n_pts, hfa, (fftw_complex *)hfa, FFTW_ESTIMATE); + fftw_execute(plan_forward); + fftw_destroy_plan(plan_forward); + + plan_forward = fftw_plan_dft_r2c_1d(n_pts, wfa, (fftw_complex *)wfa, FFTW_ESTIMATE); + fftw_execute(plan_forward); + fftw_destroy_plan(plan_forward); + + out = fftw_malloc(sizeof(fftw_complex) * (unsigned int) (n_pts/2 + 1)); + /* multiply the two complex vectors */ + for (i = 0; i < n_pts/2 + 1; i++) { + out[i][0] = hfa[i]*wfa[i] - hfa[i+1]*wfa[i+1]; + out[i][1] = hfa[i]*wfa[i+1] + hfa[i+1]*wfa[i]; + } + + /* inverse transform */ + plan_backward = fftw_plan_dft_c2r_1d(n_pts, out, X, FFTW_ESTIMATE); + fftw_execute(plan_backward); + fftw_destroy_plan(plan_backward); + + for (i = 0; i < n_pts; i++) { + X[i] = X[i] / (double) n_pts; + } + + fftw_free(out); + +#else /* Green's FFT */ /* perform the discrete Fourier transform */ fftInit(n_exp); @@ -50,9 +99,12 @@ f_alpha(int n_pts, int n_exp, double X[], double Q_d, double alpha) /* multiply the two complex vectors */ rspectprod(hfa, wfa, X, n_pts); + /* inverse transform */ riffts(X, n_exp, 1); +#endif + free(hfa); free(wfa); /* fft tables will be freed in vsrcaccept.c and isrcaccept.c @@ -80,11 +132,16 @@ trnoise_state_gen(struct trnoise_state *this, CKTcircuit *ckt) size_t newsteps = 1; int newexp = 0; +#ifdef HAVE_LIBFFTW3 + newsteps = nosteps; + newexp = 1; +#else // generate number of steps as power of 2 while (newsteps < nosteps) { newsteps <<= 1; newexp++; } +#endif this->oneof = TMALLOC(double, newsteps); this->oneof_length = newsteps; diff --git a/src/spicelib/devices/isrc/isrcacct.c b/src/spicelib/devices/isrc/isrcacct.c index 33f656d6b..c56617569 100644 --- a/src/spicelib/devices/isrc/isrcacct.c +++ b/src/spicelib/devices/isrc/isrcacct.c @@ -12,9 +12,9 @@ Author: 1985 Thomas L. Quarles #include "ngspice/missing_math.h" #include "ngspice/1-f-code.h" -extern int fftInit(long M); +#ifndef HAVE_LIBFFTW3 extern void fftFree(void); -extern void rffts(float *data, long M, long Rows); +#endif extern bool ft_ngdebug; /* some additional debug info printed */ @@ -197,12 +197,14 @@ ISRCaccept(CKTcircuit *ckt, GENmodel *inModel) if ((TS == 0.0) && (RTSAM == 0.0)) // no further breakpoint if value not given break; +#ifndef HAVE_LIBFFTW3 /* FIXME, dont' want this here, over to aof_get or somesuch */ if (ckt->CKTtime == 0.0) { if (ft_ngdebug) printf("VSRC: free fft tables\n"); fftFree(); } +#endif if(ckt->CKTbreak) { diff --git a/src/spicelib/devices/vsrc/vsrcacct.c b/src/spicelib/devices/vsrc/vsrcacct.c index a4c40af6f..d1aacacc6 100644 --- a/src/spicelib/devices/vsrc/vsrcacct.c +++ b/src/spicelib/devices/vsrc/vsrcacct.c @@ -12,9 +12,9 @@ Author: 1985 Thomas L. Quarles #include "ngspice/missing_math.h" #include "ngspice/1-f-code.h" -extern int fftInit(long M); +#ifndef HAVE_LIBFFTW3 extern void fftFree(void); -extern void rffts(float *data, long M, long Rows); +#endif extern bool ft_ngdebug; /* some additional debug info printed */ @@ -197,12 +197,14 @@ VSRCaccept(CKTcircuit *ckt, GENmodel *inModel) if ((TS == 0.0) && (RTSAM == 0.0)) // no further breakpoint if value not given break; +#ifndef HAVE_LIBFFTW3 /* FIXME, dont' want this here, over to aof_get or somesuch */ if (ckt->CKTtime == 0.0) { if(ft_ngdebug) printf("VSRC: free fft tables\n"); fftFree(); } +#endif if(ckt->CKTbreak) {