signalprocessing.java

基于java的3d开发库。对坐java3d的朋友有很大的帮助。

//===========================================================================
//= References:                                                             =
//=-------------------------------------------------------------------------=
//= Module history:                                                         =
//= - June 8 2005 - Oscar Chavarro: Original base version                   =
//=-------------------------------------------------------------------------=
//= [.wPRIN2007] Princeton University, "Introduction to programming in      =
//=     Java", course notes available at                                    =
//=     http://www.cs.princeton.edu/introcs/97data/                         =
//=     last accessed, june 8 2007.                                         =
//===========================================================================

package vsdk.toolkit.processing;

import vsdk.toolkit.common.Complex;

public class SignalProcessing extends ProcessingElement {

    /**
    Compute the FFT of x[], assuming its length is a power of 2.
    Bare bones implementation that runs in O(N log N) time. Design goal
    is to optimize the clarity of the code, rather than performance.
    Not the most memory efficient algorithm (because it uses
    an object type for representing complex numbers and because
    it re-allocates memory for the subarray, instead of doing
    in-place or reusing a single temporary array)
    Current implementation of complex number based on (borrowed from, jeje)
    sample in [.wPRIN2007].9 course notes.
    */
    public static Complex[] fft(Complex[] x) {
        int N = x.length;

        // Base case
        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;
    }


    /**
    Compute the inverse FFT of x[], assuming its length is a power of 2.
    Bare bones implementation that runs in O(N log N) time. Design goal
    is to optimize the clarity of the code, rather than performance.
    Not the most memory efficient algorithm (because it uses
    an object type for representing complex numbers and because
    it re-allocates memory for the subarray, instead of doing
    in-place or reusing a single temporary array)
    Current implementation of complex number based on (borrowed from, jeje)
    sample in [.wPRIN2007].9 course notes.
    */
    public static Complex[] ifft(Complex[] x) {
        int N = x.length;
        Complex[] y = new Complex[N];

        // Take conjugate
        for (int i = 0; i < N; i++) {
            y[i] = x[i].conjugate();
        }

        // Compute forward FFT
        y = fft(y);

        // Take conjugate again
        for (int i = 0; i < N; i++) {
            y[i] = y[i].conjugate();
        }

        // Divide by N
        for (int i = 0; i < N; i++) {
            y[i] = y[i].times(1.0 / N);
        }

        return y;
    }

    /**
    Compute the circular convolution of x and y.
    Current implementation of complex number based on (borrowed from, jeje)
    sample in [.wPRIN2007].9 course notes.
    */
    public static Complex[] circularConvolve(Complex[] x, Complex[] y) {

        // Should probably pad x and y with 0s so that they have same length
        // and are powers of 2
        if (x.length != y.length) {
            throw new RuntimeException("Dimensions don't agree");
        }

        int N = x.length;

        // Compute FFT of each sequence
        Complex[] a = fft(x);
        Complex[] b = fft(y);

        // Point-wise multiply
        Complex[] c = new Complex[N];
        for (int i = 0; i < N; i++) {
            c[i] = a[i].times(b[i]);
        }

        // Compute inverse FFT
        return ifft(c);
    }


    /**
    Compute the linear convolution of x and y.
    Current implementation of complex number based on (borrowed from, jeje)
    sample in [.wPRIN2007].9 course notes.
    */
    public static Complex[] linearConvolve(Complex[] x, Complex[] y) {
        Complex ZERO = new Complex(0, 0);

        Complex[] a = new Complex[2*x.length];
        for (int i = 0;        i <   x.length; i++) a[i] = x[i];
        for (int i = x.length; i < 2*x.length; i++) a[i] = ZERO;

        Complex[] b = new Complex[2*y.length];
        for (int i = 0;        i <   y.length; i++) b[i] = y[i];
        for (int i = y.length; i < 2*y.length; i++) b[i] = ZERO;

        return circularConvolve(a, b);
    }

}

//===========================================================================
//= EOF                                                                     =
//===========================================================================