From 655d8aea02a86f13340200f0e49b58ff55b6436e Mon Sep 17 00:00:00 2001 From: h_vogt Date: Sat, 24 May 2008 18:06:33 +0000 Subject: [PATCH] Fast fourier transform for transient data analysis --- .cvsignore | 3 + ChangeLog | 3 + src/frontend/Makefile.am | 2 + src/frontend/com_fft.c | 328 +++++++++++++++++++++++++++++++++++++++ src/frontend/com_fft.h | 13 ++ src/frontend/commands.c | 4 + src/include/fteext.h | 5 +- 7 files changed, 356 insertions(+), 2 deletions(-) create mode 100644 src/frontend/com_fft.c create mode 100644 src/frontend/com_fft.h diff --git a/.cvsignore b/.cvsignore index 52af2a9f9..6c8606a9e 100644 --- a/.cvsignore +++ b/.cvsignore @@ -22,3 +22,6 @@ ylwrap *.tar.gz +visualc +autogen-new.log +make.log \ No newline at end of file diff --git a/ChangeLog b/ChangeLog index 2341e0da3..95ed802f4 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,6 @@ +2008-05-24 Holger Vogt + * src/frontend/com_fft.c: Fast fourier transform added for transient data analysis. + 2008-05-18 Dietmar Warning * Small changes to compile under Sun Studio 11 for Solaris - may be useful in other configurations too diff --git a/src/frontend/Makefile.am b/src/frontend/Makefile.am index e4ba64f98..0cb67396c 100644 --- a/src/frontend/Makefile.am +++ b/src/frontend/Makefile.am @@ -55,6 +55,8 @@ libfte_a_SOURCES = \ com_shell.h \ com_shift.c \ com_shift.h \ + com_fft.c \ + com_fft.h \ com_state.c \ com_state.h \ com_strcmp.c \ diff --git a/src/frontend/com_fft.c b/src/frontend/com_fft.c new file mode 100644 index 000000000..4cd2bca60 --- /dev/null +++ b/src/frontend/com_fft.c @@ -0,0 +1,328 @@ +/********** +Copyright 2008 Holger Vogt. All rights reserved. +Author: 2008 Holger Vogt +**********/ + +/* + * Code to do fast fourier transform on data. + */ + +#include "ngspice.h" +#include "ftedefs.h" +#include "dvec.h" +#include "sim.h" + +#include "com_fft.h" +#include "variable.h" +#include "missing_math.h" + +void +com_fft(wordlist *wl) +{ + complex **fdvec; + double **tdvec; + double *freq, *win, *time; + double delta_t, span; + int fpts, i, j, tlen, ngood; + struct dvec *f, *vlist, *lv, *vec; + struct pnode *names, *first_name; + + float *reald, *imagd; + int size, sign, isreal; + float scaling; + + if (!plot_cur || !plot_cur->pl_scale) { + fprintf(cp_err, "Error: no vectors loaded.\n"); + return; + } + if (!isreal(plot_cur->pl_scale) || + ((plot_cur->pl_scale)->v_type != SV_TIME)) { + fprintf(cp_err, "Error: fft needs real time scale\n"); + return; + } + + + tlen = (plot_cur->pl_scale)->v_length; + time = (plot_cur->pl_scale)->v_realdata; + span = time[tlen-1] - time[0]; + delta_t = span/(tlen - 1); + + // size of input vector is power of two and larger than spice vector + size = 1; + while (size < tlen) + size *= 2; + + // output vector has length of size/2 + fpts = size/2; + + // window function + win = (double *) tmalloc(tlen * sizeof (double)); + { + char window[BSIZE_SP]; + double maxt = time[tlen-1]; + if (!cp_getvar("specwindow", VT_STRING, window)) + strcpy(window,"blackman"); + if (eq(window, "none")) + for(i=0; i span) { + win[i] = 0; + } else { + win[i] = 1; + } + } + else if (eq(window, "hanning") || eq(window, "cosine")) + for(i=0; i span) { + win[i] = 0; + } else { + win[i] = 1 - cos(2*M_PI*(time[i]-maxt)/span); + } + } + else if (eq(window, "hamming")) + for(i=0; i span) { + win[i] = 0; + } else { + win[i] = 1 - 0.92/1.08*cos(2*M_PI*(time[i]-maxt)/span); + } + } + else if (eq(window, "triangle") || eq(window, "bartlet")) + for(i=0; i span) { + win[i] = 0; + } else { + win[i] = 2 - fabs(2+4*(time[i]-maxt)/span); + } + } + else if (eq(window, "blackman")) { + int order; + if (!cp_getvar("specwindoworder", VT_NUM, &order)) order = 2; + if (order < 2) order = 2; /* only order 2 supported here */ + for(i=0; i span) { + win[i] = 0; + } else { + win[i] = 1; + win[i] -= 0.50/0.42*cos(2*M_PI*(time[i]-maxt)/span); + win[i] += 0.08/0.42*cos(4*M_PI*(time[i]-maxt)/span); + } + } + } else if (eq(window, "gaussian")) { + int order; + double scale; + extern double erfc(); + if (!cp_getvar("specwindoworder", VT_NUM, &order)) order = 2; + if (order < 2) order = 2; + scale = pow(2*M_PI/order,0.5)*(0.5-erfc(pow(order,0.5))); + for(i=0; i span) { + win[i] = 0; + } else { + win[i] = exp(-0.5*order*(1-2*(maxt-time[i])/span) + *(1-2*(maxt-time[i])/span))/scale; + } + } + } else { + fprintf(cp_err, "Warning: unknown window type %s\n", window); + tfree(win); + return; + } + } + + names = ft_getpnames(wl, TRUE); + first_name = names; + vlist = NULL; + ngood = 0; + while (names) { + vec = ft_evaluate(names); + names = names->pn_next; + while (vec) { + if (vec->v_length != tlen) { + fprintf(cp_err, "Error: lengths of %s vectors don't match: %d, %d\n", + vec->v_name, vec->v_length, tlen); + vec = vec->v_link2; + continue; + } + if (!isreal(vec)) { + fprintf(cp_err, "Error: %s isn't real!\n", + vec->v_name); + vec = vec->v_link2; + continue; + } + if (vec->v_type == SV_TIME) { + vec = vec->v_link2; + continue; + } + if (!vlist) + vlist = vec; + else + lv->v_link2 = vec; + lv = vec; + vec = vec->v_link2; + ngood++; + } + } + free_pnode(first_name); + if (!ngood) { + return; + } + + plot_cur = plot_alloc("spectrum"); + plot_cur->pl_next = plot_list; + plot_list = plot_cur; + plot_cur->pl_title = copy((plot_cur->pl_next)->pl_title); + plot_cur->pl_name = copy("Spectrum"); + plot_cur->pl_date = copy(datestring( )); + + freq = (double *) tmalloc(fpts * sizeof(double)); + f = alloc(struct dvec); + ZERO(f, struct dvec); + f->v_name = copy("frequency"); + f->v_type = SV_FREQUENCY; + f->v_flags = (VF_REAL | VF_PERMANENT | VF_PRINT); + f->v_length = fpts; + f->v_realdata = freq; + vec_new(f); + + for (i = 0; iv_realdata; /* real input data */ + fdvec[i] = (complex *) tmalloc(fpts * sizeof(complex)); /* complex output data */ + f = alloc(struct dvec); + ZERO(f, struct dvec); + f->v_name = vec_basename(vec); + f->v_type = SV_NOTYPE; //vec->v_type; + f->v_flags = (VF_COMPLEX | VF_PERMANENT); + f->v_length = fpts; + f->v_compdata = fdvec[i]; + vec_new(f); + vec = vec->v_link2; + } + + + sign = 1; + isreal = 1; + + reald = (float*)malloc(size*sizeof(float)); + imagd = (float*)malloc(size*sizeof(float)); + + printf("CPU: Delta Freq %f Hz, input length %d, output length %d\n", 1./span*tlen/size, size, fpts); + + for (i = 0; i> 1; + j = 0; + for (i=0;i>= 1; + } + j += k; + } + + /* Compute the FFT */ + c1 = -1.0; + c2 = 0.0; + l2 = 1; + for (l=0;l