CODEWITHSUNDEEP
c++matlabaudiofft
Phlox Midas
Phlox Midas

Reputation: 4349

Getting values for specific frequencies in a short time fourier transform

I'm trying to use C++ to recreate the spectrogram function used by Matlab. The function uses a Short Time Fourier Transform (STFT). I found some C++ code here that performs a STFT. The code seems to work perfectly for all frequencies but I only want a few. I found this post for a similar question with the following answer:

Just take the inner product of your data with a complex exponential at the frequency of interest. If g is your data, then just substitute for f the value of the frequency you want (e.g., 1, 3, 10, ...)

Having no background in mathematics, I can't figure out how to do this. The inner product part seems simple enough from the Wikipedia page but I have absolutely no idea what he means by (with regard to the formula for a DFT)

a complex exponential at frequency of interest

Could someone explain how I might be able to do this? My data structure after the STFT is a matrix filled with complex numbers. I just don't know how to extract my desired frequencies.

Relevant function, where window is Hamming, and vector of desired frequencies isn't yet an input because I don't know what to do with them:

Matrix<complex<double>> ShortTimeFourierTransform::Calculate(const vector<double> &signal,
    const vector<double> &window, int windowSize, int hopSize)
{
    int signalLength = signal.size();
    int nOverlap = hopSize;
    int cols = (signal.size() - nOverlap) / (windowSize - nOverlap);
    Matrix<complex<double>> results(window.size(), cols);

    int chunkPosition = 0;
    int readIndex;
    // Should we stop reading in chunks? 
    bool shouldStop = false;
    int numChunksCompleted = 0;
    int i;
    // Process each chunk of the signal
    while (chunkPosition < signalLength && !shouldStop)
    {
        // Copy the chunk into our buffer
        for (i = 0; i < windowSize; i++)
        {
            readIndex = chunkPosition + i;
            if (readIndex < signalLength)
            {
                // Note the windowing! 
                data[i][0] = signal[readIndex] * window[i];
                data[i][1] = 0.0;
            }
            else
            {
                // we have read beyond the signal, so zero-pad it!
                data[i][0] = 0.0;
                data[i][1] = 0.0;
                shouldStop = true;
            }
        }
        // Perform the FFT on our chunk
        fftw_execute(plan_forward);

        // Copy the first (windowSize/2 + 1) data points into your spectrogram.
        // We do this because the FFT output is mirrored about the nyquist 
        // frequency, so the second half of the data is redundant. This is how
        // Matlab's spectrogram routine works.
        for (i = 0; i < windowSize / 2 + 1; i++)
        {               
            double real = fft_result[i][0];
            double imaginary = fft_result[i][1];
            results(i, numChunksCompleted) = complex<double>(real, imaginary);
        }
        chunkPosition += hopSize;
        numChunksCompleted++;
    } // Excuse the formatting, the while ends here.
    return results;
}

Upvotes: 1

Views: 3540

Answers (1)

hotpaw2
hotpaw2

Reputation: 70693

Look up the Goertzel algorithm or filter for example code that uses the computational equivalent of an inner product against a complex exponential to measure the presence or magnitude of a specific stationary sinusoidal frequency in a signal. Performance or resolution will depend on the length of the filter and your signal.

Upvotes: 1

Related Questions