/**
 *   
 *   @author  HarveyD
 *   Dan Harvey - Professor of Computer Science
 *   Southern Oregon University, 1250 Siskiyou Blvd., Ashland, OR 97520-5028
 *   harveyd@sou.edu
 *   @version 1.00
 *
 *   Copyright 2010, all rights reserved
 *
 * This software is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This software 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
 * Lesser General Public License for more details.
 *
 * To receive a copy of the GNU Lesser General Public write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/** Note: Re-sample the data to 10 KHZ to normalize between different record rates.
 * 
 *    The first try used the algorithm in ResampleAudio. It worked fine for
 * 	  many conversions (like 16000 KHZ to 1000 10 KHZ). However, it fails when 
 * 	  converting from 22050 KHZ to 10000 KHZ. We then utilized the built in converters 
 *    provided by the Java Sound System. This seems to work better. The original code 
 *    follows:
 * 
 * 	  double[] samples = timeDomain.getTimeDomainFromAudio(-1, -1);
 *	  samples = ResampleAudio.apply(samples,  (int)frameRate);
 *	  audio.setFrameRate(frameRate);
 *
 *    The above works better after some debugging. The correct version is in the 
 *    Android mobile app.
 *
 */

package cs415;

import java.awt.Point;
import java.io.Serializable;

import org.acorns.audio.NormalizeFrames;
import org.acorns.audio.Pitch;
import org.acorns.audio.SoundDefaults;
import org.acorns.audio.TimeDomain;
import org.acorns.audio.frequencydomain.Cepstrum;
import org.acorns.audio.frequencydomain.FastFourierTransform;
import org.acorns.audio.frequencydomain.HarmonicProductSpectrum;
import org.acorns.audio.frequencydomain.MFCC;
import org.acorns.audio.frequencydomain.MelFilterBank;
import org.acorns.audio.frequencydomain.RastaPLP;
import org.acorns.audio.timedomain.Butterworth;
import org.acorns.audio.timedomain.Filter;
import org.acorns.audio.timedomain.LinearPrediction;
import org.acorns.audio.timedomain.ResampleAudio;
import org.acorns.audio.timedomain.Yin;
import org.acorns.data.SoundData;

public class FeatureData implements Serializable
{
    /** Java serial file version */
    private static final long serialVersionUID = 1;
    
	private static String[] optionText = 
	{   "LPErr", "LPCSum", "Energy",  "0Cross", "AutoCor", "Entropy", 
		"Higuchi", "Katz", "Box", "YIN",  "Harmonic", "Cepstral", "SpecFlux", "MelFlux", 
	};
	
	private static String[] diffText = 
	{   "MFCC", "LPC", "LErr", "LSum", "Ener", "0C", 
	    "Acor", "Ent", "Hig", "Katz", "Box", "Yin", "HAR", "CEP", 
	    "L0", "L1", "M0", "M1", "M2", "M3",    
	};
	
	// The number of LPC coefficients
	private static final int P = SoundDefaults.getLPCCoefficients();
	// The number of CEPSTRAL coefficients to model the vocal track (normally: f + sample rate/1000)
    private static final int C = SoundDefaults.getCepstrumLength();
    // The number of Diff coefficients
    private static final int DIFF = diffText.length;
    
	// The temporal filtering algorithm,
	private static final int TEMPORAL_FILTER = RastaPLP.NONE;
	
	// Number of iterations for convergence of variance and skew statistics of features
	private static final int NORMALIZE_LOOP = 3; 
	
	 /*     
	 *  The following are parameters based on Rabiner's end point algorithm
	 *      
	 *  QUARTER_SEC_FRAMES, ZCROSS_COUNT, STD_MUTIPLE are constants that Rabiner seemingly
	 *  arbitrarily picks. Literature seems to vary in how to pick their values.
	 */
	
	/** Number of high zero crossing frames to decide non-voiced sound */
	private static final int ZCROSS_COUNT = 4; 
	/** Standard deviation multiples */
	private static final int STD_MULTIPLE = 2;
	/** The minimum energy above nose for a sound to be perceived as voiced */
	private static final int DELTA_ENERGY = 10;
	
	// Define the parameter to control degree of dynamic feature linear regression curve fitting
	private final static int D = 1;   // Curve fitting loop goes from -D to +D
	
	/** Option to indicate if mean normalization should be done */
	public static int CMN = 1;
	/** Option to indicate if variance and skew normalization should be done */
	public static int CVN = 2;
	/** Option to convert LPC parameters to cepstrals */
	public static int CEP = 4;
	
	/** Bit to indicate if speech is present in a frame */
	public static int SPEECH = 1;    
	
	/** Bit to indicate if a frame is voiced */
	public static int VOICED = 2;
	
	/** Bit to indicate if a frame is silence */
	public static int SILENCE = 4;
	
	/** Bit to indicate if a frame is a phoneme boundary */
	public static int PHONEME = 8;
	
	/** Starting offset to array of MFCC CEPSTRAL coefficients */
	public static final int CEPSTRAL_COEFFICIENTS = 0;
	
	/* Easier symbol access to specific CEPSTRAL_COEFFICIENTS */
	public static final int MFCC0 = CEPSTRAL_COEFFICIENTS;
	public static final int MFCC1 = CEPSTRAL_COEFFICIENTS + 1;
	public static final int MFCC2 = CEPSTRAL_COEFFICIENTS + 2;
	public static final int MFCC3 = CEPSTRAL_COEFFICIENTS + 3;
	public static final int MFCC4 = CEPSTRAL_COEFFICIENTS + 4;
	public static final int MFCC5 = CEPSTRAL_COEFFICIENTS + 5;
	public static final int MFCC6 = CEPSTRAL_COEFFICIENTS + 6;
	public static final int MFCC7 = CEPSTRAL_COEFFICIENTS + 7;
	public static final int MFCC8 = CEPSTRAL_COEFFICIENTS + 8;
	public static final int MFCC9 = CEPSTRAL_COEFFICIENTS + 9;
	public static final int MFCC10 = CEPSTRAL_COEFFICIENTS + 10;
	public static final int MFCC11 = CEPSTRAL_COEFFICIENTS + 11;
	public static final int MFCC12 = CEPSTRAL_COEFFICIENTS + 12;
	

	/** Starting offset to array of Linear Prediction coefficients */
	public static final int LPC_COEFFICIENTS = SoundDefaults.getCepstrumLength();

	/* Easier symbol access to specific LPC_COEFFICIENTS */
	public static final int LPC0 = LPC_COEFFICIENTS;
	public static final int LPC1 = LPC_COEFFICIENTS + 1;
	public static final int LPC2 = LPC_COEFFICIENTS + 2;
	public static final int LPC3 = LPC_COEFFICIENTS + 3;
	public static final int LPC4 = LPC_COEFFICIENTS + 4;
	public static final int LPC5 = LPC_COEFFICIENTS + 5;
	public static final int LPC6 = LPC_COEFFICIENTS + 6;
	public static final int LPC7 = LPC_COEFFICIENTS + 7;

