FFT Frequency Shifter Extravaganza I've been writing a frequency shifter in C++ for an Audio Unit plugin. The best way for me to do this while keeping a large amount of other algorithms working in the same plugin was by writing a separate header file with a full class definition that could handle all things related to the frequency shift. After a few problems with my code were cleared up in this question, I decided to go look for some optimizations and better functionality. I didn't get very far, so here we are. For any people answering my question, I ask that you address the following:


Please keep the arguments to the processBlock function as they are. I am not good enough to in any way revise how audio is supplied to the processing function.
If you've got a better FFT algorithm (you probably do), please implement it! I used an elementary function from Rosetta Code that didn't even utilize all the fancy bitshifting/reversing that it could be. I also know that Apple offers various vDSP functions for FFT, but I am even more clueless on how to implement those.
Of course, optimizations would be nice since the function is called several times a second.
I always appreciate a great lecture on code conventions, especially in C++.


Thanks in advance for any answers. The following is from FrequencyShifter.h:

#include &lt;complex&gt;
#include &lt;valarray&gt;

typedef std::complex&lt;double&gt; Complex;
typedef std::valarray&lt;Complex&gt; CArray;

class FrequencyShifter {
float sampleRate;
public:
    FrequencyShifter() {

    }
    void setSampleRate(float inSampleRate) {
        sampleRate = inSampleRate;
    }
    void fft(CArray&amp; x)
    {
        const size_t N = x.size();
        if (N &lt;= 1) return;

        // divide
        CArray even = x[std::slice(0, N/2, 2)];
        CArray  odd = x[std::slice(1, N/2, 2)];

        // conquer
        fft(even);
        fft(odd);

        // combine
        for (size_t k = 0; k &lt; N/2; ++k)
        {
            Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
            x[k    ] = even[k] + t;
            x[k+N/2] = even[k] - t;
        }
    }
    void processBlock(float *inBlock, const int inFramesToProcess, float scale) {
        //inFramesToProcess is the amount of samples in inBlock
        Complex *copy = new Complex[inFramesToProcess];
        for (int frame = 0; frame&lt;inFramesToProcess; frame++) {
            copy[frame] = Complex((double)inBlock[frame], 0.0);
        }
        CArray data(copy, inFramesToProcess);
        fft(data);
        const float freqoffsets = 2.0 * PI/ inFramesToProcess;
        const float normfactor  = 2.0 / inFramesToProcess;
        for (int frame = 0; frame&lt;inFramesToProcess; frame++) {
            inBlock[frame] = 0.5 * data[0].real();
            for (int x = 1; x &lt; data.size() / 2; x++) {
                float arg = freqoffsets * x *f rame*scale;
                inBlock[frame] += data[x].real()*cos(arg) - data[x].imag()*sin(arg);
            }
            inBlock[frame] *= normfactor;
        }
    }
};
The main thing that I feel you should change is the heap allocations inside the processBlock(). With realtime processing like this audio application you never want to call new/delete within time critical code, either yourself (or indirectly like adding a vector, vector.resize() etc. etc.) The reason it's unpredictable how long a new/delete will take.

You have two lines that call new either directly/indirectly

Complex *copy = new Complex[inFramesToProcess];
CArray data(copy, inFramesToProcess); // indirect call of new


The solution is to create work areas beforehand as a member variable. This is created to be the right size at setup time before playback then the processBlock() can write to it as it wants.

So data (now renamed fftData) is outside as a member variable...



Edit: better solution below!



#include &lt;complex&gt;
#include &lt;valarray&gt;

typedef std::complex&lt;double&gt; Complex;
typedef std::valarray&lt;Complex&gt; CArray;

class FrequencyShifter {

float sampleRate;
const int guessInFramesToProcess;
CArray fftData;

public:
    FrequencyShifter()
        : guessInFramesToProcess {512}, // best guess at what the num samples per block will be
                                        // or is there a prepareToPlay function that you can
                                        // send in the blocksize before user hits play?
          fftData {Complex(0.0, 0.0), inFramesToProcess}
    {
    }
    void setSampleRate(float inSampleRate) { // something like this but for blocksize???
        sampleRate = inSampleRate;
    }
    void fft(CArray&amp; x)
    {
        const size_t N = x.size();
        if (N &lt;= 1) return;

        // divide
        CArray even = x[std::slice(0, N/2, 2)];
        CArray  odd = x[std::slice(1, N/2, 2)];

        // conquer
        fft(even);
        fft(odd);

        // combine
        for (size_t k = 0; k &lt; N/2; ++k)
        {
            Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
            x[k    ] = even[k] + t;
            x[k+N/2] = even[k] - t;
        }
    }
    void processBlock(float *inBlock, const int inFramesToProcess, float scale) {

        /* now there are no heap allocations in the time-critical processBlock */

        fillFftData(inBlock, inFramesToProcess, fftData);       // function that declares intent

        fft(fftData);

        const float freqoffsets = 2.0 * PI/ inFramesToProcess;  // below is quite hard to figure out
        const float normfactor  = 2.0 / inFramesToProcess;      // your intent. maybe functions
                                                                // like below?
        for (int frame = 0; frame&lt;inFramesToProcess; frame++) {
            inBlock[frame] = 0.5 * fftData[0].real();
            for (int x = 1; x &lt; fftData.size() / 2; x++) {
                float arg = freqoffsets * x *frame*scale;
                inBlock[frame] += fftData[x].real()*cos(arg) - fftData[x].imag()*sin(arg);
            }
            inBlock[frame] *= normfactor;
        }
    }
private:
    // probably good to make little functions like this in your processBlock
    // so you can easier debug later. here's an example
    void fillFftData(float* inBlock, const int inFramesToProcess, CArray&amp; data) {
        for (int frame = 0; frame&lt;inFramesToProcess; frame++) {
            data[frame] = Complex(static_cast&lt;double&gt; (inBlock[frame]), 0.0);
        }                       // static cast ugly but preferred to c-style casts
    }                           // for better compiler checking
};


