/*
* Grace - GRaphing, Advanced Computation and Exploration of data
*
* Home page: http://plasma-gate.weizmann.ac.il/Grace/
*
* Copyright (c) 1991-1995 Paul J Turner, Portland, OR
* Copyright (c) 1996-2000 Grace Development Team
*
* Maintained by Evgeny Stambulchik
*
*
* All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
/*
*
* DFT by definition and FFT
*/
/* This module completely rewritten February 1998 by
Marcus H. Mendenhall,
Vanderbilt University Physics Department
Nashville, Tennessee, USA
email: marcus.h.mendenhall@vanderbilt.edu
The conventional DFT and FFT are reimplemented to avoid massive gratuitous
recalculation of trig functions. Also, the signs and magnitudes are
rationalized so that a forward DFT and a forward FFT give the same result.
Also, if the variable HAVE_FFTW is defined (probably from ./configure),
it uses the FFTW libraries to do all transforms, making the DFT and FFT
entirely equivalent. See the FFTW home site
http://theory.lcs.mit.edu/~fftw/
for copyright and usage information and documentation on the FFTW libraries.
I strongly recommend using them... they work very well.
*/
#include <config.h>
#include <cmath.h>
#include <stdio.h>
#include <stdlib.h>
#include "defines.h"
#include "utils.h"
#include "protos.h"
#ifndef HAVE_FFTW
static int bit_swap(int i, int nu);
/*
DFT by definition
*/
void dft(double *jr, double *ji, int n, int iflag)
{
int i, j, sgn;
double sumr, sumi, tpi, w, *xr, *xi, on = 1.0 / n;
double *cov, *siv, co, si;
int iwrap;
sgn = iflag ? -1 : 1;
tpi = 2.0 * M_PI;
xr = xcalloc(n, SIZEOF_DOUBLE);
xi = xcalloc(n, SIZEOF_DOUBLE);
cov = xcalloc(n, SIZEOF_DOUBLE);
siv = xcalloc(n, SIZEOF_DOUBLE);
if (xr == NULL || xi == NULL || cov == NULL || siv == NULL) {
XCFREE(xr);
XCFREE(xi);
XCFREE(cov);
XCFREE(siv);
return;
}
for (i = 0; i < n; i++) {
w = tpi * i * on;
cov[i] = cos(w);
siv[i] = sin(w)*sgn;
xr[i] = jr[i];
xi[i] = ji[i];
}
for (j = 0; j < n; j++) {
sumr = 0.0;
sumi = 0.0;
for (i = 0, iwrap=0; i < n; i++, iwrap += j) {
if(iwrap >= n) iwrap -= n;
co = cov[iwrap];
si = siv[iwrap];
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 = 1.0 * on;
} else {
on = 1.0;
}
for (i = 0; i < n; i++) {
jr[i] = jr[i] * on;
ji[i] = ji[i] * on;
}
XCFREE(xr);
XCFREE(xi);
XCFREE(cov);
XCFREE(siv);
}
/*
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, i, ib ,mm, k;
int sgn, tstep;
double tr, ti, arg; /* intermediate values in calcs. */
double c, s; /* cosine & sine components of Fourier trans. */
static double *sintab = NULL;
static int last_n = 0;
n2 = n_pts / 2;
if(n_pts==0) {
if(sintab) XCFREE(sintab);
sintab=NULL;
last_n=0;
return; /* just deallocate memory if called with zero points */
} else if (n_pts != last_n) { /* allocate new sin table */
arg=2*M_PI/n_pts;
last_n=0;
if(sintab) XCFREE(sintab);
sintab=(double *) xcalloc(n_pts,SIZEOF_DOUBLE);
if(sintab == NULL) return; /* out of memory! */
for(i=0; i<n_pts; i++) sintab[i] = sin(arg*i);
last_n=n_pts;
}
/*
* sign change for inverse transform
*/
sgn = inv ? -1 : 1;
/* do bit reversal of data in advance */
for (k = 0; k != n_pts; k++) {
ib = bit_swap(k, nu);
if (ib > k) {
fswap((real_data + k), (real_data + ib));
fswap((imag_data + k), (imag_data + ib));
}
}
/*
* Calculate the componets of the Fourier series of the function
*/
tstep=n2;
for (mm = 1; mm < n_pts; mm*=2) {
int sinidx=0, cosidx=n_pts/4;
for(k=0; k<mm; k++) {
c = sintab[cosidx];
s = sgn * sintab[sinidx];
sinidx += tstep; cosidx += tstep;
if(sinidx >= n_pts) sinidx -= n_pts;
if(cosidx >= n_pts) cosidx -= n_pts;
for (i = k; i < n_pts; i+=mm*2) {
double re1, re2, im1, im2;
re1=real_data[i]; re2=real_data[i+mm];
im1=imag_data[i]; im2=imag_data[i+mm];
tr = re2 * c + im2 * s;
ti = im2 * c - re2 * s;
real_data[i+mm] = re1 - tr;
imag_data[i+mm] = im1 - ti;
real_data[i] = re1 + tr;
imag_data[i] = im1 + ti;
}
}
tstep /= 2;
}
/*
* If calculating the inverse transform, must divide the data by the number of
* data points.
*/
if (inv) {
double fac = 1.0 / n_pts;
for (k = 0; k != n_pts; k++) {
*(real_data + k) *= fac;
*(imag_data + k) *= fac;
}
}
}
/*
* Bit swapping 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);
}
#else
/* Start of new FFTW-based transforms by Marcus H. Mendenhall */
#include <fftw.h>
#include <string.h>
static char *wisdom_file=0;
static char *initial_wisdom=0;
static int using_wisdom=0;
static void save_wisdom(void){
FILE *wf;
char *final_wisdom;
final_wisdom=fftw_export_wisdom_to_string();
if(!initial_wisdom || strcmp(initial_wisdom, final_wisdom)) {
wf=fopen(wisdom_file,"w");
if(wf) {
fftw_export_wisdom_to_file(wf);
fclose(wf);
}
}
fftw_free(final_wisdom);
if(initial_wisdom) fftw_free(initial_wisdom);
}
void dft(double *jr, double *ji, int n, int iflag)
{
fftw_plan plan;
int i;
double ninv;
FFTW_COMPLEX *cbuf;
static int wisdom_inited=0;
char *ram_cache_wisdom;
int plan_flags;
if(!wisdom_inited) {
wisdom_inited=1;
wisdom_file=getenv("GRACE_FFTW_WISDOM_FILE");
ram_cache_wisdom=getenv("GRACE_FFTW_RAM_WISDOM");
if(ram_cache_wisdom) sscanf(ram_cache_wisdom, "%d", &using_wisdom);
/* turn on wisdom if it is requested even without persistent storage */
if(wisdom_file && wisdom_file[0] ) {
/* if a file was specified in GRACE_FFTW_WISDOM_FILE, try to read it */
FILE *wf;
fftw_status fstat;
wf=fopen(wisdom_file,"r");
if(wf) {
fstat=fftw_import_wisdom_from_file(wf);
fclose(wf);
initial_wisdom=fftw_export_wisdom_to_string();
} else initial_wisdom=0;
atexit(save_wisdom);
using_wisdom=1; /* if a file is specified, always use wisdom */
}
}
plan_flags=using_wisdom? (FFTW_USE_WISDOM | FFTW_MEASURE) : FFTW_ESTIMATE;
plan=fftw_create_plan(n, iflag?FFTW_BACKWARD:FFTW_FORWARD,
plan_flags | FFTW_IN_PLACE);
cbuf=xcalloc(n, sizeof(*cbuf));
if(!cbuf) return;
for(i=0; i<n; i++) {
cbuf[i].re=jr[i]; cbuf[i].im=ji[i];
}
fftw(plan, 1, cbuf, 1, 1, 0, 1, 1);
fftw_destroy_plan(plan);
if(!iflag) {
ninv=1.0/n;
for(i=0; i<n; i++) {
jr[i]=cbuf[i].re*ninv; ji[i]=cbuf[i].im*ninv;
}
} else {
for(i=0; i<n; i++) {
jr[i]=cbuf[i].re; ji[i]=cbuf[i].im;
}
}
XCFREE(cbuf);
}
void fft(double *real_data, double *imag_data, int n_pts, int nu, int inv)
{
/* let FFTW handle DFT's and FFT's identically */
dft(real_data, imag_data, n_pts, inv);
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1