Неожиданный результат дискретного анализа спектра Фурье синусоидальной пробы

Question

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

Мой код на самом деле просто копия вышеупомянутого:

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);

    }
}

Когда я запускаю его на этом примере: http://vigtig.it/sin440.wav

Iполучить такой результат:

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

Я бы ожидал самого высокого отклика при 440 Гц, но это не так.Может кто-то увидеть ошибку или просветить меня о том, как я неверно истолковываю результаты?

РЕДАКТИРОВАТЬ

После просмотра преобразования byte / int я изменил сценарий наиспользуйте вместо этого ByteBuffer.Кажется, работает как задумано сейчас.Вот рабочая копия:

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

Преобразование пар байтов в целые числа со знаком кажется неправильным. Расчет частоты кажется неправильным. Возможно, значение длины тоже плохое.

Невозможно интерпретировать результаты DFT, если ввод плохой. Попробуйте нарисовать входной сигнал DFT (форму волны во временной области) и посмотрите, будет ли это выглядеть уместно вначале.

Answer 2

Вам необходимо применить оконную функцию до БПФ, в противном случае вы увидите эффекты спектральной утечки , которая обычно приводит к размытию амплитудного спектра.