	public static final int DIFF_COEFFICIENTS = LPC_COEFFICIENTS + P;
	
	/* Easier symbol access to specific DIFF_COEFFICIENTS */
	public static final int FLUX_CEPSTRAL_COEFFICIENTS = DIFF_COEFFICIENTS;
	public static final int FLUX_LPC_COEFFICIENTS = DIFF_COEFFICIENTS + 1;
	public static final int DIFF_LPC_ERROR = DIFF_COEFFICIENTS + 2;
	public static final int DIFF_LPC_SUM = DIFF_COEFFICIENTS + 3;
	public static final int DIFF_ENERGY = DIFF_COEFFICIENTS + 4;
	public static final int DIFF_ZERO_CROSS = DIFF_COEFFICIENTS + 5;
	public static final int DIFF_AUTOCORRELATION_COEFFICIENT = DIFF_COEFFICIENTS + 6;
	public static final int DIFF_ENTROPY = DIFF_COEFFICIENTS + 7;
	public static final int DIFF_HIGUCHI_FRACTAL_DIMENSION = DIFF_COEFFICIENTS + 8;
	public static final int DIFF_KATZ_FRACTAL_DIMENSION = DIFF_COEFFICIENTS + 9;
	public static final int DIFF_BOX_FRACTAL_DIMENSION = DIFF_COEFFICIENTS + 10;
	public static final int DIFF_YIN_PITCH = DIFF_COEFFICIENTS + 11;
	public static final int DIFF_HARMONIC_PITCH = DIFF_COEFFICIENTS + 12;
	public static final int DIFF_CEPSTRAL_PITCH = DIFF_COEFFICIENTS + 13;
	public static final int DIFF_LPC0 = DIFF_COEFFICIENTS + 14;
	public static final int DIFF_LPC1 = DIFF_COEFFICIENTS + 15;
	public static final int DIFF_MFCC0 = DIFF_COEFFICIENTS + 16;
	public static final int DIFF_MFCC1 = DIFF_COEFFICIENTS + 17;
	public static final int DIFF_MFCC2 = DIFF_COEFFICIENTS + 18;
	public static final int DIFF_MFCC3 = DIFF_COEFFICIENTS + 19;

	/** Starting offset to array of Linear Prediction coefficients */
	public static final int LPC_ERROR = DIFF_COEFFICIENTS + DIFF;
	
	/** Offset to the sum of the LPC coefficients */
	public static final int LPC_SUM = LPC_ERROR + 1;
	
	/** Starting offset to array of Linear Prediction coefficients */
	public static final int ENERGY = LPC_SUM + 1;
	
	/** Starting offset to zero crossing feature */
	public static final int ZERO_CROSS = ENERGY + 1;
	
	/** Starting offset to audio correlation feature with delta = 1 */
	public static final int AUTOCORRELATION_COEFFICIENT = ZERO_CROSS + 1;
	
	/** Starting offset to entropy feature */
	public static final int ENTROPY = AUTOCORRELATION_COEFFICIENT + 1;
	
	/** Starting offset to fractal dimension feature using the HIGUCHI algorithm */
	public static final int HIGUCHI_FRACTAL_DIMENSION = ENTROPY + 1;
	
	/** Starting offset to fractal dimension feature using the KATZ algorithm */
	public static final int KATZ_FRACTAL_DIMENSION = HIGUCHI_FRACTAL_DIMENSION + 1;
	
	/** Starting offset to fractal dimension feature using the BOX counting algorithm */
	public static final int BOX_FRACTAL_DIMENSION = KATZ_FRACTAL_DIMENSION + 1;
	
	/** Starting offset to pitch estimate using the YIN algorithm */
	public static final int YIN_PITCH = BOX_FRACTAL_DIMENSION + 1;
	
	/** Starting offset to pitch estimate using the HARMONIC product spectrum */
	public static final int HARMONIC_PITCH = YIN_PITCH + 1;
	
	/** Starting offset to pitch estimate using CEPSTRALS */
	public static final int CEPSTRAL_PITCH = HARMONIC_PITCH + 1;
	
	/** Starting offset to Spectral flux */
	public static final int SPECTRAL_FLUX = CEPSTRAL_PITCH + 1;
	
	/** Starting offset to Mel filter flux */
	public static final int MEL_FLUX = SPECTRAL_FLUX + 1;
	
	/** Number of audio features */
	public static final int FEATURES_LENGTH = MEL_FLUX + 1;
	
	/** Number of audio features including delta and delta delta values */
	public static final int FEATURE_ARRAY_LENGTH = FEATURES_LENGTH * 3;

	/** Features where statistics are needed */
	public static final int[] SPEECH_AVERAGES_FEATURES =
	{   ENERGY, ENERGY + FEATURES_LENGTH, ENERGY + 2*FEATURES_LENGTH,
		ZERO_CROSS, ZERO_CROSS + FEATURES_LENGTH, ZERO_CROSS + 2*FEATURES_LENGTH,
		LPC_COEFFICIENTS, LPC_COEFFICIENTS + FEATURES_LENGTH, LPC_COEFFICIENTS + 2*FEATURES_LENGTH,
		LPC_SUM, LPC_SUM + FEATURES_LENGTH, LPC_SUM + 2*FEATURES_LENGTH,
		LPC_ERROR, LPC_ERROR + FEATURES_LENGTH, LPC_ERROR + 2*FEATURES_LENGTH,

		AUTOCORRELATION_COEFFICIENT, AUTOCORRELATION_COEFFICIENT + FEATURES_LENGTH, AUTOCORRELATION_COEFFICIENT + 2*FEATURES_LENGTH,
		ENTROPY, ENTROPY + FEATURES_LENGTH, ENTROPY + 2*FEATURES_LENGTH,
		HIGUCHI_FRACTAL_DIMENSION, HIGUCHI_FRACTAL_DIMENSION + FEATURES_LENGTH, HIGUCHI_FRACTAL_DIMENSION + 2*FEATURES_LENGTH,
		KATZ_FRACTAL_DIMENSION, KATZ_FRACTAL_DIMENSION + FEATURES_LENGTH, KATZ_FRACTAL_DIMENSION + 2*FEATURES_LENGTH,
		BOX_FRACTAL_DIMENSION, BOX_FRACTAL_DIMENSION + FEATURES_LENGTH, BOX_FRACTAL_DIMENSION + 2*FEATURES_LENGTH,

		YIN_PITCH, YIN_PITCH + FEATURES_LENGTH, YIN_PITCH + 2*FEATURES_LENGTH,
		HARMONIC_PITCH, HARMONIC_PITCH + FEATURES_LENGTH, HARMONIC_PITCH + 2*FEATURES_LENGTH,
		CEPSTRAL_PITCH, CEPSTRAL_PITCH + FEATURES_LENGTH, CEPSTRAL_PITCH + 2*FEATURES_LENGTH,
		SPECTRAL_FLUX, SPECTRAL_FLUX + FEATURES_LENGTH, SPECTRAL_FLUX + 2*FEATURES_LENGTH,
		MEL_FLUX, MEL_FLUX + FEATURES_LENGTH, MEL_FLUX + 2*FEATURES_LENGTH,
	};
	
