Unexpected result from discrete fourier spectrum analysis of sine-sample

Question

http://jvalentino2.tripod.com/dft/index.html

My code is really just a copy of the above:

package it.vigtig.realtime.fourier;

import java.io.File;
import java.io.IOException;

import javax.sound.sampled.AudioFormat;
import javax.sound.sampled.AudioInputStream;
import javax.sound.sampled.AudioSystem;
import javax.sound.sampled.DataLine;
import javax.sound.sampled.LineUnavailableException;
import javax.sound.sampled.SourceDataLine;

public class Fourier {
    // Create a global buffer size
    private static final int EXTERNAL_BUFFER_SIZE = 128000;

    public static void main(String[] args) {
        /*
         * This code is based on the example found at:
         * http://www.jsresources.org/examples/SimpleAudioPlayer.java.html
         */

        // Get the location of the sound file
        File soundFile = new File("res/sin440.wav");

        // Load the Audio Input Stream from the file
        AudioInputStream audioInputStream = null;
        try {
            audioInputStream = AudioSystem.getAudioInputStream(soundFile);
        } catch (Exception e) {
            e.printStackTrace();
            System.exit(1);
        }

        // Get Audio Format information
        AudioFormat audioFormat = audioInputStream.getFormat();

        // Handle opening the line
        SourceDataLine line = null;
        DataLine.Info info = new DataLine.Info(SourceDataLine.class,
                audioFormat);
        try {
            line = (SourceDataLine) AudioSystem.getLine(info);
            line.open(audioFormat);
        } catch (LineUnavailableException e) {
            e.printStackTrace();
            System.exit(1);
        } catch (Exception e) {
            e.printStackTrace();
            System.exit(1);
        }

        // Start playing the sound
        line.start();

        // Write the sound to an array of bytes
        int nBytesRead = 0;
        byte[] abData = new byte[EXTERNAL_BUFFER_SIZE];
        while (nBytesRead != -1) {
            try {
                nBytesRead = audioInputStream.read(abData, 0, abData.length);

            } catch (IOException e) {
                e.printStackTrace();
            }
            if (nBytesRead >= 0) {
                int nBytesWritten = line.write(abData, 0, nBytesRead);
            }

        }

        // close the line
        line.drain();
        line.close();

        // Calculate the sample rate
        float sample_rate = audioFormat.getSampleRate();
        System.out.println("sample rate = " + sample_rate);

        // Calculate the length in seconds of the sample
        float T = audioInputStream.getFrameLength()
                / audioFormat.getFrameRate();
        System.out
                .println("T = " + T + " (length of sampled sound in seconds)");

        // Calculate the number of equidistant points in time
        int n = (int) (T * sample_rate) / 2;
        System.out.println("n = " + n + " (number of equidistant points)");

        // Calculate the time interval at each equidistant point
        float h = (T / n);
        System.out.println("h = " + h
                + " (length of each time interval in seconds)");

        float fourierFreq = (sample_rate / ((float) n / 2f));
        System.out.println("Fourier frequency is:" + fourierFreq);

        // Determine the original Endian encoding format
        boolean isBigEndian = audioFormat.isBigEndian();

        // this array is the value of the signal at time i*h
        int x[] = new int[n];

        // convert each pair of byte values from the byte array to an Endian
        // value
        for (int i = 0; i < n * 2; i += 2) {
            int b1 = abData[i];
            int b2 = abData[i + 1];
            if (b1 < 0)
                b1 += 0x100;
            if (b2 < 0)
                b2 += 0x100;

            int value;

            // Store the data based on the original Endian encoding format
            if (!isBigEndian)
                value = (b1 << 8) + b2;
            else
                value = b1 + (b2 << 8);
            x[i / 2] = value;
        }

        // do the DFT for each value of x sub j and store as f sub j
        double maxAmp = 0.0;
        double f[] = new double[n / 2];
        for (int j = 1; j < n / 2; j++) {

            double firstSummation = 0;
            double secondSummation = 0;

            for (int k = 0; k < n; k++) {
                double twoPInjk = ((2 * Math.PI) / n) * (j * k);
                firstSummation += x[k] * Math.cos(twoPInjk);
                secondSummation += x[k] * Math.sin(twoPInjk);
            }

            f[j] = Math.abs(Math.sqrt(Math.pow(firstSummation, 2)
                    + Math.pow(secondSummation, 2)));

            double amplitude = 2 * f[j] / n;
            double frequency = j * h / T * sample_rate;

            if (amplitude > maxAmp) {
                maxAmp = amplitude;
                System.out.println("frequency = " + frequency + ", amp = "
                        + amplitude);
            }
        }
        // System.out.println(maxAmp + "," + maxFreq + "," + maxIndex);

    }
}

When I run it on this sample: http://vigtig.it/sin440.wav

I get this result:

sample rate = 8000.0
T = 0.999875 (length of sampled sound in seconds)
n = 3999 (number of equidistant points)
h = 2.5003127E-4 (length of each time interval in seconds)
Fourier frequency is:4.0010004
frequency = 2.000500202178955, amp = 130.77640790523128
frequency = 4.00100040435791, amp = 168.77080135041228
frequency = 6.001501083374023, amp = 291.55653027302816
frequency = 26.006502151489258, amp = 326.4618004521384
frequency = 40.01000213623047, amp = 2265.126299970012
frequency = 200.05003356933594, amp = 3310.905259926063
frequency = 360.09002685546875, amp = 9452.570363111812

