• Jeferson Andrade <geojaca bol.com.br> Jan 04 2008
• Matti Niemenmaa <see_signature for.real.address> Jan 04 2008
• bearophile <bearophileHUGS lycos.com> Jan 04 2008

Jeferson Andrade <geojaca bol.com.br> writes:

Hi,
I have this code in Java :

public static Complex[] fft(Complex[] x) {
        int N = x.length;
        System.out.println("-------------");
        System.out.println(x.length);
        // caso comum
        if (N == 1) return new Complex[] { x[0] };

        // radix 2 Cooley-Tukey FFT
        if (N % 2 != 0) { throw new RuntimeException("N is not a power of 2"); }

        // fft of even terms
        Complex[] even = new Complex[N/2];
        for (int k = 0; k < N/2; k++) {
            even[k] = x[2*k];
        }
        Complex[] q = fft(even);

        // fft of odd terms
        Complex[] odd  = even;  // reuse the array
        for (int k = 0; k < N/2; k++) {
            odd[k] = x[2*k + 1];
        }
        Complex[] r = fft(odd);

        // combine
        Complex[] y = new Complex[N];
        for (int k = 0; k < N/2; k++) {
            double kth = -2 * k * Math.PI / N;
            Complex wk = new Complex(Math.cos(kth), Math.sin(kth));
            y[k]       = q[k].plus(wk.times(r[k]));
            y[k + N/2] = q[k].minus(wk.times(r[k]));
        }
        return y;
    }

I changed Complex[] for creal[] and I made this code in D, but it contains
error:

static creal[] fft(creal[] x)
{
        int N = x.length;
        writefln("-------------");
        writefln(x.length);
        // common case
        if (N == 1) return new creal[] { x[0]; }; //the error is in this line

        // radix 2 Cooley-Tukey FFT
        if (N % 2 != 0) { throw new Exception("N nao e uma potencia de 2"); }

        // fft of even terms
        creal[] even = new creal[N/2];
        for (int k = 0; k < N/2; k++) {
            even[k] = x[2*k];
        }
        creal[] q = fft(even);

        // fft of odd terms
        creal[] odd  = even;  // reuse the array
        for (int k = 0; k < N/2; k++) {
            odd[k] = x[2*k + 1];
        }
        creal[] r = fft(odd);

        // combine
        creal[] y = new creal[N];
        creal wk;
        real wk1;
        ireal wk2;
        for (int k = 0; k < N/2; k++) {
            real kth = -2 * k * PI / N;
            //= new creal(Math.cos(kth), Math.sin(kth));
            wk1 = cos(kth);
            wk2 = cast(ireal) sin(kth);
            wk = wk1 + wk2;
            y[k] = q[k] + (wk * r[k]);
            y[k + N/2] = q[k] - (wk * r[k]);
        }

        return y;
}

Somebody can gives a clue to me?
Thank you.

Matti Niemenmaa <see_signature for.real.address> writes:

Jeferson Andrade wrote:
 Hi,
 I have this code in Java :


This kind of question should be posted to the digitalmars.D.learn newsgroup. 
Followup set.

<snip>
         if (N == 1) return new creal[] { x[0]; }; //the error is in this line



If my memory serves me correctly, that creates an array containing only one 
element, whose value is x[0].

A simple and idiomatic D way of doing this is using an array literal:

if (N == 1) return [x[0]];

-- 
E-mail address: matti.niemenmaa+news, domain is iki (DOT) fi

bearophile <bearophileHUGS lycos.com> writes:

Jeferson Andrade:
             //= new creal(Math.cos(kth), Math.sin(kth));
             wk1 = cos(kth);
             wk2 = cast(ireal) sin(kth);
             wk = wk1 + wk2;


D complex numbers are tricky...

import std.stdio, std.math;
void main() {
    int kth = 3;
    auto wk1 = cos(kth);
    auto wk2 = cast(ireal) sin(kth);
    auto wk = wk1 + wk2;
    writefln(wk);

    auto wk3 = cos(kth) + sin(kth) * 1i;
    writefln(wk3);
}

It prints:
-0.989992+0i
-0.989992+0.14112i

Bye,
bearophile