	/** 
	 * 	This section defines combinations of features that could be
	 *  useful for detecting phoneme boundaries. To add additional features
	 *  or combination of features, adjust the following table and update
	 *  the diffText table at the top of this listing with an appropriate header
	 */
     int[][] distanceFeatures = 
     {  
    	{
    		// Non energy components of the cepstrals
    		CEPSTRAL_COEFFICIENTS+1,
    		CEPSTRAL_COEFFICIENTS+2,  CEPSTRAL_COEFFICIENTS+3, CEPSTRAL_COEFFICIENTS+4,
    		CEPSTRAL_COEFFICIENTS+5,  CEPSTRAL_COEFFICIENTS+6, CEPSTRAL_COEFFICIENTS+7, 
    		CEPSTRAL_COEFFICIENTS+8,  CEPSTRAL_COEFFICIENTS+9, CEPSTRAL_COEFFICIENTS+10,
    		CEPSTRAL_COEFFICIENTS+11, CEPSTRAL_COEFFICIENTS+12,
    	},
    		
    	{
    		// Non energy components of the linear prediction
           	LPC_COEFFICIENTS + 1, 	
           	LPC_COEFFICIENTS + 2, 
           	LPC_COEFFICIENTS + 3, 
           	LPC_COEFFICIENTS + 4, 
           	LPC_COEFFICIENTS + 5, 
           	LPC_COEFFICIENTS + 6, 
           	LPC_COEFFICIENTS + 7, 
    	
    	}, 
    		
    	{  	LPC_ERROR, },
    	{  	LPC_SUM, },
    	{  	ENERGY, },
    	{   ZERO_CROSS, },
    	{	AUTOCORRELATION_COEFFICIENT,	},
    	{   ENTROPY, },
    	{  	HIGUCHI_FRACTAL_DIMENSION, },
    	{   KATZ_FRACTAL_DIMENSION, },
    	{   BOX_FRACTAL_DIMENSION, },
    	{  	YIN_PITCH, },
    	{   HARMONIC_PITCH, },
    	{   CEPSTRAL_PITCH, },
    	{	LPC_COEFFICIENTS, },
    	{	LPC_COEFFICIENTS + 1, },
    	{	CEPSTRAL_COEFFICIENTS, },
    	{	CEPSTRAL_COEFFICIENTS + 1, },
    	{	CEPSTRAL_COEFFICIENTS + 2, },
    	{	CEPSTRAL_COEFFICIENTS + 3, },
    };
    
	
	/** offsets into the statistics array */
	private static int ENERGY_STATS = 0;      
	private static int ZERO_CROSS_STATS = 3;
	
	/** statistics array row for holding mean */
	private int MEAN = 0;
	/** statistics array row for holding standard deviation */
	private int STD = 1;
	/** statistics array row for holding variance */
	private int VARIANCE = 2;
	/** statistics array row for holding skew */
	private int SKEW = 3;
	/** statistics array row for holding kirtosis */
	private int KIRTOSIS = 4;
	
	/** toString() Descriptions */
	private static String[] toStringText; 

	private double[][] features;  // Array to hold computed features for each frame
	private int[] frameType;      // Array to categorize frame types
	
	/** Spectrum and Mel coefficients from previous frame */
	private double[] prevSpectrum;
	private double[] prevMelSpectrum;
	
	private float frameRate; // Frame rate
	private int wStep;  // Window step size
	private int wSize;  // Window size
	private int FFT_Size; // FFT bins;
	private int harmonic_FFT_Size; // FFT for computing pitch
	
	/** Constructor to create a feature list of the audio signal with mean normalization 
	 * 
	 * @param audio Object containing the audio signal
	 * @param options 
	 * 		NORM for mean normalization, 
	 * 		CVN for variance and skew normalization
	 * 		CEP for converting LPC parameters to CEPSTRALS
	 * @param frameRate the desired frame rate 
	 */
	public FeatureData(SoundData audio, float frameRate, int options)
	{
    	this(audio, frameRate, (options & CMN) != 0, (options & CVN) != 0, (options & CEP) != 0);
	}
	
	/** Constructor to create a feature list of the audio signal with mean normalization 
	 * 
	 * @param audio Object containing the audio signal
	 * @param options 
	 * 		NORM for mean normalization, 
	 * 		CVN for variance and skew normalization
	 * 		CEP for converting LPC parameters to CEPSTRALS 
	 */
	public FeatureData(SoundData audio, int options)
	{
    	this(audio, SoundDefaults.getFrameRate(), (options & CMN) != 0, (options & CVN) != 0, (options & CEP) != 0);
	}
	
	/** Constructor to create a feature list of the audio signal 
	 * 
	 * @param audio Object containing the audio signal
	 * @param norm true to perform mean normalization
	 * @param cvn true to perform CVN and skew normalization, in addition to CMN
	 * @param cep true if to convert LPC parameters to CEPSTRALS
	 * @param frameRate The desired frame rate to use
	 */
    public FeatureData(SoundData audio, float frameRate, boolean norm, boolean cvn, boolean cep)
    {
    	
     	// Initialize text for the toString output
    	toStringText = new String[FEATURE_ARRAY_LENGTH];
    	int index;
    	for (int i=0; i<LPC_COEFFICIENTS; i++)
    	{
    		index = i+1;
    		toStringText[i] = "MFCC" + (index-1);
    		toStringText[i+FEATURES_LENGTH] = "DMFCC" + (index-1);
    		toStringText[i+2*FEATURES_LENGTH] = "DDMFCC" + (index-1);
    	}
    	for (int i = LPC_COEFFICIENTS; i< DIFF_COEFFICIENTS; i++)
    	{
    		index = i - LPC_COEFFICIENTS + 1;
    		toStringText[i] = "LPC" + (index-1);
    		toStringText[i+FEATURES_LENGTH] = "DLPC" + (index-1);
    		toStringText[i+2*FEATURES_LENGTH] = "DDLPC" + (index-1);
    	}
    	
    	String text;
    	for (int i= DIFF_COEFFICIENTS; i<LPC_ERROR; i++)
    	{
    		index = i - DIFF_COEFFICIENTS + 1;
    		text = (index<=2) ? "Flux" : "Diff";
    		toStringText[i] = diffText[index-1] + text;
    		toStringText[i+FEATURES_LENGTH] = "D" + diffText[index-1] + text;
    		toStringText[i+2*FEATURES_LENGTH] = "DD" + diffText[index-1] + text;
    	}
    	
    	for (int i=LPC_ERROR; i<FEATURES_LENGTH; i++)
    	{
    		index = i-LPC_ERROR;
    		toStringText[i] = optionText[index];
    		toStringText[i+FEATURES_LENGTH] = "D" + optionText[index];
    		toStringText[i+2*FEATURES_LENGTH] = "DD" + optionText[index];
    	}
    	
    	// Normalized feature frame rate
    	this.frameRate = frameRate;

    	// Window (frame) step (overlap) - 10 ms
    	wStep = (int)(frameRate * SoundDefaults.getWindowShift()/1000); 
    	
    	// Window (frame) size - 25.625 ms
    	wSize = (int)(frameRate * SoundDefaults.getWindowSize()/1000); 
    	// FFT size for the frame rate
    	FFT_Size =  1<<(32 - Integer.numberOfLeadingZeros(wSize - 1));
    	// FFT size for harmonic product spectrum
    	harmonic_FFT_Size = (FFT_Size<2048) ? 2048 : FFT_Size;   
    	
    	updateFeatures(audio, norm, cvn, cep);
    }
    
