FFT true frequency estimation

Trying to compute the true frequency from individual fft data bins. I delved into fft theory but the results of the computation for frequencies which deviate from the bin´s frequency are wrong. On the basis of juce´s SimpleFFT Tutorial I implemented a method computeTrueFreqs(float* inputOutputData). The input of this method is the result from forwardFFT.performRealOnlyForwardTransform without any overlap. The input signal is one sine wave e.g. frequency 2416.03 Hz with Sampel rate 44.1 KHz, fft size 2048. As true frequncy of bin 112 I get 2420.33 Hz with delta phase 2.51 Radian. The right delta phase would have been 1.25. As far as I could figure out it seems to be the phase diff between adjacent fft frames which is wrong but I cannot see what I do wrong in my code. I would be grateful if someone could give me a hint?

void computeTrueFreqs(float* inputOutputData)
{
    for (int i = 0; i < fftSize; ++i)
    {
        //magnitude
        std::complex<float> z(inputOutputData[i], inputOutputData[i+1]);
        float zm = std::abs(z);

        //check minimum magnitude
        if (zm > magTresh)
        {
            phase[count] = std::arg(z);

            //compute phase diff for each bin between adjacent fft frames
            float deltaPhase = (phase[count] - phaseLast[count]);

            // Wrap delta phase into [-Pi, Pi) interval
            float deltaPhase = 0;
            deltaPhase = deltaPhase - MathConstants<float>::twoPi * floor(deltaPhase / MathConstants<float>::twoPi + .5f);

            //add frequency shift to the center frequency of each bin
            freqsCorrected[count] = binFreqs[count] + deltaPhase * sr_2PI_WindowSize; //where binFreqs[i] = (float)i * (float)sampleRate / (float)fftSize and
            //sr_2PI_WindowSize = SampleRate / (MathConstants<float>::twoPi * fftSize);

            phaseLast[count] = phase[count];
        }
    }
}

After many tests I figured the following code variation of Simpel fft tutorial out which seems to do the true frequency estimation right. I acutally did not quite understand why the before mentioned errors do not anymore occure but I suppose cleaning up most of the unnecessary parts (even visuals) helped to rerarrange the code more efficiently and possibly to avoid delays in the audio data-architecture. So I give my solution (SimpleFFTTutorial_01.h) here for notice. The following tutorial was very helpful for understanding the concept of phase differentiation: http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

/*
  ==============================================================================

   This file is part of the JUCE tutorials.
   Copyright (c) 2020 - Raw Material Software Limited

   The code included in this file is provided under the terms of the ISC license
   http://www.isc.org/downloads/software-support-policy/isc-license. Permission
   To use, copy, modify, and/or distribute this software for any purpose with or
   without fee is hereby granted provided that the above copyright notice and
   this permission notice appear in all copies.

   THE SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES,
   WHETHER EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR
   PURPOSE, ARE DISCLAIMED.

  ==============================================================================
*/

/*******************************************************************************
 The block below describes the properties of this PIP. A PIP is a short snippet
 of code that can be read by the Projucer and used to generate a JUCE project.

 BEGIN_JUCE_PIP_METADATA

 name:             SimpleFFTTutorial
 version:          1.0.0
 vendor:           JUCE
 website:          http://juce.com
 description:      Displays an FFT spectrogram.

 dependencies:     juce_audio_basics, juce_audio_devices, juce_audio_formats,
                   juce_audio_processors, juce_audio_utils, juce_core,
                   juce_data_structures, juce_dsp, juce_events, juce_graphics,
                   juce_gui_basics, juce_gui_extra
 exporters:        xcode_mac, vs2019, linux_make

 type:             Component
 mainClass:        SpectrogramComponent

 useLocalCopy:     1

 END_JUCE_PIP_METADATA

*******************************************************************************/


#pragma once


//==============================================================================
class SpectrogramComponent   : public juce::AudioAppComponent,  //AnalyserComponent
                               private juce::Timer
{
public:
    SpectrogramComponent()
        : forwardFFT(fftOrder)
    {        
        setAudioChannels (2, 0);  // we want a couple of input channels but no outputs
        startTimerHz (60);
        setSize (700, 500);
    }

    ~SpectrogramComponent() override
    {
        shutdownAudio();
    }

    //==============================================================================
    void prepareToPlay (int, double sampleRate) override
    {
        //precompute Bin Freqs    
        for (int i = 0; i < fftSize; ++i)
        {
            binFreqs[i] = (float)i * (float)sampleRate / (float)fftSize;
        }

        sr_2PI_WindowSize = (float)sampleRate / (MathConstants<float>::twoPi * (float)fftSize);        
        
        for (auto& v : phase)
        {
            v = 0;
        }
        for (auto& v : phaseLast)
        {
            v = 0;
        }

    }

    void releaseResources() override          {}

    void getNextAudioBlock (const juce::AudioSourceChannelInfo& bufferToFill) override 
    {
        if (bufferToFill.buffer->getNumChannels() > 0)
        {
            auto* channelData = bufferToFill.buffer->getReadPointer(0, bufferToFill.startSample);

            for (auto i = 0; i < bufferToFill.numSamples; ++i)
                pushNextSampleIntoFifo(channelData[i]);
        }
    }

    //==============================================================================
    void paint (juce::Graphics& g) override
    {
        
    }

    void timerCallback() override 
    {
        if (nextFFTBlockReady)
        {
            forwardFFT.performRealOnlyForwardTransform(fftData.data());
            computeTrueFreqs(fftData.data());
            nextFFTBlockReady = false;
            //repaint();
        }
    }

    void pushNextSampleIntoFifo (float sample) noexcept
    {
        // if the fifo contains enough data, set a flag to say
        // that the next line should now be rendered..
        if (fifoIndex == fftSize)       // [8]
        {
            if (!nextFFTBlockReady)    // [9]
            {
                std::fill(fftData.begin(), fftData.end(), 0.0f);

                //Windowing
                dsp::WindowingFunction <float>(fftSize, dsp::WindowingFunction<float>::hamming, false, 0.0f).multiplyWithWindowingTable(&fifo[0], fftSize);

                std::copy(fifo.begin(), fifo.end(), fftData.begin());
                nextFFTBlockReady = true;
            }

            fifoIndex = 0;
        }

        fifo[(size_t)fifoIndex++] = sample; // [9]
    }   
    
    void computeTrueFreqs(float* inputOutputData)    {        
        //clear old freq values
        std::fill(freqsCorrected.begin(), freqsCorrected.end(), 0);              
        for (int i = 0; i < fftSize; ++i)
        {
            //magnitude berechnen
            std::complex<float> z(inputOutputData[i*2], inputOutputData[i*2+1]);
            float zm = std::abs(z);
            
            //the following procedure can be under the condition of a certain minimum amplitude/magnitude
            if (zm > mag_Tresh)
            {                  
                phase[i] = std::arg(z) ;              

                //compute phase diff for each bin between adjacent fft frames               
                float deltaPhase = (phase[i] - phaseLast[i]);
                phaseLast[i] = phase[i];
                
                // Wrap delta phase into [-Pi, Pi) interval //phase must be inside [pi -pi]
                float deltaPhase_ = 0;
                deltaPhase_ = deltaPhase - MathConstants<float>::twoPi * floor(deltaPhase / MathConstants<float>::twoPi + .5f);    

                //compute centerfrequnecy of bin[i]                
                binFreqs[i] = (float)i * 44100.f / (float)fftSize;

                //add frequency shift to the center frequency of each bin
                freqsCorrected[i] = binFreqs[i] + deltaPhase_ * sr_2PI_WindowSize; 

                float x = 0;
            }         
        }      
    }   
   

    /*[1]: The FFT order designates the size of the FFT window and the number of points on which it will operate corresponds to 2 to the power of the order. In this case, let's use an order of 10 which will produce an FFT with 2 ^ 11 = 2048 points.
      [2]: To calculate the corresponding FFT size, we use the left bit shift operator which produces 2048 as binary number 100000000000.*/
    static constexpr auto fftOrder = 11;               // [1]
    static constexpr auto fftSize = 1 << fftOrder;     // [2]

private:
    /*[3]: Declare a dsp::FFT object to perform the forward FFT on.
    [4]: The fifo float array of size 2048 will contain our incoming audio data in samples.
    [5]: The fftData float array of size 2048 will contain the results of our FFT calculations.
    [6]: This temporary index keeps count of the amount of samples in the fifo.
    [7]: This temporary boolean tells us whether the next FFT block is ready to be rendered.*/
    juce::dsp::FFT forwardFFT;                          // [3]
    juce::Image spectrogramImage;

    std::array<float, fftSize> fifo;                    // [4]
    std::array<float, fftSize * 2> fftData;             // [5]
    int fifoIndex = 0;                                  // [6]
    bool nextFFTBlockReady = false;                     // [7]

    std::array<float, fftSize> binFreqs;
    std::array<float, fftSize> phase;
    std::array<float, fftSize> phaseLast;             
    float sr_2PI_WindowSize = 44100.0f / (MathConstants<float>::twoPi * fftSize);
    std::array<float, fftSize> freqsCorrected;
    float mag_Tresh = 2;
    int sizeDoubled = (int)fftSize * 2;

    JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (SpectrogramComponent)
};

1 Like