Fourier transforms (coefficient calculation)...

[Matthew]

Hello,

For a few days now, I have tried a number of methods that are supposed
to provide the a(k) and b(k) coefficients for a Fourier series. These
methods are listed at the end of this post (in Java). However, not
one of these methods seems to provide the correct coefficients for the
following function;

f(x) = 2 - 2 * cos(x)

....I never see a '-2' or '2' in the resulting generated arrays. I get
a(0) = 4 which is correct, but there is no sign of a(1) anywhere (when
the expression is changed to f(x) = 2 - 2 * sin(x), things seem to
work).

Does someone have an algorithm that works for the above? I have tried
all the Google searches (the source of the methods below), but to no
avail. Please post some example code either in Java or C++ as all the
'try searching for ...' suggestions seen in other postings have been
fruitless. In the meanwhile, perhaps the stuff below will be of use to
someone.

Thanks,

Matthew


============ Method 1 ===================

public static void computeFFT(double ar[], double ai[])
{
if (ar.length != ai.length)
{
System.err.println("array dimensions must match");
}
else if (!isPowerOfTwo(ar.length))
{
System.err.println("array dimensions must be multiple of 2");
}
else
{
computeFFT(1, ar.length, ar, ai);
if (ai[0] > EPSI)
{
System.err.println("imaginary part of constant not zero");
}
}
}

public static void computeFFT(int sign, int n, double ar[], double
ai[])
{
double scale = 2.0 / (double) n;
int i, j;
for (i = j = 0; i < n; ++i)
{
if (j >= i)
{
double tempr = ar[j] * scale;
double tempi = ai[j] * scale;
ar[j] = ar * scale;
ai[j] = ai * scale;
ar = tempr;
ai = tempi;
}
int m = n / 2;
while ((m >= 1) && (j >= m))
{
j -= m;
m /= 2;
}
j += m;
}

int mmax, istep;
for (mmax = 1, istep = 2 * mmax;
mmax < n;
mmax = istep, istep = 2 * mmax)
{
double delta = sign * Math.PI / (double) mmax;
for (int m = 0; m < mmax; ++m)
{
double w = m * delta;
double wr = Math.cos(w);
double wi = Math.sin(w);
for (i = m; i < n; i += istep)
{
j = i + mmax;
double tr = wr * ar[j] - wi * ai[j];
double ti = wr * ai[j] + wi * ar[j];
ar[j] = ar - tr;
ai[j] = ai - ti;
ar += tr;
ai += ti;
}
}
mmax = istep;
}
}

============= Methd 2 =================