    /** Get frames of features extracted from raw data */
    public double[][] getFrames()
    {
    	return features;
    }
    
    /** extract features from signal
     * 
     * @param audio SoundData audio object
     * @param norm true to perform mean normalization
     * @param cvn true to perform CVN and skew normalization, in addition to cmn
	 * @param cep true if to convert LPC parameters to CEPSTRALS
     */
    private void updateFeatures(SoundData audio, boolean norm, boolean cvn, boolean cep)
    {
    	// Cannot extract features from an audio with nothing recorded
    	if (audio.isRecorded())
    	{  // Remove DC
			// Normalize the sample rate to accommodate recordings at differing rates
	   		ResampleAudio.apply(frameRate, audio);

	   		// Get time domain data from the audio object
	        TimeDomain timeDomain = new TimeDomain(audio);
	        double[] samples = timeDomain.getTimeDomainFromAudio();   

	        if (timeDomain!=null)
	        {
	            TimeDomain.removeDC(samples);
	            timeDomain.saveTimeDomainIntoAudio(samples);
	        }
	        extractFeatures(samples, norm, cvn, cep);
	        normalizePitch();
    	}
    	
    }
    
	/** Method to extract features from a time domain array in complex format
	 * 
	 * @param array of samples (even indices - real, odd indices - imaginary)
	 * @param norm true if to perform mean normalization
	 * @param cvn true if to perform variance and skew normalization, in addition to mean normalization
	 * @param cep true if to convert LPC parameters to CEPSTRALS
	 * @return feature list
	 * 
	 */
	private double[][] extractFeatures(double[] samples, boolean norm, boolean cvn, boolean cep)
	{
        // Create band pass filter for pitch extraction
        Butterworth pitchFilter =  new Butterworth
        		(1.0*Pitch.MIN_PITCH_FREQ/frameRate, 1.0*Pitch.MAX_PITCH_FREQ/frameRate, true);
        double[] pitchSamples = pitchFilter.applyFilter(samples, true);

        // Get the band pass filter for other features
        double minFreq = Math.min(SoundDefaults.getMinFreq(), frameRate/2);
        Butterworth btw = new Butterworth(minFreq/frameRate, false);
        samples = btw.applyFilter(samples, true);
        
        // Real (non complex and not filtered) representation of the signal
        double[] complex = new double[samples.length*2];
        for (int s=0; s<samples.length; s++)  {  complex[s*2] = samples[s]; }
        
        // Create MEL Filter Bank
        MelFilterBank melFilter = new MelFilterBank(frameRate, FFT_Size);

        // Compute windowing parameters
        // Window size, step size between windows, FFT size, number of windows 
        float frameRate = SoundDefaults.getFrameRate();
 
        // Create windowing function
        Filter filter = new Filter();
        double[] hammingFilter = filter.createHammingWindow(wSize);
        double[] harmonicFilter = filter.createHammingWindow(harmonic_FFT_Size);
        
        // Create Pitch detection objects
        // Create objects to perform digital signal algorithms 
        Pitch yin = new Yin(frameRate, wSize);
        Pitch harmonic  = new HarmonicProductSpectrum(frameRate, harmonic_FFT_Size);
        Pitch cepstrum = new Cepstrum(frameRate, harmonic_FFT_Size);

        FastFourierTransform fourier = new FastFourierTransform(FFT_Size);
        double[] fftWindow = new double[2*FFT_Size];		// Frequency domain window
        double[] spectrum = new double[FFT_Size];		// For computing spectral flux
        double[] filtered = new double[wSize];          // Band pass filtered time domain window
        double[] pitchWindow = new double[wSize];       // Band pass filtered for pitch
        int startFrame, endFrame;

        FastFourierTransform harmonicFourier = new FastFourierTransform(harmonic_FFT_Size);
        double[] harmonicWindow = new double[2*harmonic_FFT_Size];
 
        // Instantiate the object for computing the RASTA PLP of each frame
        RastaPLP rasta = new RastaPLP(frameRate, FFT_Size, TEMPORAL_FILTER);
        
        // Create array references
        double[] power, rastaCoefficients, melSpectrum, plpWindow;
        
        int frames = samples.length / wStep;
        // Note that we might have to subtract twice, because there could be
        // two partial frames. However this computation agrees with Sphinx 4.
        if (frames*wStep + wSize > samples.length) frames--;
        
        features = new double[frames][FEATURE_ARRAY_LENGTH];
        double prior = 0, nextPrior, energy;
        if (frames==0) return features;

        for (int frame = 0; frame< frames; frame++)
        {
           	startFrame = 2*frame*wStep;
        	endFrame = startFrame + wSize*2;
        	if (endFrame>complex.length) 
        	{
        		endFrame = complex.length; 
        	}
        	
        	/* Clear data left over from the previous frame */
        	if (endFrame - startFrame != wSize*2)
        	{
        		for (int s=0; s<2*wSize; s+=2)
        		{
        			filtered[s/2] = pitchWindow[s/2] = fftWindow[s] = fftWindow[s+1] = 0;	
        		}
        	}
    		for (int s=wSize; s<2*FFT_Size; s++) fftWindow[s] = 0;

        	
        	/* Load windows for processing the next frame */
        	System.arraycopy(complex, startFrame, fftWindow, 0, endFrame - startFrame);
            System.arraycopy(samples,  startFrame/2, filtered, 0, (endFrame - startFrame)/2);
            System.arraycopy(pitchSamples, startFrame/2, pitchWindow, 0, (endFrame - startFrame)/2);
           	for (int i=0; i<wSize; i++)	{ harmonicWindow[i*2] = pitchWindow[i];	}
     
            // Compute the number of zero crossings per window step size
            features[frame][ZERO_CROSS]   	     
       	    		= TimeDomain.zeroCrossings(filtered, 0, wStep, 1);
                    	
            // Get the fractal dimension of the window
            features[frame][HIGUCHI_FRACTAL_DIMENSION] 
            		= filter.getHiguchiFractalDimension(filtered);
            features[frame][KATZ_FRACTAL_DIMENSION] 
            		= filter.getKatzFractalDimension(filtered);
            
    		int B = (int)(Math.log(filtered.length/4)/Math.log(2));
            features[frame][BOX_FRACTAL_DIMENSION] 
            		= filter.getBoxCountingFractalDimension(filtered, B);

            // Get the normalized auto correlation coefficient
            features[frame][AUTOCORRELATION_COEFFICIENT] = filter.autoCorrelate(filtered);
             
           	features[frame][YIN_PITCH] = yin.getPitch(pitchWindow);

            // Calculate the pitch for the current window
            System.arraycopy(fftWindow, 0, harmonicWindow, 0, wSize*2);
            for (int i=wSize*2; i<harmonic_FFT_Size*2; i++) harmonicWindow[i] = 0;
            
        	harmonicWindow = filter.applyWindow(wSize, harmonicFilter, harmonicWindow);
        	harmonicWindow = harmonicFourier.fft(harmonicWindow);
        	
        	features[frame][HARMONIC_PITCH] = harmonic.getPitch(harmonicWindow);
        	features[frame][CEPSTRAL_PITCH] = cepstrum.getPitch(harmonicWindow);
        	
        	// Perform the perceptual linear prediction
        	plpWindow = fftWindow.clone();

        	// Apply the window to the frame
        	plpWindow = filter.applyWindow(wSize, hammingFilter, plpWindow);

        	// Apply the Fast Fourier Transform
        	plpWindow = fourier.fft(plpWindow);

        	// Get spectrum power (square of amplitudes)
        	power = FastFourierTransform.powerSpectrum(plpWindow, null);
        	
        	// Compute spectral flux
        	for  (int i=0; i<power.length; i++)
        	{
        		spectrum[i] = Math.sqrt(power[i]);
        	}
            features[frame][SPECTRAL_FLUX] = computeFlux(spectrum, prevSpectrum);
            prevSpectrum = spectrum.clone();
        	
        	energy = computeShortTermSpectralEnergy(power);
        	
        	// The entire FFT spectrum scales the amplitudes by  window size / 2
        	// The lower half of the FFT spectrum scales amplitudes by window size / 4
        	energy -= 10 * Math.log10((endFrame - startFrame)/4);
        	features[frame][ENERGY] = energy;

        	// Compute the Rasta PLP coefficients
            // Linear prediction coefficients for LPC
        	rastaCoefficients = rasta.computeRastaPLP(power, (C==P) ? P-1 : P);
            if (rastaCoefficients !=null)
            {
            	features[frame][LPC_SUM] = 0;
            	features[frame][LPC_COEFFICIENTS] = rastaCoefficients[0];
            	for (int c = 1; c<P; c++)
            	{
            		features[frame][LPC_SUM] += Math.pow(rastaCoefficients[c] - rastaCoefficients[c-1], 2);
            		features[frame][LPC_COEFFICIENTS + c] = rastaCoefficients[c];
            	}
            	features[frame][LPC_SUM] = Math.sqrt(features[frame][LPC_SUM]);

                // Compute projected linear prediction error of original window 
                features[frame][LPC_ERROR] = rasta.getNormalizedError();
            }
            
            if (cep)
            {
               	rastaCoefficients = LinearPrediction.lpcToCepstral(P-1, P, rastaCoefficients, rasta.getGain());
            	for (int c = 0; c<P; c++)
            	{
            		features[frame][LPC_COEFFICIENTS + c] = rastaCoefficients[c];
            	}
            }
        	
            // Perform the MFCC algorithm
            
            // Preemphasis to emphasize the higher frequency components
            // Emulates how human hearing perceives frequencies (higher and sharper sound)
        	nextPrior = fftWindow[wStep*2-2]; // Last value of the frame step size
        	Filter.preEmphasizer(fftWindow, wSize, prior); // Emphasize higher frequencies
        	prior = nextPrior;
        	
        	// Apply the window to the frame
        	fftWindow = filter.applyWindow(wSize, hammingFilter, fftWindow);

            // Apply the Fast Fourier Transform
        	fftWindow = fourier.fft(fftWindow);

        	// Get spectrum power (square of amplitudes)
        	power = FastFourierTransform.powerSpectrum(fftWindow, null);

        	// Retrieve the MEL Filter Coefficients
            melSpectrum = melFilter.applyMelFilterBank(power);

        	// Compute mel filter flux
            features[frame][MEL_FLUX] = computeFlux(melSpectrum, prevMelSpectrum);
            prevMelSpectrum = melSpectrum.clone();

            features[frame][ENTROPY] = fourier.getEntropy(melSpectrum, 0, melSpectrum.length);
            features[frame] =  MFCC.getMFCCs(melSpectrum, features[frame] );
        }
        
        for (int i=0; i<DIFF; i++)
        {
        	computeDistances(DIFF_COEFFICIENTS+i, distanceFeatures[i]);
        }
        
        // Calculate delta values
        calculateDynamicFeatures(CEPSTRAL_COEFFICIENTS);
        
        // Calculate delta delta values
        calculateDynamicFeatures(CEPSTRAL_COEFFICIENTS + FEATURES_LENGTH);
        
        // Remove channel distortion
        // Normalize cepstrum mean (0), variance (1), skew (0)
    	int start = CEPSTRAL_COEFFICIENTS;
    	int end = FEATURES_LENGTH; 
    	if (norm) normalizeFeatures(start, end, cvn);

    	// Normalize delta mean (0), variance (1), skew (0)
      	start += FEATURES_LENGTH;
      	end += FEATURES_LENGTH;
    	if (norm) normalizeFeatures(start, end, cvn);

    	// Normalize delta-delta mean (0), variance (1), skew (0)
      	start += FEATURES_LENGTH;
      	end += FEATURES_LENGTH;
    	if (norm) normalizeFeatures(start, end, cvn);
       
 	    return features;
	}
	