I would expect the highest response at 440 hz, but this is not the case. Can someone see a bug or enlighten me as to how I am misinterpreting the results?

EDIT

After looking through the byte/int conversion, I changed the script to use ByteBuffer instead. It seems to work as intended now. Here is the working copy:

package it.vigtig.realtime.fourier;

import java.io.File;
import java.io.IOException;
import java.nio.ByteBuffer;
import java.nio.ShortBuffer;

import javax.sound.sampled.AudioFormat;
import javax.sound.sampled.AudioInputStream;
import javax.sound.sampled.AudioSystem;
import javax.sound.sampled.DataLine;
import javax.sound.sampled.LineUnavailableException;
import javax.sound.sampled.SourceDataLine;

public class Fourier {
    // Create a global buffer size
    private static final int EXTERNAL_BUFFER_SIZE = 16000*16;

    public static void main(String[] args) {
        /*
         * This code is based on the example found at:
         * http://www.jsresources.org/examples/SimpleAudioPlayer.java.html
         */

        // Get the location of the sound file
        File soundFile = new File("res/saw880.wav");

        // Load the Audio Input Stream from the file
        AudioInputStream audioInputStream = null;
        try {
            audioInputStream = AudioSystem.getAudioInputStream(soundFile);
        } catch (Exception e) {
            e.printStackTrace();
            System.exit(1);
        }

        // Get Audio Format information
        AudioFormat audioFormat = audioInputStream.getFormat();

        // Handle opening the line
        SourceDataLine line = null;
        DataLine.Info info = new DataLine.Info(SourceDataLine.class,
                audioFormat);
        try {
            line = (SourceDataLine) AudioSystem.getLine(info);
            line.open(audioFormat);
        } catch (LineUnavailableException e) {
            e.printStackTrace();
            System.exit(1);
        } catch (Exception e) {
            e.printStackTrace();
            System.exit(1);
        }

        // Start playing the sound
        line.start();

        // Write the sound to an array of bytes
        int nBytesRead = 0;
        byte[] abData = new byte[EXTERNAL_BUFFER_SIZE];
        while (nBytesRead != -1) {
            try {
                nBytesRead = audioInputStream.read(abData, 0, abData.length);

            } catch (IOException e) {
                e.printStackTrace();
            }
            if (nBytesRead >= 0) {
                int nBytesWritten = line.write(abData, 0, nBytesRead);
            }

        }

        // close the line
        line.drain();
        line.close();

        // Calculate the sample rate
        float sample_rate = audioFormat.getSampleRate();
        System.out.println("sample rate = " + sample_rate);

        // Calculate the length in seconds of the sample
        float T = audioInputStream.getFrameLength()
                / audioFormat.getFrameRate();
        System.out
                .println("T = " + T + " (length of sampled sound in seconds)");

        // Calculate the number of equidistant points in time
        int n = (int) (T * sample_rate) / 2;
        System.out.println("n = " + n + " (number of equidistant points)");

        // Calculate the time interval at each equidistant point
        float h = (T / n);
        System.out.println("h = " + h
                + " (length of each time interval in seconds)");

        float fourierFreq = (sample_rate / ((float) n / 2f));
        System.out.println("Fourier frequency is:" + fourierFreq);

        // Determine the original Endian encoding format
        boolean isBigEndian = audioFormat.isBigEndian();

        // this array is the value of the signal at time i*h
        int x[] = new int[n];

        ByteBuffer bb = ByteBuffer.allocate(n * 2);

        for (int i = 0; i < n * 2; i++)
            bb.put(abData[i]);


        // do the DFT for each value of x sub j and store as f sub j
        double maxAmp = 0.0;
        double f[] = new double[n / 2];
        for (int j = 1; j < n / 2; j++) {

            double firstSummation = 0;
            double secondSummation = 0;

            for (int k = 0; k < n; k++) {
                double twoPInjk = ((2 * Math.PI) / n) * (j * k);
                firstSummation += bb.getShort(k) * Math.cos(twoPInjk);
                secondSummation += bb.getShort(k) * Math.sin(twoPInjk);
            }

            f[j] = Math.abs(Math.sqrt(Math.pow(firstSummation, 2)
                    + Math.pow(secondSummation, 2)));

            double amplitude = 2 * f[j] / n;
            double frequency = j * h / T * sample_rate;

            if (amplitude > maxAmp) {
                maxAmp = amplitude;
                System.out.println("frequency = " + frequency*2 + ", amp = "
                        + amplitude);
            }
        }
//       System.out.println(maxAmp + "," + maxFreq + "," + maxIndex);

    }
}

Answer 1

The conversion of byte pairs to signed integers seems wrong. The calculation of frequency seems wrong. Perhaps the length value is bad as well.

Can't interpret the DFT results if the input is bad. Try plotting the DFT input (the time domain waveform) and see if that looks appropriate first.

Answer 2

You need to apply a window function prior to the FFT, otherwise you will see the effects of spectral leakage, which typically results in a smearing of the amplitude spectrum.