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

improve amplitude accuracy

This commit is contained in:
dwarning 2010-02-28 17:52:09 +00:00
parent 651e749a03
commit 6d1ab1e505
2 changed files with 91 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

178
src/frontend/com_fft.c
View File

@ -14,7 +14,6 @@ Author: 2008 Holger Vogt
#include "com_fft.h" #include "com_fft.h"
#include "variable.h" #include "variable.h"
#include "missing_math.h"
void void
com_fft(wordlist *wl) com_fft(wordlist *wl)
@ -27,106 +26,111 @@ com_fft(wordlist *wl)
struct dvec *f, *vlist, *lv = NULL, *vec; struct dvec *f, *vlist, *lv = NULL, *vec;
struct pnode *names, *first_name; struct pnode *names, *first_name;
float *reald, *imagd; double *reald, *imagd;
int size, sign, isreal; int size, sign, order;
float scaling; double scale, sigma;
if (!plot_cur || !plot_cur->pl_scale) { if (!plot_cur || !plot_cur->pl_scale) {
fprintf(cp_err, "Error: no vectors loaded.\n"); fprintf(cp_err, "Error: no vectors loaded.\n");
return; return;
} }
if (!isreal(plot_cur->pl_scale) || if (!isreal(plot_cur->pl_scale) ||
((plot_cur->pl_scale)->v_type != SV_TIME)) { ((plot_cur->pl_scale)->v_type != SV_TIME)) {
fprintf(cp_err, "Error: fft needs real time scale\n"); fprintf(cp_err, "Error: fft needs real time scale\n");
return; return;
} }
tlen = (plot_cur->pl_scale)->v_length; tlen = (plot_cur->pl_scale)->v_length;
time = (plot_cur->pl_scale)->v_realdata; time = (plot_cur->pl_scale)->v_realdata;
span = time[tlen-1] - time[0]; span = time[tlen-1] - time[0];
delta_t = span/(tlen - 1); 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; size = 1;
while (size < tlen) while (size < tlen)
size *= 2; size *= 2;
// output vector has length of size/2 /* output vector has length of size/2 */
fpts = size/2; fpts = size/2;
// window function /* window functions - should have an average of one */
win = (double *) tmalloc(tlen * sizeof (double)); win = (double *) tmalloc(tlen * sizeof (double));
{ {
char window[BSIZE_SP]; char window[BSIZE_SP];
double maxt = time[tlen-1]; double maxt = time[tlen-1];
if (!cp_getvar("specwindow", VT_STRING, window)) if (!cp_getvar("specwindow", VT_STRING, window))
strcpy(window,"blackman"); strcpy(window,"blackman");
if (eq(window, "none")) if (eq(window, "none"))
for(i=0; i<tlen; i++) { for(i=0; i<tlen; i++) {
win[i] = 1; win[i] = 1.0;
} }
else if (eq(window, "rectangular")) else if (eq(window, "rectangular"))
for(i=0; i<tlen; i++) {
if (maxt-time[i] > span) {
win[i] = 0;
} else {
win[i] = 1;
}
}
else if (eq(window, "hanning") || eq(window, "cosine"))
for(i=0; i<tlen; i++) { for(i=0; i<tlen; i++) {
if (maxt-time[i] > span) { if (maxt-time[i] > span) {
win[i] = 0; win[i] = 0.0;
} else { } else {
win[i] = 1 - cos(2*M_PI*(time[i]-maxt)/span); 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")) else if (eq(window, "hamming"))
for(i=0; i<tlen; i++) { for(i=0; i<tlen; i++) {
if (maxt-time[i] > span) { if (maxt-time[i] > span) {
win[i] = 0; win[i] = 0.0;
} else { } else {
win[i] = 1 - 0.92/1.08*cos(2*M_PI*(time[i]-maxt)/span); win[i] = 1.0 - 0.46/0.54*cos(2*M_PI*(time[i]-maxt)/span);
} }
} }
else if (eq(window, "triangle") || eq(window, "bartlet")) else if (eq(window, "blackman"))
for(i=0; i<tlen; i++) { for(i=0; i<tlen; i++) {
if (maxt-time[i] > span) { if (maxt-time[i] > span) {
win[i] = 0; win[i] = 0;
} else { } else {
win[i] = 2 - fabs(2+4*(time[i]-maxt)/span); win[i] = 1.0;
}
}
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<tlen; i++) {
if (maxt-time[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.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); win[i] += 0.08/0.42*cos(4*M_PI*(time[i]-maxt)/span);
} }
} }
} else if (eq(window, "gaussian")) { else if (eq(window, "flattop"))
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<tlen; i++) { for(i=0; i<tlen; i++) {
if (maxt-time[i] > span) { if (maxt-time[i] > span) {
win[i] = 0; win[i] = 0;
} else { } else {
win[i] = exp(-0.5*order*(1-2*(maxt-time[i])/span) win[i] = 1.0;
*(1-2*(maxt-time[i])/span))/scale; 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 { else if (eq(window, "gaussian")) {
if (!cp_getvar("specwindoworder", VT_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); fprintf(cp_err, "Warning: unknown window type %s\n", window);
tfree(win); tfree(win);
return; return;
@ -148,7 +152,7 @@ com_fft(wordlist *wl)
continue; continue;
} }
if (!isreal(vec)) { 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->v_name);
vec = vec->v_link2; vec = vec->v_link2;
continue; continue;
@ -166,12 +170,12 @@ com_fft(wordlist *wl)
ngood++; ngood++;
} }
} }
free_pnode_o(first_name); /* h_vogt 081206 */ free_pnode_o(first_name);
if (!ngood) { if (!ngood) {
tfree(win); tfree(win);
return; return;
} }
plot_cur = plot_alloc("spectrum"); plot_cur = plot_alloc("spectrum");
plot_cur->pl_next = plot_list; plot_cur->pl_next = plot_list;
plot_list = plot_cur; plot_list = plot_cur;
@ -189,7 +193,7 @@ com_fft(wordlist *wl)
f->v_realdata = freq; f->v_realdata = freq;
vec_new(f); vec_new(f);
for (i = 0; i<fpts; i++) freq[i] = i*1./span*tlen/size; for (i = 0; i<fpts; i++) freq[i] = i*1.0/span*tlen/size;
tdvec = (double **) tmalloc(ngood * sizeof(double *)); tdvec = (double **) tmalloc(ngood * sizeof(double *));
@ -200,7 +204,7 @@ com_fft(wordlist *wl)
f = alloc(struct dvec); f = alloc(struct dvec);
ZERO(f, struct dvec); ZERO(f, struct dvec);
f->v_name = vec_basename(vec); f->v_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_flags = (VF_COMPLEX | VF_PERMANENT);
f->v_length = fpts; f->v_length = fpts;
f->v_compdata = fdvec[i]; f->v_compdata = fdvec[i];
@ -208,63 +212,61 @@ com_fft(wordlist *wl)
vec = vec->v_link2; vec = vec->v_link2;
} }
sign = 1; sign = 1;
isreal = 1;
printf("FFT: Time span: %g s, input length: %d, zero padding: %d\n", span, size, size-tlen);
reald = (float*)tmalloc(size*sizeof(float)); printf("FFT: Freq. resolution: %g Hz, output length: %d\n", 1.0/span*tlen/size, fpts);
imagd = (float*)tmalloc(size*sizeof(float));
reald = (double*)tmalloc(size*sizeof(double));
printf("CPU: Delta Freq %f Hz, input length %d, output length %d\n", 1./span*tlen/size, size, fpts); imagd = (double*)tmalloc(size*sizeof(double));
for (i = 0; i<ngood; i++) { for (i = 0; i<ngood; i++) {
for (j = 0; j < tlen; j++){ for (j = 0; j < tlen; j++){
reald[j] = tdvec[i][j]*win[j]; reald[j] = tdvec[i][j]*win[j];
imagd[j] = 0; imagd[j] = 0.0;
} }
for (j = tlen; j < size; j++){ for (j = tlen; j < size; j++){
reald[j] = 0; reald[j] = 0.0;
imagd[j] = 0; imagd[j] = 0.0;
} }
fftext(reald, imagd, size, sign); fftext(reald, imagd, size, tlen, sign);
scaling = 0.3; scale = 0.66;
for (j=0;j<fpts;j++){ for (j=0;j<fpts;j++){
fdvec[i][j].cx_real = reald[j]/scaling; fdvec[i][j].cx_real = reald[j]/scale;
fdvec[i][j].cx_imag = imagd[j]/scaling; fdvec[i][j].cx_imag = imagd[j]/scale;
} }
} }
tfree(reald); tfree(reald);
tfree(imagd); tfree(imagd);
tfree(tdvec); tfree(tdvec);
tfree(fdvec); tfree(fdvec);
tfree(win); tfree(win);
} }
static void fftext(float* x, float* y, long int n, int dir) static void fftext(double* x, double* y, long int n, long int nn, int dir)
{ {
/* /*
http://local.wasp.uwa.edu.au/~pbourke/other/dft/ http://local.wasp.uwa.edu.au/~pbourke/other/dft/
download 22.05.08 download 22.05.08
Used with permission from the author Paul Bourke Used with permission from the author Paul Bourke
*/ */
/* /*
This computes an in-place complex-to-complex FFT This computes an in-place complex-to-complex FFT
x and y are the real and imaginary arrays x and y are the real and imaginary arrays
n is the number of points, has to be to the power of 2 n is the number of points, has to be to the power of 2
nn is the number of points w/o zero padded values
dir = 1 gives forward transform dir = 1 gives forward transform
dir = -1 gives reverse transform dir = -1 gives reverse transform
*/ */
long i,i1,j,k,i2,l,l1,l2; long i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z; double c1,c2,tx,ty,t1,t2,u1,u2,z;
int m=0, mm=1; int m=0, mm=1;
/* get the exponent to the base of 2 from the number of points */ /* get the exponent to the base of 2 from the number of points */
while (mm < n) { while (mm < n) {
mm *= 2; mm *= 2;
@ -292,20 +294,20 @@ static void fftext(float* x, float* y, long int n, int dir)
} }
/* Compute the FFT */ /* Compute the FFT */
c1 = -1.0; c1 = -1.0;
c2 = 0.0; c2 = 0.0;
l2 = 1; l2 = 1;
for (l=0;l<m;l++) { for (l=0;l<m;l++) {
l1 = l2; l1 = l2;
l2 <<= 1; l2 <<= 1;
u1 = 1.0; u1 = 1.0;
u2 = 0.0; u2 = 0.0;
for (j=0;j<l1;j++) { for (j=0;j<l1;j++) {
for (i=j;i<n;i+=l2) { for (i=j;i<n;i+=l2) {
i1 = i + l1; i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1]; t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1]; t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1; x[i1] = x[i] - t1;
y[i1] = y[i] - t2; y[i1] = y[i] - t2;
x[i] += t1; x[i] += t1;
y[i] += t2; y[i] += t2;
@ -315,7 +317,7 @@ static void fftext(float* x, float* y, long int n, int dir)
u1 = z; u1 = z;
} }
c2 = sqrt((1.0 - c1) / 2.0); c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1) if (dir == 1)
c2 = -c2; c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0); c1 = sqrt((1.0 + c1) / 2.0);
} }
@ -323,8 +325,8 @@ static void fftext(float* x, float* y, long int n, int dir)
/* Scaling for forward transform */ /* Scaling for forward transform */
if (dir == 1) { if (dir == 1) {
for (i=0;i<n;i++) { for (i=0;i<n;i++) {
x[i] /= n; x[i] /= nn; /* don't consider zero padded values */
y[i] /= n; y[i] /= nn;
} }
} }
} }

2
src/frontend/com_fft.h
View File

@ -10,6 +10,6 @@ extern void free_pnode_o(struct pnode *t);
void com_fft(wordlist *wl); void com_fft(wordlist *wl);
static void fftext(float*, float*, long int, int); static void fftext(double*, double*, long int, long int, int);
#endif #endif