	/** Categorize frame types as speech, voiced, silence, using a variant of
	 * Rabiner's energy, zero cross algorithm.
	 *
	 *  Rabiner describes the following constants, ITU, ITR, IF in his
	 *   Digital Speech Processing book
	 *    
	 *  In his book, the value for ITU is -20 for ITU, and -30 for IRU, and 
	 *  IF for 35 zero crossings. However these constants seem arbitrary and
	 *  actual implementations seem to vary. In this implementation, we 
	 *  use samples of frames and their means and standard deviations to set the
	 *  values. This should be a more robust approach.
	 *  
	 *  We keep several constants as described below:
	 *  
	 *  DELTA_ENERGY: The minimum decibels above nose that is perceived to be voiced
	 *  STD_MULTIPLE: The number of standard deviation units to use in the calculations
	 *  sampleFrames: The size of a representative sample of frames
	 *  ZERO_CROSS_COUNT: High zero cross frames in sample around voiced sound
	 *  	 *  
	 */
	public Point updateFrameTypes()
	{
		if (features.length==0) return null;
		
		// Frames in a quarter of a second
		final int sampleFrames = (int)frameRate / (4 * wStep); 

		double[][] averages = computeAverages(SPEECH_AVERAGES_FEATURES);
        double[][] noise = computeFrameFeatureData(ENERGY, false);
        
		int[] frames = new int[features.length];
		for (int frame=0; frame<features.length; frame++)
		{
			frames[frame] = frame;
		}
        
        double noiseEnergyMean = noise[MEAN][ENERGY_STATS];
        double noiseEnergySTD = noise[STD][ENERGY_STATS];
        double zCrossMean = averages[MEAN][ZERO_CROSS_STATS];
        double zCrossSTD = averages[STD][ZERO_CROSS_STATS];
        
        double voicedEnergyMean = averages[MEAN][ENERGY_STATS];
        double energy;
       
        double itu = Math.floor(voicedEnergyMean);
        double itr = Math.min(itu-DELTA_ENERGY, Math.floor(noiseEnergyMean) + STD_MULTIPLE * Math.floor(noiseEnergySTD));
        double izct = Math.ceil(zCrossMean) + Math.ceil(zCrossSTD);
        
        frameType = new int[features.length];
        
        boolean foundStartSpeech = false;
        int startSpeech = -1, endSpeech = -1;

        // Compute start of speech (store in startSpeech)
        for (int frame=0; frame<features.length && !foundStartSpeech; frame++)
        {
        	energy = features[frame][ENERGY];
        	if (energy < itr)
        		frameType[frame] &= ~SPEECH;
        	else
        	{
        		int forwardFrame = startSpeech = frame;
        		while (++forwardFrame < features.length)
        		{
        			energy = features[forwardFrame][ENERGY];
        			if (energy < itr) break;
        			else
        			{
        				if (energy < itu) continue;

        				int count = 0;
        				for (int tempFrame = forwardFrame; 
        						tempFrame>Math.max(forwardFrame - sampleFrames, 0); 
        						tempFrame--)
        				{
        					if (features[tempFrame][ZERO_CROSS]>izct) 
        					{
        						if (++count>=ZCROSS_COUNT) startSpeech = tempFrame;
        					}
        				}
        				foundStartSpeech = true;
        				break;
        			}    // end back search for 0 crossings > izct
        		}        // end while search forward for frame > itu
        	}            // end else
        }                // end for search for start of speech
        
        boolean foundEndSpeech = false;
        for (int frame=features.length-1; frame>startSpeech && !foundEndSpeech; frame--)
        {
        	energy = features[frame][ENERGY];
        	if (energy < itr)
        		frameType[frame] &= ~SPEECH;
        	else
        	{
        		int backwardsFrame = endSpeech = frame;
        		while (--backwardsFrame > startSpeech)
        		{
        			energy = features[backwardsFrame][ENERGY];
        			if (energy < itr) break;
        			else
        			{
        				if (energy < itu) continue;
        				
        				int count = 0;
        				int endFrame = Math.min(backwardsFrame+sampleFrames
        									, features.length-1);
        				for (int tempFrame = backwardsFrame;tempFrame<endFrame; tempFrame++)
        				{
        					if (features[tempFrame][ZERO_CROSS]>izct) 
        					{
        						if (++count>=ZCROSS_COUNT) endSpeech = tempFrame;
        					}
        				}
        				foundEndSpeech = true;
        				break;
        			}    // end forward search for 0 crossings > izct
        		}        // end while search backwards for frame > itu
        	}            // end else
        }                // end for search for end of speech
        
        // Mark all the other frames to be speech and to further classify their types
        if (startSpeech<0) startSpeech = 0;
        
    	for (int frame=startSpeech; frame<=endSpeech; frame++)
    	{
    		frameType[frame] |= SPEECH;
    		
    		// designate voiced frames
           	energy = features[frame][ENERGY];
        	if (energy > itu)
        		frameType[frame] |= VOICED;
    		
    	}	// end for    	
    	
    	return new Point(startSpeech, endSpeech);
	}   // end of updateFrameTypes
	
