Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Timescale FFT Issues
```
#include <Bela.h>
#include <libraries/Fft/Fft.h>
#include <libraries/Scope/Scope.h>
#include <libraries/Gui/Gui.h>
#include <libraries/GuiController/GuiController.h>
#include <cmath>
#include <cstring>
#include <vector>
#include <algorithm>
#include "MonoFilePlayer.h"
// ---- NOTES TO SELF ---- //
//==========================================//
// - outputRead and Write buffers and hop size attention for time scale changes
// - record and playback make separate functions to call in render -- attention to setup of both
// ---- DECLARATIONS ---- //
//====================================//
// FFT-related variables
Fft gFft; // FFT processing object
const int gFftSize = 1024; // FFT window size in samples
float gTimeScale = 2.0; // Scale time and phase difference
const int gHopSize = 128; // How often we calculate a window -- analyis hopsize
int gSynthHopSize = gHopSize * gTimeScale; // synthesis hopSize
float gScaleFactor = 0.5; // How much to scale the output, based on window type and overlap
float gPitchShift = 1.0; // Ratio of output to input frequency
// Window of Samples Circular Buffer
const int gBufferSize = 16384; // arbitrary large number for maximum audio data needed
std::vector<float> gInputBuffer;
int gInputBufferPointer = 0;
int gHopCounter = 0;
// Overlap-add Output Circular Buffer
std::vector<float> gOutputBuffer;
// Start the write pointer ahead of the read pointer by at least window + hop, with some margin
int gOutputBufferWritePointer = gFftSize + 2*gHopSize; // NOPE
int gOutputBufferReadPointer = 0;
// FFT Window Buffers
std::vector<float> gAnalysisWindowBuffer;
std::vector<float> gSynthesisWindowBuffer;
// Name of the sound file (in project folder) -- needs modified for real-time recording
std::string gFilename = "voice.wav";
// Play mono sound from buffer
MonoFilePlayer gPlayer;
// Thread for FFT processing -- lower priority
AuxiliaryTask gFftTask;
int gCachedInputBufferPointer = 0;
void process_fft_background(void *);
// Declare Scope + GUI
Scope gScope;
Gui gGui;
GuiController gGuiController;
// ---- SETUP ---- //
//====================================//
bool setup(BelaContext *context, void *userData)
{
// Load audio file + check for error
if(!gPlayer.setup(gFilename)) {
rt_printf("Error loading audio file '%s'\n", gFilename.c_str());
return false;
}
// Print audio file info
rt_printf("Loaded the audio file '%s' with %d frames (%.1f seconds)\n",
gFilename.c_str(), gPlayer.size(),
gPlayer.size() / context->audioSampleRate);
// Set up FFT and its buffers
gFft.setup(gFftSize);
gInputBuffer.resize(gBufferSize);
gOutputBuffer.resize(gBufferSize);
// Calculate FFT windows
gAnalysisWindowBuffer.resize(gFftSize);
gSynthesisWindowBuffer.resize(gFftSize);
for(int n = 0; n < gFftSize; n++) {
// Generate Hann window and copy to both window buffers -- TRY other windows
gAnalysisWindowBuffer[n] = 0.5f * (1.0f - cosf(2.0 * M_PI * n / (float)(gFftSize - 1)));
gSynthesisWindowBuffer[n] = gAnalysisWindowBuffer[n];
}
// Init the oscilloscope
gScope.setup(2, context->audioSampleRate);
// Set up GUI
gGui.setup(context->projectName);
gGuiController.setup(&gGui, "Pitch Shift Controller");
// Arguments: name, default value, minimum, maximum, increment
gGuiController.addSlider("Shift", 0, -12, 12, 0);
// Set up the thread for the FFT
gFftTask = Bela_createAuxiliaryTask(process_fft_background, 50, "bela-process-fft");
return true;
}
// ---- WRAP PHASE FUNCTION ---- //
//==========================================//
// Wrap the phase to the range -pi to pi
float wrapPhase(float phaseIn)
{
if (phaseIn >= 0)
return fmodf(phaseIn + M_PI, 2.0 * M_PI) - M_PI;
else
return fmodf(phaseIn - M_PI, -2.0 * M_PI) + M_PI;
}
// ---- FFT PROCESS FUNCTION ---- //
//==========================================//
void process_fft(std::vector<float> const& inBuffer, unsigned int inPointer, std::vector<float>& outBuffer, unsigned int outPointer)
{
static std::vector<float> unwrappedBuffer(gFftSize); // Container to hold the unwrapped time-domain values
static std::vector<float> lastInputPhases(gFftSize); // Hold the phases from the previous hop of input signal
static std::vector<float> lastOutputPhases(gFftSize); // and output (synthesised) signal
// Hold data that has been converted from magnitude-phase to magnitude-frequency -- Pitch Shifting purposes
static std::vector<float> analysisMagnitudes(gFftSize / 2 + 1);
static std::vector<float> analysisFrequencies(gFftSize / 2 + 1);
static std::vector<float> synthesisMagnitudes(gFftSize / 2 + 1);
static std::vector<float> synthesisFrequencies(gFftSize / 2 + 1);
// Copy inBuffer into FFT input
for(int n = 0; n < gFftSize; n++) {
// Calculate index of inBuffer, apply window and store in unwrappedBuffer for FFT processing
int circularBufferIndex = (inPointer + n - gFftSize + gBufferSize) % gBufferSize;
unwrappedBuffer[n] = inBuffer[circularBufferIndex] * gAnalysisWindowBuffer[n];
}
// Run the FFT based on the time domain input
gFft.fft(unwrappedBuffer);
// ---- FFT ANALYSIS ---- //
//==========================================//
// Analyse lower half of the spectrum -- upper half = complex conjugate (identical apart from sign)
for(int n = 0; n <= gFftSize/2; n++) {
// Convert real and imaginary components to amplitude and phase values
float amplitude = gFft.fda(n);
float phase = atan2f(gFft.fdi(n), gFft.fdr(n));
// Calculate exact frequency by taking the phase difference for this bin between the last and current hop -- Causal Process
float phaseDiff = phase - lastInputPhases[n];
// Subtract the amount of phase increment we'd expect to see based
// on the centre frequency of this bin (2*pi*n/gFftSize) for this
// hop size, then wrap to the range -pi to pi
float binCentreFrequency = 2.0 * M_PI * (float)n / (float)gFftSize;
phaseDiff = wrapPhase(phaseDiff - binCentreFrequency * gHopSize);
// Find deviation in (fractional) number of bins from the centre frequency
float binDeviation = phaseDiff * (float)gFftSize / (float)gHopSize / (2.0 * M_PI);
// Add the original bin number to get the fractional bin where this partial belongs
analysisFrequencies[n] = (float)n + binDeviation;
// Save the magnitude for later and for the GUI
analysisMagnitudes[n] = amplitude;
// Save the phase for next hop
lastInputPhases[n] = phase;
}
// ---- FFT SYNTHESIS ---- //
//==========================================//
// Zero out the synthesis bins, ready for new data
for(int n = 0; n <= gFftSize/2; n++) {
synthesisMagnitudes[n] = synthesisFrequencies[n] = 0;
}
// remove pitch shifting for now by simply copying data from anaylsis process
for(int n = 0; n <= gFftSize/2; n++) {
synthesisMagnitudes[n] = analysisMagnitudes[n];
synthesisFrequencies[n] = analysisFrequencies[n];
}
//PITCH SHIFT PROCESSING -- plan to reincorporate later on
// // Handle the pitch shift, storing frequencies into new bins
// for(int n = 0; n <= gFftSize/2; n++) {
// // TODO: find the nearest bin to the shifted frequency
// int newBin = floorf(n * gPitchShift + 0.5);
// // Ignore any bins that have shifted above Nyquist
// if(newBin <= gFftSize / 2) {
// synthesisMagnitudes[newBin] += analysisMagnitudes[n];
// // Scale the frequency by the pitch shift ratio
// synthesisFrequencies[newBin] = analysisFrequencies[n] * gPitchShift;
// }
// }
// Synthesise frequencies into new magnitude and phase values for FFT bins
for(int n = 0; n <= gFftSize / 2; n++) {
float amplitude = synthesisMagnitudes[n];
// Get the fractional offset from the bin centre frequency
float binDeviation = synthesisFrequencies[n] - n;
// Multiply to get back to a phase value
float phaseDiff = binDeviation * 2.0 * M_PI * (float)gSynthHopSize / (float)gFftSize;
// Add the expected phase increment based on the bin centre frequency
float binCentreFrequency = 2.0 * M_PI * (float)n / (float)gFftSize;
phaseDiff += binCentreFrequency * gSynthHopSize;
// Advance the phase from the previous hop
float phase = wrapPhase(lastOutputPhases[n] + phaseDiff);
// float outPhase = inPhase * gTimeScale;
// Now convert magnitude and phase back to real and imaginary components
gFft.fdr(n) = amplitude * cosf(phase);
gFft.fdi(n) = amplitude * sinf(phase);
// Store the complex conjugate in the upper half of the spectrum
if(n > 0 && n < gFftSize / 2) {
gFft.fdr(gFftSize - n) = gFft.fdr(n);
gFft.fdi(gFftSize - n) = -gFft.fdi(n);
}
// Save the phase for the next hop
lastOutputPhases[n] = phase;
}
// Run the inverse FFT
gFft.ifft();
// Add timeDomainOut into the output buffer
for(int n = 0; n < gFftSize; n++) {
int circularBufferIndex = (outPointer + n - gFftSize + gBufferSize) % gBufferSize;
// circularBufferIndex = circularBufferIndex * gTimeScale;
outBuffer[circularBufferIndex] += gFft.td(n) * gSynthesisWindowBuffer[n];
}
}
// ---- CALL FFT PROCESS ---- //
//==========================================//
void process_fft_background(void *)
{
process_fft(gInputBuffer, gCachedInputBufferPointer, gOutputBuffer, gOutputBufferWritePointer);
// Update the output buffer write pointer to start at the next hop
gOutputBufferWritePointer = (gOutputBufferWritePointer + gHopSize) % gBufferSize;
}
// ---- MAIN FUNCTION ---- //
//==========================================//
void render(BelaContext *context, void *userData)
{
// Get the pitch shift in semitones from the GUI slider and convert to ratio
float pitchShiftSemitones = gGuiController.getSliderValue(0);
gPitchShift = powf(2.0, pitchShiftSemitones / 12.0);
for(unsigned int n = 0; n < context->audioFrames; n++) {
// Read the next sample from the buffer
float in = gPlayer.process();
// Store the input sample in a buffer for the FFT
// Increment the pointer and when full window has been
// assembled, call process_fft()
gInputBuffer[gInputBufferPointer++] = in;
if(gInputBufferPointer >= gBufferSize) {
// Wrap the circular buffer
// Notice: this is not the condition for starting a new FFT
gInputBufferPointer = 0;
}
// Get the output sample from the output buffer
float out = gOutputBuffer[gOutputBufferReadPointer];
// Then clear the output sample in the buffer so it is ready for the next overlap-add
gOutputBuffer[gOutputBufferReadPointer] = 0;
// Scale the output down by the scale factor, compensating for the overlap
out *= gScaleFactor;
// Increment the read pointer in the output cicular buffer
gOutputBufferReadPointer++;
if(gOutputBufferReadPointer >= gBufferSize)
gOutputBufferReadPointer = 0;
// Increment the hop counter and start a new FFT if we've reached the hop size -- WHICH hop size?
if(++gHopCounter >= gHopSize) {
gHopCounter = 0;
gCachedInputBufferPointer = gInputBufferPointer;
Bela_scheduleAuxiliaryTask(gFftTask);
}
// Write the audio to the output
for(unsigned int channel = 0; channel < context->audioOutChannels; channel++) {
audioWrite(context, n, channel, out);
}
// Log to the Scope
gScope.log(in, out);
}
}
void cleanup(BelaContext *context, void *userData)
{
}
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment