来自正弦样本的离散傅立叶频谱分析的意外结果

Posted 2023-02-25

【中文标题】来自正弦样本的离散傅立叶频谱分析的意外结果

【英文标题】：Unexpected result from discrete fourier spectrum analysis of sine-sample

【发布时间】：2011-12-30 00:36:17

【问题描述】：

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

我得到这个结果：

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 赫兹的响应最高，但事实并非如此。有人可以看到错误或启发我如何误解结果吗？

编辑

查看完字节/整数转换后，我将脚本改为使用 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);

    

【参考方案1】：

字节对到有符号整数的转换似乎是错误的。频率的计算似乎是错误的。也许长度值也很糟糕。

如果输入错误，则无法解释 DFT 结果。尝试绘制 DFT 输入（时域波形）并首先查看是否合适。

【讨论】：

看来这确实是一个糟糕的转换。我改用了 ByteBuffer，现在它似乎可以工作了 :) 谢谢！ （我将工作副本添加到问题中）

【参考方案2】：

您需要在 FFT 之前应用window function，否则您将看到spectral leakage 的效果，这通常会导致幅度谱的拖尾。