	/** Get the statistical averages for all frames of the audio
	 * 
	 * @param averages The statistical averages in deviation units of the speech features
	 */
	public double[][] getStatisticsSTD(double[][] averages)
	{
		Point bounds = new Point(0, features.length - 1);
		return getStatisticsSTD(bounds, averages);
	}
	
	public double[][] getStatisticsSTD(Point bounds, double[][] features, double[][] averages)
	{
		if (bounds==null)
			bounds = new Point(0, features.length - 1);
		
		int startSpeech = bounds.x;
		int endSpeech = bounds.y;
		
		double[][] results = new double[endSpeech -startSpeech + 1][averages[0].length + 2];
		
		double value, mean, std, units;
		int    feature;
		for (int frame=startSpeech; frame<=endSpeech; frame++)
		{
			results[frame-startSpeech][0] = frame*wStep;
			for (int stat=0; stat<averages[0].length; stat++)
			{
				feature = stat;
				value = features[frame][feature];
	
				mean = averages[MEAN][stat];
				std = averages[STD][stat];
				units = (value - mean) / std;
				results[frame-startSpeech][stat+1] = units;
			}
		}
		return results;
		
	}
	
	/** Get the statistical averages for the selected audio frames
	 * 
	 * @param bounds - bounds.x = starting speech frame, bounds.y = The ending speech frame
	 * @param averages The statistical averages in deviation units of the speech features
	 */
	public double[][] getStatisticsSTD(Point bounds, double[][] averages)
	{
		return getStatisticsSTD(bounds, features, averages);
	}
	
	
    /** Smooth the pitch estimates between frames */
	private void normalizePitch()
	{
	    // Smooth the pitch estimates between frames
        medianOfFive(YIN_PITCH);
        linearSmoother(YIN_PITCH);
        medianOfFive(HARMONIC_PITCH);
        linearSmoother(HARMONIC_PITCH);
        medianOfFive(CEPSTRAL_PITCH);
        linearSmoother(CEPSTRAL_PITCH);
	}
	
	/** Compute the short term energy of a frequency spectrum
	 * 
	 * Don't include low frequencies containing DC offset and hum
	 * 
	 * @param spectrum The FFT spectrum
	 * @return The short term energy in decibels
	 */
	private double computeShortTermSpectralEnergy(double[] spectrum)
	{
		double energy = 0.0;
		
        double minFreq = SoundDefaults.getMinFreq();
        double maxFreq = Math.min(SoundDefaults.getMaxFreq(), frameRate/2);
		double binSize = frameRate / FFT_Size;
		int startBin = (int)(minFreq / binSize);
		int endBin = (int)(maxFreq / binSize);
				
		for (int bin=startBin; bin<endBin; bin++)
		{
			energy += spectrum[bin];
		}
		
		return 10 * Math.log10(energy/(endBin-startBin+1));
	}
	