Now the issue is with how do you tell in advance what the blocksize is going to be i.e. inFramesToProcess ?

The way the code stands now, it may not work because we need some way to let the host tell us the blocksize before we actually start processing the audio in processBlock().

We need this because now the fftData needs to be fully resized the same as the blocksize (and filled with zeros) before during the setup time.

Is there a function that you can hook into that will tell your plugin what the expected blocksize will be?

(Sorry I'm a JUCE user I don't know the specifics of what framework you're using but for example in JUCE there's a function called prepareToPlay(int sampleRate, int expectedBlockSize) that gets called before processing the audio with processBlock(), this would be a perfect place to resize and zero fftData).

Hopefully you can figure this bit out!



Some smaller points


Not sure you need to copy into the buffer then into the fftData? Seems unnecessary. (If you need it, again, allocate the buffer beforehand).
If you split your code into small functions with good names that state the purpose, then it's easier to come back to your code in a few months time :)
Try and avoid C-style casts in C++ as the compiler 'knows' what's going on better with a static_cast.




The main point is really the heap allocations in the processBlock(). Only after it might be time to look at implementing and benchmarking the FFT performance vs. Apple Accelerate vDSP :)

There are other things you should avoid in time critical code too, for more info see this section of this talk...

CppCon 2015: Timur Doumler “C++ in the Audio Industry”



EDIT: Just overcompensate the FFT data buffer.



Sorry I wasn't thinking clearly about the issue previously! The simplest way around this issue is to make fftData larger than you could possibly need i.e. larger than the largest possible block size, then if a smaller block size is actually dealt with in processBlock(), it will only access the start of the buffer. Note how we are now sending the block size through to fft() and now iterating over currentInFramesToProcess (since fftData.size() is overcompensated). 

#include &lt;complex&gt;
#include &lt;valarray&gt;

typedef std::complex&lt;double&gt; Complex;
typedef std::valarray&lt;Complex&gt; CArray;

class FrequencyShifter {

float sampleRate;
const int largestPossibleInFramesToProcess;
CArray fftData;

public:
    FrequencyShifter()
        : largestPossibleInFramesToProcess {2048}, // best guess, could add more for safety
          fftData {Complex(0.0, 0.0), largestPossibleInFramesToProcess}
    {
    }
    void setSampleRate(float inSampleRate) {
        sampleRate = inSampleRate;
    }
    void fft(CArray&amp; x, int currentInFramesToProcess)
    {
        const size_t N = currentInFramesToProcess;
        if (N &lt;= 1) return;

        // divide
        CArray even = x[std::slice(0, N/2, 2)];
        CArray  odd = x[std::slice(1, N/2, 2)];

        // conquer
        fft(even);
        fft(odd);

        // combine
        for (size_t k = 0; k &lt; N/2; ++k)
        {
            Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
            x[k    ] = even[k] + t;
            x[k+N/2] = even[k] - t;
        }
    }
    void processBlock(float *inBlock, const int inFramesToProcess, float scale) {

        /* now there are no heap allocations in the time-critical processBlock */

        fillFftData(inBlock, inFramesToProcess, fftData);       // function that declares intent

        fft(fftData, inFramesToProcess);

        const float freqoffsets = 2.0 * PI/ inFramesToProcess;  // below is quite hard to figure out
        const float normfactor  = 2.0 / inFramesToProcess;      // your intent. maybe functions
                                                                // like below?
        for (int frame = 0; frame &lt; inFramesToProcess; frame++) {
            inBlock[frame] = 0.5 * fftData[0].real();
            for (int x = 1; x &lt; inFramesToProcess / 2; x++) {
                float arg = freqoffsets * x * frame * scale;
                inBlock[frame] += fftData[x].real()*cos(arg) - fftData[x].imag()*sin(arg);
            }
            inBlock[frame] *= normfactor;
        }
    }
private:
    // probably good to make little functions like this in your processBlock
    // so you can easier debug later. here's an example
    void fillFftData(float* inBlock, const int inFramesToProcess, CArray&amp; data) {
        for (int frame = 0; frame &lt; inFramesToProcess; frame++) {
            data[frame] = Complex(static_cast&lt;double&gt; (inBlock[frame]), 0.0);
        }                       // static cast ugly but preferred to c-style casts
    }                           // for better compiler checking
};


This would get you by for now.

The 'best practise' would be to make your FFT size independent of whatever the processBlock() size is. That way if the user changes the soundcard block settings, they will still hear the same thing. i.e. lock FFT to size 1024, soundcard block size processBlock() variable 64-512

This is accomplished with a circular buffer of size 1024, you fill this up sample by sample with input, then when it's full you call fft() on those 1024 samples.

Hope these are all helpful tips to get started!