/*______________________________________________________________________________
 * 
 * Copyright 2003 Paul Cantrell
 * http://innig.net/util/
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * (1) Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *
 * (2) Redistributions in binary form must reproduce the above copyright
 *     notice, this list of conditions and the following disclaimer in
 *     the documentation and/or other materials provided with the
 *     distribution. 
 *
 * (3) The name of the author may not be used to endorse or promote
 *     products derived from this software without specific prior
 *     written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
 * OF THE POSSIBILITY OF SUCH DAMAGE.
 *_______________________________________________________________________________
 */

package net.innig.math;

/**
    Simple one-dimensional FFT, with some very basic optimizations.  Not particularly high-tech.
*/
public class Fourier1D
    {
    /** Factory method allows future etpansion to optimized subclasses (e.g. for non-powers of two). */
    public static Fourier1D create(int size)
        { return new Fourier1D(size); }
    
    private Fourier1D(int size)
        {
        if(size-1 != ((size ^ (size-1)) >> 1))  // There's probably a better way.
            throw new IllegalArgumentException(size + " is not a power of two");
        
        this.size = size;
        evena  = new double[size];
        evenb  = new double[size];
        odda   = new double[size];
        oddb   = new double[size];
        cosTable = new double[size];
        sinTable = new double[size];
        for(int n = 0; n < size; n++)
            {
            double theta = 2 * Math.PI * n / size;
            cosTable[n] = Math.cos(theta);
            sinTable[n] = Math.sin(theta);
            }
        }

    public void transform(
            double[] fromReal, double[] fromImag,
            double[] toReal,   double[] toImag)
        {
        checkSize(fromReal);
        checkSize(fromImag);
        checkSize(toReal);
        checkSize(toImag);
        
        fft(fromReal, fromImag, 0,
            toReal, toImag, 0,
            evena, evenb, 0,
            odda, oddb, 0,
            1, size);
        }
    
    private void checkSize(double[] a)
        {
        if(a.length != size)
            throw new IllegalArgumentException("wrong size: expected " + size + "; got " + a.length);
        }

    public void inverseTransform(
            double[] fromReal, double[] fromImag,
            double[] toReal,   double[] toImag)
        {
        double sizeInv = 1.0 / size;
        double sizeInvNeg = -sizeInv;

        for(int n = 0; n < size; n++)
            fromImag[n] = -fromImag[n];

        transform(fromReal, fromImag, toReal, toImag);

        for(int n = 0; n < size; n++)
            {
            fromImag[n] = -fromImag[n];
            toReal[n] *= sizeInv;
            toImag[n] *= sizeInvNeg;
            }
        }

    private void fft(
            double[] sa,    double[] sb,    int st,
            double[] da,    double[] db,    int dt,
            double[] evena, double[] evenb, int event,
            double[] odda,  double[] oddb,  int oddt,
            int step, int size)
        {
        switch(size)
            {
            case 1:
                {
                da[st] = sa[dt];
                db[dt] = sb[dt];
                return;
                }
    
            case 2:
                {
                da[dt  ] = sa[st] + sa[st + step];
                db[dt  ] = sb[st] + sb[st + step];
                da[dt+1] = sa[st] - sa[st + step];
                db[dt+1] = sb[st] - sb[st + step];
                return;
                }

            case 4:
                {
                double
                    s1a = sa[st           ], a1b = sb[st           ],
                    s2a = sa[st + step    ], s2b = sb[st + step    ],
                    s3a = sa[st + step * 2], s3b = sb[st + step * 2],
                    s4a = sa[st + step * 3], s4b = sb[st + step * 3];
                
                da[dt  ] = s1a + s2a + s3a + s4a;  db[dt  ] = a1b + s2b + s3b + s4b;
                da[dt+1] = s1a - s2b - s3a + s4b;  db[dt+1] = a1b + s2a - s3b - s4a;
                da[dt+2] = s1a - s2a + s3a - s4a;  db[dt+2] = a1b - s2b + s3b - s4b;
                da[dt+3] = s1a + s2b - s3a - s4b;  db[dt+3] = a1b - s2a - s3b + s4a;
                return;
                }
    
            default:
                {
                fft(sa,    sb,    st,
                    evena, evenb, event,
                    evena, evenb, event + size / 2,
                    odda,  oddb,  oddt + size / 2,
                    step * 2, size / 2);
                    
                fft(sa,    sb,    st + step,
                    odda,  oddb,  oddt,
                    evena, evenb, event + size / 2,
                    odda,  oddb,  oddt + size / 2,
                    step * 2, size / 2);

                for(int j = 0; j < size; j++)
                    {
                    double taua = cosTable[j * step],
                           taub = sinTable[j * step];
            
                    int k = j & (size / 2 - 1);
            
                    da[dt+j] = evena[event + k] + taua * odda[oddt + k] - taub * oddb[oddt + k];
                    db[dt+j] = evenb[event + k] + taua * oddb[oddt + k] + taub * odda[oddt + k];
                    }
                }
            }
        }

    private double[] evena, evenb, odda, oddb, cosTable, sinTable;
    private int size;
    }