/* $Id: fourier.c,v 1.1.1.1 2003/07/21 16:18:41 pturner Exp $ * * DFT by definition and FFT */ #include #include #ifndef M_PI # define M_PI 3.14159265358979323846 #endif static int bit_swap(int i, int nu); static void sswap(double *x1, double *x2); /* DFT by definition */ void dft(double *jr, double *ji, int n, int iflag) { int i, j, sgn; double sumr, sumi, tpi, *w, *xr, *xi, co, si, o, on = 1.0 / n; sgn = iflag ? -1 : 1; tpi = 2.0 * M_PI; w = (double *) calloc(n, sizeof(double)); xr = (double *) calloc(n, sizeof(double)); xi = (double *) calloc(n, sizeof(double)); if (w == NULL || xr == NULL || xi == NULL) { errwin("Can't allocate temporary in DFT"); cxfree(w); cxfree(xr); cxfree(xi); return; } for (i = 0; i < n; i++) { w[i] = tpi * i * on; xr[i] = jr[i]; xi[i] = ji[i]; } for (j = 0; j < n; j++) { sumr = 0.0; sumi = 0.0; for (i = 0; i < n; i++) { o = w[j] * i; co = cos(o); si = sin(o); sumr = sumr + xr[i] * co + sgn * xi[i] * si; sumi = sumi + xi[i] * co - sgn * xr[i] * si; } jr[j] = sumr; ji[j] = sumi; } if (sgn == 1) { on = 2.0 * on; } else { on = 0.5; } for (i = 0; i < n; i++) { jr[i] = jr[i] * on; ji[i] = ji[i] * on; } cxfree(w); cxfree(xr); cxfree(xi); } /* this came off the net - if you wrote it let me know and I'll give you credit - PJT */ /* real_data ... ptr. to real part of data to be transformed imag_data ... ptr. to imag " " " " " " inv ..... Switch to flag normal or inverse transform n_pts ... Number of real data points nu ...... logarithm in base 2 of n_pts e.g. nu = 5 if n_pts = 32. */ void fft(double *real_data, double *imag_data, int n_pts, int nu, int inv) { int n2, j, l, i, ib, k, k1, k2; int sgn; double tr, ti, arg, nu1; /* intermediate values in calcs. */ double c, s; /* cosine & sine components of Fourier trans. */ double fac; n2 = n_pts / 2; nu1 = nu - 1.0; k = 0; /* * sign change for inverse transform */ sgn = inv ? -1 : 1; /* * Calculate the componets of the Fourier series of the function */ for (l = 0; l != nu; l++) { do { for (i = 0; i != n2; i++) { j = k / (int) (pow(2.0, nu1)); ib = bit_swap(j, nu); arg = 2.0 * M_PI * ib / n_pts; c = cos(arg); s = sgn * sin(arg); k1 = k; k2 = k1 + n2; tr = *(real_data + k2) * c + *(imag_data + k2) * s; ti = *(imag_data + k2) * c - *(real_data + k2) * s; *(real_data + k2) = *(real_data + k1) - tr; *(imag_data + k2) = *(imag_data + k1) - ti; *(real_data + k1) = *(real_data + k1) + tr; *(imag_data + k1) = *(imag_data + k1) + ti; k++; } k += n2; } while (k < n_pts - 1); k = 0; nu1 -= 1.0; n2 /= 2; } for (k = 0; k != n_pts; k++) { ib = bit_swap(k, nu); if (ib > k) { sswap((real_data + k), (real_data + ib)); sswap((imag_data + k), (imag_data + ib)); } } /* * If calculating the inverse transform, must divide the data by the number of * data points. */ if (inv) fac = 2.0 / n_pts; else fac = 0.5; for (k = 0; k != n_pts; k++) { *(real_data + k) *= fac; *(imag_data + k) *= fac; } } /* * Bit swaping routine in which the bit pattern of the integer i is reordered. * See Brigham's book for details */ static int bit_swap(int i, int nu) { int ib, i1, i2; ib = 0; for (i1 = 0; i1 != nu; i1++) { i2 = i / 2; ib = ib * 2 + i - 2 * i2; i = i2; } return (ib); } /* * Simple exchange routine where *x1 & *x2 are swapped */ static void sswap(double *x1, double *x2) { double temp_x; temp_x = *x1; *x1 = *x2; *x2 = temp_x; }