# 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
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.

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

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

#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)
{

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

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

}

void timerCallback() override
{
{
forwardFFT.performRealOnlyForwardTransform(fftData.data());
computeTrueFreqs(fftData.data());
//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]
{
{
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());
}

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