luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

move the `fft_windows()' function to fft/fftext.c

This commit is contained in:
dwarning 2013-08-03 17:35:01 +02:00 committed by rlar
parent c972744ff1
commit f3e7bf141d
3 changed files with 99 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

91
src/frontend/com_fft.c
View File

@ -36,6 +36,8 @@ com_fft(wordlist *wl)
int fpts, i, j, tlen, ngood;
struct dvec *f, *vlist, *lv = NULL, *vec;
struct pnode *pn, *names = NULL;
char window[BSIZE_SP];
double maxt;
#ifdef GREEN
int mm;
@ -43,7 +45,7 @@ com_fft(wordlist *wl)
double *reald = NULL, *imagd = NULL;
int size, sign, order;
double scale, sigma;
double scale;
if (!plot_cur || !plot_cur->pl_scale) {
fprintf(cp_err, "Error: no vectors loaded.\n");
@ -76,84 +78,17 @@ com_fft(wordlist *wl)
/* output vector has length of size/2 */
fpts = size/2;
/* window functions - should have an average of one */
win = TMALLOC(double, tlen);
{
char window[BSIZE_SP];
double maxt = time[tlen-1];
if (!cp_getvar("specwindow", CP_STRING, window))
strcpy(window, "blackman");
if (eq(window, "none"))
for (i = 0; i < tlen; i++)
win[i] = 1.0;
else if (eq(window, "rectangular"))
for (i = 0; i < tlen; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 1.0;
}
else if (eq(window, "triangle") || eq(window, "bartlet") || eq(window, "bartlett"))
for (i = 0; i < tlen; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 2.0 - fabs(2+4*(time[i]-maxt)/span);
}
else if (eq(window, "hann") || eq(window, "hanning") || eq(window, "cosine"))
for (i = 0; i < tlen; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 1.0 - cos(2*M_PI*(time[i]-maxt)/span);
}
else if (eq(window, "hamming"))
for (i = 0; i < tlen; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 1.0 - 0.46/0.54*cos(2*M_PI*(time[i]-maxt)/span);
}
else if (eq(window, "blackman"))
for (i = 0; i < tlen; i++) {
if (maxt-time[i] > span) {
win[i] = 0;
} else {
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, "flattop"))
for (i = 0; i < tlen; i++) {
if (maxt-time[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", CP_NUM, &order))
order = 2;
if (order < 2)
order = 2;
sigma = 1.0/order;
scale = 0.83/sigma;
for (i = 0; i < tlen; i++) {
if (maxt-time[i] > span)
win[i] = 0;
else
win[i] = scale*exp(-0.5 * pow((time[i]-maxt/2)/(sigma*maxt/2), 2));
}
} else {
fprintf(cp_err, "Warning: unknown window type %s\n", window);
goto done;
}
}
maxt = time[tlen-1];
if (!cp_getvar("specwindow", CP_STRING, window))
strcpy(window, "blackman");
if (!cp_getvar("specwindoworder", CP_NUM, &order))
order = 2;
if (order < 2)
order = 2;
if (fft_windows(window, win, time, tlen, maxt, span, order) == 0)
goto done;
names = ft_getpnames(wl, TRUE);
vlist = NULL;

4
src/include/ngspice/fftext.h
View File

@ -19,6 +19,10 @@ int fftInit(int M);
void fftFree(void);
// release storage for all private cosine and bit reversed tables
int
fft_windows(char *window, double *win, double *time, double length, double maxt, double span, int order);
// computes some popular window functions
void ffts(double *data, int M, int Rows);
/* Compute in-place complex fft on the rows of the input array */
/* INPUTS */

82
src/maths/fft/fftext.c
View File

@ -9,11 +9,16 @@
#define FFT(a,n) if(!fftInit(roundtol(LOG2(n)))) ffts(a,roundtol(LOG2(n)),1);else printf("fft error\n");
*******************************************************************/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include "fftlib.h"
#include "matlib.h"
#include "ngspice/fftext.h"
#include "ngspice/memory.h"
#define eq(a,b) (!strcmp((a), (b)))
// pointers to storage of Utbl's and BRLow's
static double *UtblArray[8*sizeof(int)] =
{0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0};
@ -83,6 +88,83 @@ void fftFree(void)
}
}
int
fft_windows(char *window, double *win, double *time, double length, double maxt, double span, int order)
{
int i;
double sigma, scale;
/* window functions - should have an average of one */
if (eq(window, "none"))
for (i = 0; i < length; i++)
win[i] = 1.0;
else if (eq(window, "rectangular"))
for (i = 0; i < length; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 1.0;
}
else if (eq(window, "triangle") || eq(window, "bartlet") || eq(window, "bartlett"))
for (i = 0; i < length; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 2.0 - fabs(2+4*(time[i]-maxt)/span);
}
else if (eq(window, "hann") || eq(window, "hanning") || eq(window, "cosine"))
for (i = 0; i < length; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 1.0 - cos(2*M_PI*(time[i]-maxt)/span);
}
else if (eq(window, "hamming"))
for (i = 0; i < length; i++) {
if (maxt-time[i] > span)
win[i] = 0.0;
else
win[i] = 1.0 - 0.46/0.54*cos(2*M_PI*(time[i]-maxt)/span);
}
else if (eq(window, "blackman"))
for (i = 0; i < length; i++) {
if (maxt-time[i] > span) {
win[i] = 0;
} else {
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, "flattop"))
for (i = 0; i < length; i++) {
if (maxt-time[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")) {
sigma = 1.0/order;
scale = 0.83/sigma;
for (i = 0; i < length; i++) {
if (maxt-time[i] > span)
win[i] = 0;
else
win[i] = scale*exp(-0.5 * pow((time[i]-maxt/2)/(sigma*maxt/2), 2));
}
} else {
printf( "Warning: unknown window type %s\n", window);
return 0;
}
return 1;
}
/*************************************************
The following calls are easier than calling to fftlib directly.
Just make sure fftlib has been called for each M first.