fastfouriertransform.java

来自「一个很好的LIBSVM的JAVA源码。对于要研究和改进SVM算法的学者。可以参考」· Java 代码 · 共 185 行

/*
 *  YALE - Yet Another Learning Environment
 *  Copyright (C) 2001-2004
 *      Simon Fischer, Ralf Klinkenberg, Ingo Mierswa, 
 *          Katharina Morik, Oliver Ritthoff
 *      Artificial Intelligence Unit
 *      Computer Science Department
 *      University of Dortmund
 *      44221 Dortmund,  Germany
 *  email: [e-mail]
 *  web:   http://yale.cs.uni-dortmund.de/
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License as 
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version. 
 *
 *  This program is distributed in the hope that it will be useful, but
 *  WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 *  USA.
 */
package edu.udo.cs.yale.tools.math;

import edu.udo.cs.yale.example.Example;
import edu.udo.cs.yale.example.ExampleSet;
import edu.udo.cs.yale.example.ExampleReader;
import edu.udo.cs.yale.example.Attribute;
import edu.udo.cs.yale.operator.OperatorException;

import java.util.SortedMap;
import java.util.TreeMap;
import java.util.Iterator;
import java.util.List;
import java.util.Collections;

/** <p>
 *  Performs a FastFourierTransform on an array of complex values. The runtime is {@yale.math O(n log n)}.
 *  </p>
 *  <p>
 *  The used direction simply defines used norm factors, it can be omitted but the resulting values may be 
 *  quantitative wrong. 
 *  </p>
 *
 *  @version $Id: FastFourierTransform.java,v 1.2 2004/08/27 11:57:46 ingomierswa Exp $
 */
public class FastFourierTransform {

    /** Specifies the transformation from time into frequency domain. */
    public static final int TIME2FREQUENCY = 0;
    /** Specifies the transformation from frequency into time domain. */
    public static final int FREQUENCY2TIME = 1;

    /** The window function which should be applied before calculating the FT. */
    private int windowFunctionType = WindowFunction.HANNING;

    public FastFourierTransform(int windowFunctionType) {
	this.windowFunctionType = windowFunctionType;
    }

    /** Builds the fourier transform from the values of the first attribute onto the second. */
    public Complex[] getFourierTransform(ExampleSet exampleSet, Attribute source, Attribute target) 
	throws OperatorException {
	SortedMap pairs = new TreeMap();
	ExampleReader reader = exampleSet.getExampleReader();
	while (reader.hasNext()) {
	    Example example = reader.next();
	    pairs.put(new Double(example.getValue(target)), new Double(example.getValue(source)));
	}
	Complex[] complex = new Complex[pairs.size()];
	Iterator p = pairs.values().iterator();
	int k = 0;
	while (p.hasNext()) {
	    complex[k++] = new Complex(((Double)p.next()).doubleValue(), 0.0d);
	}
	Complex[] result = null;
	try {
	    result = getFourierTransform(complex, TIME2FREQUENCY, windowFunctionType);
	} catch (Exception e) {
	    throw new OperatorException("Exception occured during FFT!", e);
	}
	return result;
    }

    /** Normalizes the frequency to the correct value. */
    public static double convertFrequency(double frequency, int nyquist, int totalLength) {
	return (double)frequency / (2.0d * Math.PI) * ((double)totalLength / (double)nyquist);
    }

    /** Performs a fourier transformation in the specified direction. The window function type defines 
     *  one of the possible window functions. */
    public Complex[] getFourierTransform(Complex[] series, int direction, int functionType) throws Exception {
	int n = getGreatestPowerOf2LessThan(series.length);
	WindowFunction filter = new WindowFunction(functionType, n);
	if (n < 2) {
	    throw new Exception("Cannot transform series: not enough attributes!");
	}
        int nu = (int)(Math.log(n)/Math.log(2));
        int n2 = n/2;
        int nu1 = nu - 1;
        double[] xre = new double[n];
        double[] xim = new double[n];
        double tr, ti, p, arg, c, s;
        for (int i = 0; i < n; i++) {
            xre[i] = filter.getFactor(i) * series[i].getReal();
            xim[i] = filter.getFactor(i) * series[i].getImaginary();
        }
        int k = 0;

        for (int l = 1; l <= nu; l++) {
            while (k < n) {
                for (int i = 1; i <= n2; i++) {
                    p = bitrev (k >> nu1, nu);
                    arg = 2 * (double) Math.PI * p / n;
                    c = (double) Math.cos (arg);
                    s = (double) Math.sin (arg);
                    tr = xre[k+n2]*c + xim[k+n2]*s;
                    ti = xim[k+n2]*c - xre[k+n2]*s;
                    xre[k+n2] = xre[k] - tr;
                    xim[k+n2] = xim[k] - ti;
                    xre[k] += tr;
                    xim[k] += ti;
                    k++;
                }
                k += n2;
            }
            k = 0;
            nu1--;
            n2 = n2/2;
        }
        k = 0;
        int r;
        while (k < n) {
            r = bitrev (k, nu);
            if (r > k) {
                tr = xre[k];
                ti = xim[k];
                xre[k] = xre[r];
                xim[k] = xim[r];
                xre[r] = tr;
                xim[r] = ti;
            }
            k++;
        }

	int nyquist = n / 2;
	Complex[] result = new Complex[nyquist];
	switch (direction) {
	    case TIME2FREQUENCY:
		for (int i = 0; i < nyquist; i++)
		    result[i] = new Complex(- xre[i], xim[i]);
		break;
	    case FREQUENCY2TIME:
		for (int i = 0; i < nyquist; i++)
		    result[i] = new Complex(- xre[i], xim[i]);
		break;
	}
	return result;
    }

    /** Calculates the greatest power of 2 which is smaller than the given number. */
    public static int getGreatestPowerOf2LessThan(int n) {
	int power = (int)(Math.log(n)/Math.log(2));
	return (int)Math.pow(2, power);
    }

    /** Calculates ... */
    private int bitrev(int j, double nu) {
        int j2;
        int j1 = j;
        int k = 0;
        for (int i = 1; i <= nu; i++) {
            j2 = j1/2;
            k  = 2*k + j1 - 2*j2;
            j1 = j2;
        }
        return k;
    }
}