diff --git a/src/frontend/com_fft.c b/src/frontend/com_fft.c index ced012d4e..3e001ba94 100644 --- a/src/frontend/com_fft.c +++ b/src/frontend/com_fft.c @@ -14,7 +14,6 @@ Author: 2008 Holger Vogt #include "com_fft.h" #include "variable.h" -#include "missing_math.h" void com_fft(wordlist *wl) @@ -27,106 +26,111 @@ com_fft(wordlist *wl) struct dvec *f, *vlist, *lv = NULL, *vec; struct pnode *names, *first_name; - float *reald, *imagd; - int size, sign, isreal; - float scaling; + double *reald, *imagd; + int size, sign, order; + double scale, sigma; if (!plot_cur || !plot_cur->pl_scale) { fprintf(cp_err, "Error: no vectors loaded.\n"); return; } - if (!isreal(plot_cur->pl_scale) || + 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 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; + /* output vector has length of size/2 */ + fpts = size/2; - // window function + /* window functions - should have an average of one */ win = (double *) tmalloc(tlen * sizeof (double)); { char window[BSIZE_SP]; double maxt = time[tlen-1]; - if (!cp_getvar("specwindow", VT_STRING, window)) + if (!cp_getvar("specwindow", VT_STRING, window)) strcpy(window,"blackman"); if (eq(window, "none")) for(i=0; i span) { - win[i] = 0; + win[i] = 0.0; } else { - win[i] = 1; + win[i] = 1.0; } } - else if (eq(window, "hanning") || eq(window, "cosine")) + else if (eq(window, "triangle") || eq(window, "bartlet") || eq(window, "bartlett")) for(i=0; i span) { - win[i] = 0; + win[i] = 0.0; } else { - win[i] = 1 - cos(2*M_PI*(time[i]-maxt)/span); + win[i] = 2.0 - fabs(2+4*(time[i]-maxt)/span); } } - else if (eq(window, "hamming")) + else if (eq(window, "hann") || eq(window, "hanning") || eq(window, "cosine")) for(i=0; i span) { - win[i] = 0; + win[i] = 0.0; } else { - win[i] = 1 - 0.92/1.08*cos(2*M_PI*(time[i]-maxt)/span); + win[i] = 1.0 - cos(2*M_PI*(time[i]-maxt)/span); } } - else if (eq(window, "triangle") || eq(window, "bartlet")) + else if (eq(window, "hamming")) for(i=0; i span) { - win[i] = 0; + win[i] = 0.0; } else { - win[i] = 2 - fabs(2+4*(time[i]-maxt)/span); + win[i] = 1.0 - 0.46/0.54*cos(2*M_PI*(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 */ + else if (eq(window, "blackman")) for(i=0; i span) { win[i] = 0; } else { - win[i] = 1; + win[i] = 1.0; 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(); + else if (eq(window, "flattop")) + for(i=0; i span) { + win[i] = 0; + } else { + win[i] = 1.0; + win[i] -= 1.93*cos(2*M_PI*(time[i]-maxt)/span); + win[i] += 1.29*cos(4*M_PI*(time[i]-maxt)/span); + win[i] -= 0.388*cos(6*M_PI*(time[i]-maxt)/span); + win[i] += 0.032*cos(8*M_PI*(time[i]-maxt)/span); + } + } + else if (eq(window, "gaussian")) { 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))); + sigma=1.0/order; + scale=0.83/sigma; 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; + win[i] = scale*exp(-0.5*pow((time[i]-maxt/2)/(sigma*maxt/2),2)); } } - } else { + } else { fprintf(cp_err, "Warning: unknown window type %s\n", window); tfree(win); return; @@ -148,7 +152,7 @@ com_fft(wordlist *wl) continue; } if (!isreal(vec)) { - fprintf(cp_err, "Error: %s isn't real!\n", + fprintf(cp_err, "Error: %s isn't real!\n", vec->v_name); vec = vec->v_link2; continue; @@ -166,12 +170,12 @@ com_fft(wordlist *wl) ngood++; } } - free_pnode_o(first_name); /* h_vogt 081206 */ + free_pnode_o(first_name); if (!ngood) { tfree(win); return; } - + plot_cur = plot_alloc("spectrum"); plot_cur->pl_next = plot_list; plot_list = plot_cur; @@ -189,7 +193,7 @@ com_fft(wordlist *wl) f->v_realdata = freq; vec_new(f); - for (i = 0; iv_name = vec_basename(vec); - f->v_type = SV_NOTYPE; //vec->v_type; + f->v_type = SV_NOTYPE; f->v_flags = (VF_COMPLEX | VF_PERMANENT); f->v_length = fpts; f->v_compdata = fdvec[i]; @@ -208,63 +212,61 @@ com_fft(wordlist *wl) vec = vec->v_link2; } - sign = 1; - isreal = 1; - - reald = (float*)tmalloc(size*sizeof(float)); - imagd = (float*)tmalloc(size*sizeof(float)); - - printf("CPU: Delta Freq %f Hz, input length %d, output length %d\n", 1./span*tlen/size, size, fpts); - + + printf("FFT: Time span: %g s, input length: %d, zero padding: %d\n", span, size, size-tlen); + printf("FFT: Freq. resolution: %g Hz, output length: %d\n", 1.0/span*tlen/size, fpts); + + reald = (double*)tmalloc(size*sizeof(double)); + imagd = (double*)tmalloc(size*sizeof(double)); for (i = 0; i