discretefouriertransform.java

It is the Speech recognition software. It is platform independent. To execute the source code,

     * and <b>N</b> is the number of points in the FFT.
     * Since <b>w</b> is complex, this is the same as</p>
     * <p><b>Re(weightFft[k]) = cos ( -2 * PI * k / N)</b></p>
     * <p><b>Im(weightFft[k]) = sin ( -2 * PI * k / N)</b></p>
     *
     * @param numberFftPoints number of points in the FFT
     * @param invert whether it's direct (false) or inverse (true) FFT
     *
     */
    private void createWeightFft(int numberFftPoints, boolean invert) {
	/**
	 * weightFFT will have numberFftPoints/2 complex elements.
	 */
        weightFft = new Complex[numberFftPoints >> 1];

	/**
	 * For the inverse FFT,
	 * w = 2 * PI / numberFftPoints;
	 */
	double w = -2 * Math.PI / numberFftPoints;
        if (invert) {
            w = -w;
        }
	
	for (int k = 0; k < (numberFftPoints >> 1); k++) {
            weightFft[k] = new Complex(Math.cos (w * k), Math.sin (w * k));
	}
    }


    /**
     * Reads the next DoubleData object, which is a data frame from which
     * we'll compute the power spectrum. Signal objects just pass through
     * unmodified.
     *
     * @return the next available power spectrum DoubleData object,
     *         or null if no Spectrum object is available
     *
     * @throws DataProcessingException if there is a processing error
     */
    public Data getData() throws DataProcessingException {

	Data input = getPredecessor().getData();
        
        getTimer().start();

        if ((input != null) && (input instanceof DoubleData)) {
            DoubleData data = (DoubleData) input;
            if (!isNumberFftPointsSet) {
                /*
                 * If numberFftPoints is not set by the user,
                 * figure out the numberFftPoints and initialize the
                 * data structures appropriately.
                 */
                if (numberDataPoints != data.getValues().length) {
                    numberDataPoints = data.getValues().length;
                    numberFftPoints = getNumberFftPoints(numberDataPoints);
                    initializeFFT();
                }
            } else {
                /*
                 * Warn if the user-set numberFftPoints is not ideal.
                 */
                if (numberDataPoints != data.getValues().length) {
                    numberDataPoints = data.getValues().length;
                    int idealFftPoints = getNumberFftPoints(numberDataPoints);
                    if (idealFftPoints != numberFftPoints) {
                        logger.warning("User set numberFftPoints (" + 
                                       numberFftPoints + ") is not ideal (" + 
                                       idealFftPoints + ")");
                    }
                }
            }
	    input = process(data);
	}
        
        // At this point - or in the call immediatelly preceding
        // this -, we should have created a cepstrum frame with
        // whatever data came last, even if we had less than
        // window size of data.
        
        getTimer().stop();
	
	return input;
    }


    /**
     * Returns the ideal number of FFT points given the number of samples.
     * The ideal number of FFT points is the closest power of 2 that is 
     * equal to or larger than the number of samples in the incoming window.
     *
     * @param numberSamples the number of samples in the incoming window
     *
     * @return the closest power of 2 that is equal to or larger than
     *         the number of samples in the incoming window
     */
    private final static int getNumberFftPoints(int numberSamples) {
        int fftPoints = 1;
        
        while (fftPoints < numberSamples) {
            fftPoints *= 2;
            if (fftPoints < 1 || fftPoints > Integer.MAX_VALUE) {
                throw new Error("Invalid # of FFT points: " + fftPoints);
            }
        }
        return fftPoints;
    }


    /**
     * Establish the recursion. The FFT computation will be 
     * computed by as a recursion. Each stage in the butterfly
     * will be fully computed during recursion. In fact, we use
     * the mechanism of recursion only because it's the simplest 
     * way of switching the "input" and "output" vectors. The 
     * output of a stage is the input to the next stage. The 
     * butterfly computes elements in place, but we still need to
     * switch the vectors. We could copy it (not very efficient...)
     * or, in C, switch the pointers. We can avoid the pointers by
     * using recursion.
     *
     * @param input input sequence
     * @param output output sequence
     * @param numberFftPoints number of points in the FFT
     * @param invert whether it's direct (false) or inverse (true) FFT
     *
     */
    private void recurseFft(Complex[] input,
                            double[] output,
                            int numberFftPoints,
                            boolean invert) {
        
	double divisor;
        
	/**
	 * The direct and inverse FFT are essentially the same
	 * algorithm, except for two difference: a scaling factor of
	 * "numberFftPoints" and the signal of the exponent in the
	 * weightFft vectors, defined in the method
	 * <code>createWeightFft</code>.
	 */

 	if (!invert) {
	    divisor = 1.0;
	} else {
	    divisor = (double) numberFftPoints;
	}

	/**
	 * Initialize the "from" and "to" variables.
	 */
	for (int i = 0; i < numberFftPoints; i++){
	    to[i].reset();
	    from[i].scaleComplex(input[i], divisor);
	}

	/**
	 * Repeat the recursion log2(numberFftPoints) times,
	 * i.e., we have log2(numberFftPoints) butterfly stages.
	 */
	butterflyStage(from, to, numberFftPoints, numberFftPoints >> 1);

	/**
	 * Compute energy ("float") for each frequency point
	 * from the fft ("complex")
	 */
	if ((this.logBase2NumberFftPoints & 1) == 0) {
	    for (int i = 0; i <= (numberFftPoints >> 1); i++){
		output[i] = from[i].squaredMagnitudeComplex(); 
            }
	} else {
	    for (int i = 0; i <= (numberFftPoints >> 1); i++){
		output[i] = to[i].squaredMagnitudeComplex();
            }
	}
	return;
    }


    /**
     * Compute one stage in the FFT butterfly. The name "butterfly"
     * appears because this method computes elements in pairs, and
     * a flowgraph of the computation (output "0" comes from input "0"
     * and "1" and output "1" comes from input "0" and "1") resembles a
     * butterfly.
     *
     * We repeat <code>butterflyStage</code> for
     * <b>log_2(numberFftPoints)</b> stages, by calling the recursion
     * with the argument <code>currentDistance</code> divided by 2 at
     * each call, and checking if it's still > 0.
     *
     * @param from the input sequence at each stage
     * @param to the output sequence
     * @param numberFftPoints the total number of points
     * @param currentDistance the "distance" between elements in the
     *     butterfly
     */
    private void butterflyStage(Complex[] from, 
				Complex[] to, 
				int numberFftPoints,
				int currentDistance)
    {
	int ndx1From;
	int ndx2From;
	int ndx1To;
	int ndx2To;
	int ndxWeightFft;

	if (currentDistance > 0) {

            int twiceCurrentDistance = 2 * currentDistance;

	    for (int s = 0; s < currentDistance; s++) {
		ndx1From = s;
		ndx2From = s + currentDistance;
		ndx1To = s;
		ndx2To = s + (numberFftPoints >> 1);
		ndxWeightFft = 0;
		while (ndxWeightFft < (numberFftPoints >> 1)) {
		    /**
		     * <b>weightFftTimesFrom2 = weightFft[k]   </b>
		     * <b>                      *from[ndx2From]</b>
		     */
		    weightFftTimesFrom2.multiplyComplex
			(weightFft[ndxWeightFft], from[ndx2From]);
		    /**
		     * <b>to[ndx1To] = from[ndx1From]       </b>
		     * <b>             + weightFftTimesFrom2</b>
		     */
		    to[ndx1To].addComplex
			(from[ndx1From], weightFftTimesFrom2);
		    /**
		     * <b>to[ndx2To] = from[ndx1From]       </b>
		     * <b>             - weightFftTimesFrom2</b>
		     */
		    to[ndx2To].subtractComplex
			(from[ndx1From], weightFftTimesFrom2);
		    ndx1From += twiceCurrentDistance;
		    ndx2From += twiceCurrentDistance;
		    ndx1To += currentDistance;
		    ndx2To += currentDistance;
		    ndxWeightFft += currentDistance;
		}
	    }

	    /**
	     * This call'd better be the last call in this block, so when
	     * it returns we go straight into the return line below.
	     *
	     * We switch the <i>to</i> and <i>from</i> variables,
	     * the total number of points remains the same, and
	     * the <i>currentDistance</i> is divided by 2.
	     */
	    butterflyStage(to, from, numberFftPoints, (currentDistance >> 1));
	}
	return;
    }
}