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(doub le ar[], double ai[])
{
if (ar.length != ai.length)
{
System.err.prin tln("array dimensions must match");
}
else if (!isPowerOfTwo( ar.length))
{
System.err.prin tln("array dimensions must be multiple of 2");
}
else
{
computeFFT(1, ar.length, ar, ai);
if (ai[0] > EPSI)
{
System.err.prin tln("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[i] * scale;
ai[j] = ai[i] * scale;
ar[i] = tempr;
ai[i] = 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[i] - tr;
ai[j] = ai[i] - ti;
ar[i] += tr;
ai[i] += 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((doub le) 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((doub le) l * Math.log(2)) + 0.5);
le1 = le / 2;
u_r = 1.0;
u_i = 0.0;
w_r = Math.cos(Math.P I / (double) le1);
w_i = -Math.sin(Math.P I / (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((doub le) 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((doub le) l * Math.log(2)) + 0.5);
le1 = le / 2;
u_r = 1.0;
u_i = 0.0;
w_r = Math.cos(Math.P I / (double) le1);
w_i = -Math.sin(Math.P I / (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[i] * scale;

ai[j] = ai[i] * scale;

ar[i] = tempr;

ai[i] = 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[i] - tr;

ai[j] = ai[i] - ti;

ar[i] += tr;

ai[i] += ti;

}

}

mmax = istep;

}

}


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

public static void otherCompute(in t 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[i] * scale;
ai[j] = ai[i] * scale;
ar[i] = tempr;
ai[i] = 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[i] - tr;
ai[j] = ai[i] - ti;
ar[i] += tr;
ai[i] += 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[i] = (float) x_re[i / 2];
data[i + 1] = (float) x_im[i / 2];
}
else
{
data[i] = 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
//multidimensiona l 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.2831853071795 9 / (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[i];
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;
}

[Chris Dams]

Re: Fourier transforms (coefficient calculation)...

Hello,

deja@webfuture. com (Matthew) writes:
[color=blue]
>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.[/color]

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[i]=2-2*Math.cos(((do uble)i)/size*2*Math.PI) ;
ai[i]=0;
}
computeFFT(ar,a i);
for(int i=0;i<size;++i)
{ System.out.prin tln("ar["+i+"]="+ar[i]);
System.out.prin tln("ai["+i+"]="+ai[i]);
}
}

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

[E. Robert Tisdale]

Off Topic in comp.lang.c and comp.lang.c++: Fourier transforms (coefficientcal culation)...

Matthew wrote:
[color=blue]
>
>
> 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.
>[/color]

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/