Browse Source

cmath4.c, com_fft.c, rename `size' --> `N' and `mm' --> `M'

pre-master-46
dwarning 13 years ago
committed by rlar
parent
commit
f989219f6f
  1. 88
      src/frontend/com_fft.c
  2. 58
      src/maths/cmaths/cmath4.c

88
src/frontend/com_fft.c

@ -40,13 +40,13 @@ com_fft(wordlist *wl)
double maxt;
#ifdef GREEN
int mm;
int M;
#else
int sign;
#endif
double *reald = NULL, *imagd = NULL;
int size, order;
int N, order;
double scale;
if (!plot_cur || !plot_cur->pl_scale) {
@ -65,20 +65,20 @@ com_fft(wordlist *wl)
#ifdef GREEN
// size of input vector is power of two and larger than spice vector
size = 1;
mm = 0;
while (size < tlen) {
size <<= 1;
mm++;
N = 1;
M = 0;
while (N < tlen) {
N <<= 1;
M++;
}
#else
/* size of input vector is power of two and larger than spice vector */
size = 1;
while (size < tlen)
size *= 2;
N = 1;
while (N < tlen)
N *= 2;
#endif
/* output vector has length of size/2 */
fpts = size/2;
/* output vector has length of N/2 */
fpts = N/2;
win = TMALLOC(double, tlen);
maxt = time[tlen-1];
@ -145,7 +145,7 @@ com_fft(wordlist *wl)
vec_new(f);
for (i = 0; i<fpts; i++)
freq[i] = i*1.0/span*tlen/size;
freq[i] = i*1.0/span*tlen/N;
tdvec = TMALLOC(double *, ngood);
fdvec = TMALLOC(ngcomplex_t *, ngood);
@ -163,26 +163,26 @@ com_fft(wordlist *wl)
vec = vec->v_link2;
}
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);
printf("FFT: Time span: %g s, input length: %d, zero padding: %d\n", span, N, N-tlen);
printf("FFT: Freq. resolution: %g Hz, output length: %d\n", 1.0/span*tlen/N, fpts);
reald = TMALLOC(double, size);
imagd = TMALLOC(double, size);
reald = TMALLOC(double, N);
imagd = TMALLOC(double, N);
for (i = 0; i<ngood; i++) {
for (j = 0; j < tlen; j++) {
reald[j] = tdvec[i][j]*win[j];
imagd[j] = 0.0;
}
for (j = tlen; j < size; j++) {
for (j = tlen; j < N; j++) {
reald[j] = 0.0;
imagd[j] = 0.0;
}
#ifdef GREEN
// Green's FFT
fftInit(mm);
rffts(reald, mm, 1);
fftInit(M);
rffts(reald, M, 1);
fftFree();
scale = size;
scale = N;
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
for (j = 0; j < fpts; j++) {
fdvec[i][j].cx_real = reald[2*j]/scale;
@ -191,7 +191,7 @@ com_fft(wordlist *wl)
fdvec[i][0].cx_imag = 0;
#else
sign = 1;
fftext(reald, imagd, size, tlen, sign);
fftext(reald, imagd, N, tlen, sign);
scale = 0.66;
for (j = 0; j < fpts; j++) {
@ -220,8 +220,8 @@ com_psd(wordlist *wl)
double **tdvec = NULL;
double *freq, *win = NULL, *time, *ave;
double span, noipower;
int mm;
int size, ngood, fpts, i, j, tlen, jj, smooth, hsmooth;
int M;
int N, ngood, fpts, i, j, tlen, jj, smooth, hsmooth;
char *s;
struct dvec *f, *vlist, *lv = NULL, *vec;
struct pnode *pn, *names = NULL;
@ -259,15 +259,15 @@ com_psd(wordlist *wl)
wl = wl->wl_next;
// size of input vector is power of two and larger than spice vector
size = 1;
mm = 0;
while (size < tlen) {
size <<= 1;
mm++;
N = 1;
M = 0;
while (N < tlen) {
N <<= 1;
M++;
}
// output vector has length of size/2
fpts = size>>1;
// output vector has length of N/2
fpts = N>>1;
win = TMALLOC(double, tlen);
maxt = time[tlen-1];
@ -334,7 +334,7 @@ com_psd(wordlist *wl)
vec_new(f);
for (i = 0; i <= fpts; i++)
freq[i] = i*1./span*tlen/size;
freq[i] = i*1./span*tlen/N;
tdvec = TMALLOC(double*, ngood);
fdvec = TMALLOC(ngcomplex_t*, ngood);
@ -352,11 +352,11 @@ com_psd(wordlist *wl)
vec = vec->v_link2;
}
printf("PSD: Time span: %g s, input length: %d, zero padding: %d\n", span, size, size-tlen);
printf("PSD: Freq. resolution: %g Hz, output length: %d\n", 1.0/span*tlen/size, fpts);
printf("PSD: Time span: %g s, input length: %d, zero padding: %d\n", span, N, N-tlen);
printf("PSD: Freq. resolution: %g Hz, output length: %d\n", 1.0/span*tlen/N, fpts);
reald = TMALLOC(double, size);
imagd = TMALLOC(double, size);
reald = TMALLOC(double, N);
imagd = TMALLOC(double, N);
// scale = 0.66;
@ -366,19 +366,19 @@ com_psd(wordlist *wl)
reald[j] = (tdvec[i][j]*win[j]);
imagd[j] = 0.;
}
for (j = tlen; j < size; j++) {
for (j = tlen; j < N; j++) {
reald[j] = 0.;
imagd[j] = 0.;
}
// Green's FFT
fftInit(mm);
rffts(reald, mm, 1);
fftInit(M);
rffts(reald, M, 1);
fftFree();
scaling = size;
scaling = N;
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
intres = (double)size * (double)size;
intres = (double)N * (double)N;
noipower = fdvec[i][0].cx_real = reald[0]*reald[0]/intres;
fdvec[i][fpts].cx_real = reald[1]*reald[1]/intres;
noipower += fdvec[i][fpts-1].cx_real;
@ -459,11 +459,11 @@ fftext(double *x, double *y, long int n, long int nn, int dir)
long i, i1, j, k, i2, l, l1, l2;
double c1, c2, tx, ty, t1, t2, u1, u2, z;
int m = 0, mm = 1;
int m = 0, M = 1;
/* get the exponent to the base of 2 from the number of points */
while (mm < n) {
mm *= 2;
while (M < n) {
M *= 2;
m++;
}

58
src/maths/cmaths/cmath4.c

@ -507,7 +507,7 @@ cx_group_delay(void *data, short int type, int length, int *newlength, short int
void *
cx_fft(void *data, short int type, int length, int *newlength, short int *newtype, struct plot *pl, struct plot *newpl, int grouping)
{
int i, size, mm, fpts, order;
int i, N, M, fpts, order;
double span, scale, maxt;
double *indata, *xscale;
double *time = NULL, *xtime = NULL, *win = NULL;
@ -530,14 +530,14 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
}
/* size of fft input vector is power of two and larger than spice vector */
size = 1;
mm = 0;
while (size < 2*length) {
size <<= 1;
mm++;
N = 1;
M = 0;
while (N < 2*length) {
N <<= 1;
M++;
}
/* output vector has length of size/2 */
fpts = size/2;
/* output vector has length of N/2 */
fpts = N/2;
*newlength = fpts;
*newtype = VF_COMPLEX;
@ -553,7 +553,7 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
span = pl->pl_scale->v_realdata[length-1] - pl->pl_scale->v_realdata[0];
for (i = 0; i<fpts; i++)
xscale[i] = i*1.0/span*(2*length)/size;
xscale[i] = i*1.0/span*(2*length)/N;
for (i = 0; i<length; i++)
time[i] = pl->pl_scale->v_realdata[i];
@ -572,7 +572,7 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
}
for (i = 0; i < length; i++)
time[i] = i*1.0/span*(2*length)/size;
time[i] = i*1.0/span*(2*length)/N;
span = time[length-1] - time[0];
@ -590,7 +590,7 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
}
reald = TMALLOC(double, size);
reald = TMALLOC(double, N);
xtime = TMALLOC(double, 2*length);
/* Now interpolate the data... */
if (!doubledouble(indata, length, reald)) {
@ -624,19 +624,19 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
sv->v_realdata = xscale;
vec_new(sv);
printf("FFT: Time span: %g s, input length: %d, zero padding: %d\n", span, length, size/2-length);
printf("FFT: Frequency resolution: %g Hz, output length: %d\n", 1.0/span*(2*length)/size, fpts);
printf("FFT: Time span: %g s, input length: %d, zero padding: %d\n", span, length, N/2-length);
printf("FFT: Frequency resolution: %g Hz, output length: %d\n", 1.0/span*(2*length)/N, fpts);
for (i = 0; i < 2*length; i++)
reald[i] = reald[i] * win[i];
for (i = 2*length; i < size; i++)
for (i = 2*length; i < N; i++)
reald[i] = 0.0;
fftInit(mm);
rffts(reald, mm, 1);
fftInit(M);
rffts(reald, M, 1);
fftFree();
scale = size;
scale = N;
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
for (i = 0; i < fpts; i++) {
outdata[i].cx_real = reald[2*i]/scale;
@ -657,7 +657,7 @@ void *
cx_ifft(void *data, short int type, int length, int *newlength, short int *newtype, struct plot *pl, struct plot *newpl, int grouping)
{
ngcomplex_t *indata = (ngcomplex_t *) data;
int i, size, mm, tpts, order;
int i, N, M, tpts, order;
double span, scale, maxt;
double *xscale, *win = NULL;
double *outdata;
@ -679,11 +679,11 @@ cx_ifft(void *data, short int type, int length, int *newlength, short int *newty
}
/* size of ifft input vector is power of two and larger or equal than spice vector */
size = 1;
mm = 0;
while (size < length) {
size <<= 1;
mm++;
N = 1;
M = 0;
while (N < length) {
N <<= 1;
M++;
}
if (pl->pl_scale->v_type == SV_TIME) { /* take the time from transient */
@ -710,7 +710,7 @@ cx_ifft(void *data, short int type, int length, int *newlength, short int *newty
span = pl->pl_scale->v_realdata[tpts-1] - pl->pl_scale->v_realdata[0];
for (i = 0; i<tpts; i++)
xscale[i] = i*1.0/span*length/size;
xscale[i] = i*1.0/span*length/N;
} else {
@ -753,22 +753,22 @@ cx_ifft(void *data, short int type, int length, int *newlength, short int *newty
if (fft_windows(window, win, xscale, tpts, maxt, span, order) == 0)
goto done;
printf("IFFT: Frequency span: %g Hz, input length: %d\n", 1/span*tpts/size*(length-1), length);
printf("IFFT: Frequency span: %g Hz, input length: %d\n", 1/span*tpts/N*(length-1), length);
printf("IFFT: Time resolution: %g s, output length: %d\n", span/(tpts-1), tpts);
reald = TMALLOC(double, 2*size);
reald = TMALLOC(double, 2*N);
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
for (i = 0; i < length; i++) {
reald[2*i] = indata[i].cx_real;
reald[2*i+1] = indata[i].cx_imag;
}
for (i = length; i < size; i++) {
for (i = length; i < N; i++) {
reald[2*i] = 0.0;
reald[2*i+1] = 0.0;
}
fftInit(mm);
riffts(reald, mm, 1);
fftInit(M);
riffts(reald, M, 1);
fftFree();
scale = length;

Loading…
Cancel
Save