	/** Estimate mean and standard deviation for a sample of minimum or maximum frames 
	 * 
	 * @param feature The frame feature to consider
	 * @param maximum true if looking for maximum, else minimum
	 */
	private double[][] computeFrameFeatureData(int feature, boolean maximum)
	{
		// Frames in a quarter of a second
		final int sampleFrames = (int)frameRate / (4 * wStep); 

		// Find the frames in the signal with the lowest energy
		int[] selectFrames = new int[sampleFrames];
		double frameValue, featureValue;
		int index;
		
		for (int frame=0; frame<features.length; frame++)
		{
			frameValue = features[frame][feature];
			index = (frame > sampleFrames) ? sampleFrames : frame;
			
			while (index > 0)
			{
				featureValue = features[selectFrames[index-1]][feature];
				if (featureValue > frameValue && maximum) break;
				if (featureValue < frameValue && !maximum) break;
				if (index<sampleFrames) 
					selectFrames[index] = selectFrames[index-1];
				index--;
			}
			if (index < selectFrames.length) selectFrames[index] = frame;
		}
		double[][] totals = computeAverages(selectFrames, SPEECH_AVERAGES_FEATURES);
		return totals;
	}
	
	/** Compute averages for all features and all frames */
	public double[][] computeAverages()
	{
	   	int[] stats = new int[FeatureData.FEATURE_ARRAY_LENGTH];
		for (int feature=0; feature<FeatureData.FEATURE_ARRAY_LENGTH; feature++)
		{
			stats[feature] = feature;
		}
		return computeAverages(stats);		
	}
	
	/** Compute averages for all frames, but selected features
	 * 
	 * @param offsets Offsets to desired features
	 * @return Computed selected averages
	 */
	public double[][] computeAverages(int[] offsets)
	{
		int[] frames = new int[features.length];
		for (int f=0; f<features.length; f++)  frames[f] = f;
		return computeAverages(frames, offsets);
	}
    
	/** Compute statistics on the features object
	 * 
	 * @param frames the frames of interest
	 * @param offsets Array of which features are of interest
	 * @return Computed averages
	 */
	public double[][] computeAverages(int[] frames, int[] offsets)
	{
		double[][] totals = new double[5][offsets.length];
		
		int N = frames.length;
		if (N==0) return totals;
		
		// Total each feature across the frames
		int frame;
		for (int index=0; index < N; index++)
		{
			frame = frames[index];
			if (frame<0 || frame>=features.length) return totals; 
			
			for (int feature=0; feature<offsets.length; feature++)
			{
				totals[MEAN][feature] += features[frame][offsets[feature]];
			}
		}
		
		// Compute the means
		for (int feature = 0; feature<offsets.length; feature++)
		{
			totals[MEAN][feature]/= N;
		}
		if (N<2) return totals;
		
		// Compute totals for variance, stdev, skew, and kirtosis
		double delta;
		for (int index=0; index < N; index++)
		{
			frame = frames[index];

			for (int feature=0; feature<offsets.length; feature++)
			{
				delta = (features[frame][offsets[feature]] - totals[MEAN][feature]);
				totals[VARIANCE][feature] += delta * delta;
				if (N>2) totals[SKEW][feature] += Math.pow(delta, 3);
				if (N>3) totals[KIRTOSIS][feature] += Math.pow(delta,  4);
			}
		}
		
		// Complete the calculations (using Excel formulas)
		double stdev, factor;
		for (int feature=0; feature<offsets.length; feature++)
		{
			totals[VARIANCE][feature] = totals[VARIANCE][feature] /= N - 1;
			totals[STD][feature] = stdev = Math.sqrt(totals[VARIANCE][feature]);
			if (N > 2 && stdev != 0)
			{
				factor = 1.0 * N / ((N-1)*(N-2));
				totals[SKEW][feature] = factor * totals[SKEW][feature] / Math.pow(stdev, 3);
			}
			if (N > 3 && stdev !=0)
			{
				factor = 1.0 * N * (N+1) / ((N-1)*(N-2)*(N-3));
				totals[KIRTOSIS][feature] = factor * totals[KIRTOSIS][feature] / Math.pow(stdev, 4);
				factor = 3.0 * (N-1)*(N-1) / ((N-2)*(N-3));
				totals[KIRTOSIS][feature] -= factor;
			}
		}
		return totals;
	}
	
	/** Create the title for toString output
	 * 
	 * @param offsets array of offsets of the desired features
	 * @param detail true if this is for detailed data rather then summary totals
	 */
	public static String title(int[] offsets, boolean detail)
	{
		StringBuilder build = new StringBuilder();
		String spaces = "        ", text;
	
        if (detail) 	build.append("#### ");
        else build.append("     ");
        
		for (int feature = 0; feature<offsets.length; feature++)
		{
			text = spaces + toStringText[offsets[feature]];
			build.append(String.format("%s ", text.substring(text.length()-10)));
		}
		build.append("\n");
		return build.toString();		
	}
	
	public String title(int[] offsets)
	{
		StringBuilder build = new StringBuilder();
	
        build.append("####,");
		for (int feature = 0; feature<offsets.length; feature++)
		{
			build.append(String.format("%s,", toStringText[offsets[feature]]));
		}
		return build.toString();		
	}
	
	/** Median of five non linear filter to eliminate discontinuities across frames
	 * 
	 * @param feature The feature to smooth
	 */
	private void medianOfFive(int feature)
	{
		if (features==null) return;
		
		double median, save, middle, out[] = new double[features.length];
		
		for (int frame=2; frame<features.length - 2; frame++)
		{
			// median is the lower of top two elements; save is the greater
			median = features[frame+2][feature];
			save = features[frame+1][feature];
			if (median > save)
			{
				median = features[frame+1][feature];
				save = features[frame+2][feature];
			}
			
			// Eliminate the high and low of the bottom two and top to
			if (features[frame-2][feature]<features[frame-1][feature])
			{
				if (features[frame-2][feature]>median) median = features[frame-2][feature];
				if (features[frame-1][feature]<save) save = features[frame-1][feature];
			}
			else
			{
				if (features[frame-1][feature]>median) median = features[frame-1][feature];
				if (features[frame-2][feature]<save) save = features[frame-2][feature];
			}
			
			// Pick the correct median to be the middle of the remaining elements
			middle = features[frame][feature];
			if ((save - middle) * (save-median) <= 0) median = save;
			if ((middle - save) * (middle - median) <= 0) median = middle;
			out[frame] = median;
		}
		
		// Move the medians to the appropriate feature locations
		for (int frame=2; frame<features.length - 2; frame++)
		{
			features[frame][feature] = out[frame];
		}
		
	}
	
	/** A linear smoother to further smooth a feature curve after median of five
	 * 
	 * @param feature The feature in question
	 */
	private void linearSmoother(int feature)
	{
		if (features==null) return;

		double[] out = new double[features.length];
		for (int frame = 2; frame<features.length; frame++)
		{
			out[frame] = features[frame][feature]/4 
			  			+ features[frame-1][feature]/2 + features[frame-2][feature]/4;
		}
		
		// Move the medians to the appropriate feature locations
		for (int frame=2; frame<features.length; frame++)
		{
			features[frame][feature] = out[frame];
		}
	}
	