public static double[][] fft_1d(double[][] array)
{

double u_r, u_i, w_r, w_i, t_r, t_i;
int ln, nv2, k, l, le, le1, j, ip, i, n;

n = array.length;
ln = (int) (Math.log((double) n) / Math.log(2) + 0.5);
nv2 = n / 2;
j = 1;
for (i = 1; i < n; i++)
{

if (i < j)
{
t_r = array[i - 1][0];
t_i = array[i - 1][1];
array[i - 1][0] = array[j - 1][0];
array[i - 1][1] = array[j - 1][1];
array[j - 1][0] = t_r;
array[j - 1][1] = t_i;
}

k = nv2;
while (k < j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}

for (l = 1; l <= ln; l++) /* loops thru stages */
{
le = (int) (Math.exp((double) l * Math.log(2)) + 0.5);
le1 = le / 2;
u_r = 1.0;
u_i = 0.0;
w_r = Math.cos(Math.PI / (double) le1);
w_i = -Math.sin(Math.PI / (double) le1);
for (j = 1;
j <= le1;
j++) /* loops thru 1/2 twiddle values per stage */
{
for (i = j;
i <= n;
i += le) /* loops thru points per 1/2 twiddle */
{
ip = i + le1;
t_r = array[ip - 1][0] * u_r - u_i * array[ip - 1][1];
t_i = array[ip - 1][1] * u_r + u_i * array[ip - 1][0];

array[ip - 1][0] = array[i - 1][0] - t_r;
array[ip - 1][1] = array[i - 1][1] - t_i;

array[i - 1][0] = array[i - 1][0] + t_r;
array[i - 1][1] = array[i - 1][1] + t_i;
}
t_r = u_r * w_r - w_i * u_i;
u_i = w_r * u_i + w_i * u_r;
u_r = t_r;
}
}
return array;
}


============ Method 3 =================

public static void fft(double[][] array)
{
double u_r, u_i, w_r, w_i, t_r, t_i;
int ln, nv2, k, l, le, le1, j, ip, i, n;

n = array.length;
ln = (int) (Math.log((double) n) / Math.log(2) + 0.5);
nv2 = n / 2;
j = 1;
for (i = 1; i < n; i++)
{
if (i < j)
{
t_r = array[i - 1][0];
t_i = array[i - 1][1];
array[i - 1][0] = array[j - 1][0];
array[i - 1][1] = array[j - 1][1];
array[j - 1][0] = t_r;
array[j - 1][1] = t_i;
}
k = nv2;
while (k < j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}

/* loops thru stages */
for (l = 1; l <= ln; l++)
{
le = (int) (Math.exp((double) l * Math.log(2)) + 0.5);
le1 = le / 2;
u_r = 1.0;
u_i = 0.0;
w_r = Math.cos(Math.PI / (double) le1);
w_i = -Math.sin(Math.PI / (double) le1);

/* loops thru 1/2 twiddle values per stage */
for (j = 1; j <= le1; j++)
{
/* loops thru points per 1/2 twiddle */
for (i = j; i <= n; i += le)
{
ip = i + le1;
t_r = array[ip - 1][0] * u_r - u_i * array[ip - 1][1];
t_i = array[ip - 1][1] * u_r + u_i * array[ip - 1][0];
array[ip - 1][0] = array[i - 1][0] - t_r;
array[ip - 1][1] = array[i - 1][1] - t_i;
array[i - 1][0] = array[i - 1][0] + t_r;
array[i - 1][1] = array[i - 1][1] + t_i;
}
t_r = u_r * w_r - w_i * u_i;
u_i = w_r * u_i + w_i * u_r;
u_r = t_r;
}
}
}


=============== Method 4 ================


public static void newCompute(int sign, int n, double ar[], double
ai[])
{
double scale = 2.0 / (double) n;

int i, j;

for (i = j = 0; i < n; ++i)
{
if (j >= i)
{
double tempr = ar[j] * scale;

double tempi = ai[j] * scale;

ar[j] = ar * scale;

ai[j] = ai * scale;

ar = tempr;

ai = tempi;

}

int m = n / 2;

while ((m >= 1) && (j >= m))
{
j -= m;

m /= 2;

}

j += m;

}

int mmax, istep;

for (mmax = 1, istep = 2 * mmax;
mmax < n;
mmax = istep, istep = 2 * mmax)
{
double delta = sign * Math.PI / (double) mmax;

for (int m = 0; m < mmax; ++m)
{
double w = m * delta;

double wr = Math.cos(w);

double wi = Math.sin(w);

for (i = m; i < n; i += istep)
{
j = i + mmax;

double tr = wr * ar[j] - wi * ai[j];

double ti = wr * ai[j] + wi * ar[j];

ar[j] = ar - tr;

ai[j] = ai - ti;

ar += tr;

ai += ti;

}

}

mmax = istep;

}

}


============= Method 5 ==============

public static void otherCompute(int sign, int n, double ar[], double
ai[])
{
double scale = 2.0 / (double) n;
int i, j;
for (i = j = 0; i < n; ++i)
{
if (j >= i)
{
double tempr = ar[j] * scale;
double tempi = ai[j] * scale;
ar[j] = ar * scale;
ai[j] = ai * scale;
ar = tempr;
ai = tempi;
}
int m = n / 2;
while ((m >= 1) && (j >= m))
{
j -= m;
m /= 2;
}
j += m;
}

int mmax, istep;
for (mmax = 1, istep = 2 * mmax;
mmax < n;
mmax = istep, istep = 2 * mmax)
{
double delta = sign * Math.PI / (double) mmax;
for (int m = 0; m < mmax; ++m)
{
double w = m * delta;
double wr = Math.cos(w);
double wi = Math.sin(w);
for (i = m; i < n; i += istep)
{
j = i + mmax;
double tr = wr * ar[j] - wi * ai[j];
double ti = wr * ai[j] + wi * ar[j];
ar[j] = ar - tr;
ai[j] = ai - ti;
ar += tr;
ai += ti;
}
}
mmax = istep;
}
}


=============== Method 6 ==================

public static Ret fourn(
double x_re[],
double x_im[],
int nn[],
int ndim,
int isign)
{

Ret ret = new Ret();

float[] data = new float[2 * nn[0]];

for (int i = 0; i < 2 * nn[0]; i = i + 2)
{
if (i < 2 * x_re.length && i < 2 * x_im.length)
{
data = (float) x_re[i / 2];
data[i + 1] = (float) x_im[i / 2];
}
else
{
data = 0;
data[i + 1] = 0;
}
}

//Replaces data by its ndim-dimensional discreate Fourier transform,
//if isign is input as 1. nn[0..ndim-1] is an integer array
containing
//the lengths of each dimension (number of complex values), which
//MUST all be power of 2. "data" is a real array of length twice
the
//product of these lengths, in which the data are stored as in a
//multidimensional complex array: real and imaginary parts of each
//element are in consecutive locations, and the right most index of
//the array inreases most rapidly as one proceeds along data. For a
//two-dimensionalbe array, this is equivalent to storing the array
by
//rows. If isign is input as -1, data is replaced by its inverse
//transform times the product of the lengths of all dimensions.
int idim;
int i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
int ibit, k1, k2, n, nprev, nrem, ntot;
float tempr, tempi;

//Double precision for the trigonometric recurrences
double wtemp, wr, wpr, wpi, wi, theta;

//Compute total number of complex values
for (ntot = 1, idim = 0; idim < ndim; idim++)
ntot *= nn[idim];
nprev = 1;
for (idim = ndim - 1; idim >= 0; idim--)
{
n = nn[idim];
nrem = ntot / (n * nprev);
ip1 = nprev << 1;
ip2 = ip1 * n;
ip3 = ip2 * nrem;
i2rev = 1;

//This is the bit-reversal section of the routine.
for (i2 = 1; i2 <= ip2; i2 += ip1)
{
if (i2 < i2rev)
{
for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2)
{
for (i3 = i1; i3 <= ip3; i3 += ip2)
{
i3rev = i2rev + i3 - i2;
float temp = data[i3 - 1];
data[i3 - 1] = data[i3rev - 1];
data[i3rev - 1] = temp;
temp = data[i3];
data[i3] = data[i3rev];
data[i3rev] = temp;
}
}
}
ibit = ip2 >>> 1;
while (ibit >= ip1 && i2rev > ibit)
{
i2rev -= ibit;
ibit >>>= 1;
}
i2rev += ibit;
}

//Here begins the Danielson-Lanczos section of the routine.
ifp1 = ip1;
while (ifp1 < ip2)
{
ifp2 = ifp1 << 1;
// Initialized for the trig. recurrence
theta = isign * 6.28318530717959 / (ifp2 / ip1);
wtemp = Math.sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = Math.sin(theta);
wr = 1.0;
wi = 0.0;
for (i3 = 1; i3 <= ifp1; i3 += ip1)
{
for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
{
for (i2 = i1; i2 <= ip3; i2 += ifp2)
{
//Danielson-Lanczos formula:
k1 = i2;
k2 = k1 + ifp1;
tempr =
(float) wr * data[k2
- 1]
- (float) wi * data[k2];
tempi =
(float) wr * data[k2]
+ (float) wi * data[k2
- 1];
data[k2 - 1] = data[k1 - 1] - tempr;
data[k2] = data[k1] - tempi;
data[k1 - 1] += tempr;
data[k1] += tempi;
}
}
//Trigonometric recurrence
wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
ifp1 = ifp2;
}
nprev *= n;
}

ret.x_re = new double[data.length / 2];
ret.x_im = new double[data.length / 2];
for (int i = 0; i < data.length; i = i + 2)
{
ret.x_re[i / 2] = data;
ret.x_im[i / 2] = data[i + 1];
}

return ret;

}

===== Class Ret (used above) ======

public class Ret
{
public double[] x_re;
public double[] x_im;
}

 

[E. Robert Tisdale]

For a few days now, I have tried a number of methods that are supposed
to provide the a(k) and b(k) coefficients for a Fourier series. These
methods are listed at the end of this post (in Java). However, not
one of these methods seems to provide the correct coefficients for the
following function;

f(x) = 2 - 2 * cos(x)

...I never see a '-2' or '2' in the resulting generated arrays. I get
a(0) = 4 which is correct, but there is no sign of a(1) anywhere (when
the expression is changed to f(x) = 2 - 2 * sin(x), things seem to
work).

Does someone have an algorithm that works for the above? I have tried
all the Google searches (the source of the methods below), but to no
avail. Please post some example code either in Java or C++ as all the
'try searching for ...' suggestions seen in other postings have been
fruitless. In the meanwhile, perhaps the stuff below will be of use to
someone.

This is Off Topic
in comp.lang.c comp.lang.c++ and probably comp.lang.java.

Try The Object-Oriented Numerics Page

http://www.oonumerics.org/oon/

GSL -- The GNU Scientific Library

http://sources.redhat.com/gsl/

Available C++ Libraries FAQ

http://www.faqs.org/faqs/C++-faq/libraries/part1/

 

[Chris Dams]

Hello,

Does someone have an algorithm that works for the above? I have tried
all the Google searches (the source of the methods below), but to no
avail. Please post some example code either in Java or C++ as all the
'try searching for ...' suggestions seen in other postings have been
fruitless. In the meanwhile, perhaps the stuff below will be of use to
someone.

This is off-topic in at least two ways here, but I feel nice today, so
I tried (after making it compilable) your method number I with the
following main function.

public static void main(String argv[])
{ final int size=8;
double ar[]=new double[size];
double ai[]=new double[size];
for(int i=0;i<size;++i)
{ ar=2-2*Math.cos(((double)i)/size*2*Math.PI);
ai=0;
}
computeFFT(ar,ai);
for(int i=0;i<size;++i)
{ System.out.println("ar["+i+"]="+ar);
System.out.println("ai["+i+"]="+ai);
}
}

The result was

ar[0]=4
ar[1]=-2
ar[7]=-2

while the others were essentially 0. This (ignoring normalization) is what
is expected from transforming 2-2*cos(x) because this is equal to
2-exp(ix)-exp(-ix).

Bye,
Chris Dams