	/** Calculate the delta value for a set of feature
	 * 
	 * @param frame The frame number containing the feature
	 * @param featureOffset The feature offset into the array of features
	 */
	private void calculateDynamicFeatures(int featureOffset)
	{
		int start, end;
		double numerator, denominator;
		
        // Update the delta and delta delta totals
        for (int frame=0; frame<features.length; frame++)
        {
        	for (int feature=0; feature<FEATURES_LENGTH; feature++)
        	{
        		numerator = denominator = 0;
            	start = (frame<D)? -frame: - D;
            	end = (frame>=features.length - D) ? features.length - frame - 1: +D;
        		for (int d= start; d<= end; d++)
        		{
            		numerator +=  d * features[frame + d][feature + featureOffset];
            		denominator += d * d;
        			
        		}
        		if (denominator !=0)
        			features[frame][feature + featureOffset + FEATURES_LENGTH] 
        		          = numerator / denominator;
        	}
        }
	}
	

	/** Compute distance features 
	 * 
	 * @param distanceFeature The feature in questions (LPC_DISTANCE, CEPSTRAL_DISTANCE, MISC_DISTANCE)
	 * @param Array of feature indices for which to apply distances
	 */
	private void computeDistances(int distanceFeature, int[] distanceFeatures)
	{
		double total, previous, diff, current;
		int feature;
		for (int frame=0; frame<features.length; frame++)
		{
			total =  0;
			for (int index=0; index < distanceFeatures.length; index++)
			{
				feature = distanceFeatures[index];
				previous = (frame==0)  ? 0 : features[frame-1][feature];
				current = features[frame][feature];
				diff = current - previous;
				total += diff * diff;
			}
		
			// Return the average value
			features[frame][distanceFeature] = Math.sqrt(total);
		}
	}
	
	/** absolute difference between frames (1 Norm)
	 * 
	 * @param current Data from current frame
	 * @param previous Data from previous frame (or null)
	 */
	private double computeFlux(double[] current, double[] previous)
	{
		double total = 0.0;
		for (int i=0; i<current.length; i++)
		{
			total += (previous==null) ? Math.abs(current[i]) : Math.abs(current[i]-previous[i]);
		}
		return total;
	}
	

	/** Method to scale the features and normalize the skew
	 *    Because variance and skew cannot be perfectly normalized
	 *      the method iterates a few times to converge on a better result
	 * 
	 * @param start starting feature to normalize
	 * @param end  ending feature to normalize
	 * @param cvn true to perform variance and skew normalization, in addition to mean normalization
	 */
	public void normalizeFeatures(int start, int end, boolean cvn)
	{
		// Iterate to converge mean to 0, variance to 1, and skew to 0
		for (int n=0; n<NORMALIZE_LOOP; n++)
		{
			// Normalize mean of featues to 0
	       features = NormalizeFrames.meanNormalization(features, start, end);
	       if (!cvn) return;
	        
	        // Normalize variances of features unity for better low SNR recognition
	        features = NormalizeFrames.varianceNormalization(features, start, end);
	        // Normalize skew of features to 0 for better Normal Distribution fit
	      	features = NormalizeFrames.skewNormalization(features, start, end);
		}	
	}

	/** Determine if speech in frame is voiced 
	 * 
	 * @param frame Frame number
	 * @return true if yes
	 */
	public boolean isVoiced(int frame)
	{
		if (frameType==null || frame<0 || frame>=frameType.length) return false; 
		return (frameType[frame] & VOICED) != 0;
	}
	
	/** Determine if a particular frame contains speech 
	 * 
	 * @param frame Frame number
	 * @param state true if previouse frame contains speech
	 * @return true if yes
	 */
	public boolean isSpeech(int frame, boolean state)
	{
		if (frameType==null || frame<0 || frame>=frameType.length) return false; 
		return (frameType[frame] & SPEECH) != 0;
	}
	
	/** Get the number of frames in this audio signal */
	public int getSize()
	{
		if (features==null) return 0;
		return features.length;
	}
	
	/** Create string representation of audio signal between starting and ending samples
	 * 
	 * @param bounds The starting and ending offsets (-1,-1) means entire signal
	 * @param offsets The features of interest
	 * @return String representation of the frames in question
	 */
	public String toString(Point bounds, int[] offsets)
	{
		
		if (bounds == null) bounds = new Point(-1,-1);
		
		int start = (bounds.x<0) ? 0 : bounds.x / wStep;
		int end = (bounds.y + wStep - 1)/wStep;
		if (end > features.length)  end = features.length;
		
		if (end<start) return "No frames within the specified bounds";

		int[] frames = new int[end-start];
		for (int frame=start; frame<end; frame++)
		{
			frames[frame-start] = frame;
		}
		return toString(frames, offsets);
		
	}
	
	/** Create a string representation of the selected features
	 * 
	 * @param offsets An array of feature offsets
	 * @param frames array of frames of interest (if null, all frames)
	 * @return The String result
	 * 
	 */
	public String toString(int[] frames, int[] offsets)
	{  
		StringBuilder build = new StringBuilder();
		build.append(title(offsets, true));
	
		int frame;
		double number;
		for (int index=0; index<frames.length; index++)
		{
			frame = frames[index];
			if (frame<0 || frame>=features.length) continue;
			
			build.append(
			String.format("%4d", frame));
			for (int feature=0; feature<offsets.length; feature++)
			{
				number = features[frame][offsets[feature]];
				build.append(String.format("%11.3f", number)); 
			}
			build.append("\n");
		}
		
		// Create the statistical totals
		build.append("\n");
		build.append(title(offsets, false));
		double[][] totals = computeAverages(frames, offsets);
		for (int line=0; line<totals.length; line++)
		{
			build.append("    ");
			for (int feature=0; feature<offsets.length; feature++)
			{
				number =  totals[line][feature];
				build.append(String.format("%11.3f", number)); 
			}
			build.append("\n");
			
		}
		return build.toString();
	}
	
	/** Construct an array of only the designated audio features */
	public double[][] getValidFeatures
	   (int[] validFeatures, double[][] audioFeatures)
	{
		int len = validFeatures.length;
		
		double[][] stats = new double[audioFeatures.length][len];
		for (int frame=0; frame<audioFeatures.length; frame++)
		{	for (int feature=0; feature<len; feature++)
			{
				stats[frame][feature] = audioFeatures[frame][validFeatures[feature]];
			}
		}
		return stats;
	}
	
	@Override
	public String toString()
	{
		int[] offsets = new int[FEATURE_ARRAY_LENGTH];
		for (int feature=0; feature<FEATURE_ARRAY_LENGTH; feature++)
		{
			offsets[feature] = feature;
		}
		
		return toString(new Point(-1,-1), offsets);
	}

}   // End of FeatureData class