A small IFFT addition to bela/Fft library

Tobias

Hi!
I am working on a project where I would like to use the usual FFT/IFFT methodology to transform an audio signal to frequency domain, do some filtering on the magnitude spectrum, recombine with original phase and go back to the time domain.
However, I found the (otherwise very useful) Fft library lacking, since there is no direct way to do an IFFT on an arbitrary ne10_fft_cpx_float32_t* frequencyDomain object. So I have made a small proposal for a function that does just that:

// Header
void ifft(const std::vector<float>& reInput, const std::vector<float>& imInput);

// Source
void Fft::ifft(const std::vector<float>& reInput, const std::vector<float>& imInput)
{
	for (int i = 0; i < length/2 + 1; i++) {
		frequencyDomain[i].r = reInput[i];
		frequencyDomain[i].i = imInput[i];
	}
	ifft();
}

It should be perfectly compatible with the existing library. Using a loop is probably not the most efficient way, but I could not make a memcpy version work. It assumes that ne10_fft_cpx_float32_t* frequencyDomain has (NFFT / 2) + 1 elements (<https://projectne10.github.io/Ne10/doc/groupR2CFFT__IFFT.html>).
What do you think? And feel free to add it to the standard library if you'd like.

giuliomoro

you won't be able to do that with memcpy(), because the ne10_fft_cpx_float32_t array that is pointed to by frequencyDomain is effectively an interleaved array, so the loop is needed to interleave the two non-interleaved input vectors. You may want to add a check for the length of the input arrays before trying to access their elements. It's will be a relatively very cheap check compared to the CPU load of the whole method.
You could use fdr() and fdi() to access the elements. It's going to compile to exactly the same assembly code, but it will be slightly more maintainable if the internal representation of frequencyDomain ever changes (e.g.: in case we move away from ne10 - which I am not planning to do anytime soon).

Tobias

Yeah that's what I thought too after reading about the interleaved array in the documentation. Definitely needs some safety checks.
But the reason I shared this here in the first palce was just as an idea for a small but maybe useful addition to the library since doing an IFFT now requires accessing private memebers of the Fft object which kind of defeats the purpose of a wrapper 😉
Nothing more to see here!

giuliomoro

Tobias useful addition to the library

it is a useful addition: a nice complement to fft(const std::vector<float>& input);.

Tobias now requires accessing private memebers of the Fft object

it doesn't. The members are private, but they have accessor methods (fdr() and fdi()) that make them available publicly. So one could currently write:

	for (int i = 0; i < length/2 + 1; i++) {
		myFft.fdr(i) = reInput[i];
		myFft.fdi(i) = imInput[i];
	}
	myFftifft();

a revised method, which checks the input size and uses fdr() and fdi() could be like this:

void Fft::ifft(const std::vector<float>& reInput, const std::vector<float>& imInput)
{
	if(reInput.size() < length / 2 + 1 || imInput.size() < length / 2 + 1)
		return;
	for (int i = 0; i < length/2 + 1; i++) {
		fdr(i) = reInput[i];
		fdi(i) = imInput[i];
	}
	ifft();
}

lobo

giuliomoro I am currently working on autocorrelation using FFT and trying to emulate the code (https://github.com/yoyololicon/ape-examples/blob/main/yin.hpp) below using Bela

        std::copy(cyclic_begin(history, history_counter),
			  cyclic_end(history, history_counter, history.size()),
			  fft_buf1.begin());
	std::fill(fft_buf1.begin() + 2 * W - 1, fft_buf1.end(), 0);

	std::copy(fft_buf1.begin(), fft_buf1.begin() + W, fft_buf2.begin());
	std::fill(fft_buf2.begin() + W, fft_buf2.end(), 0);

	fft.forward(fft_buf1);
	fft.forward(fft_buf2);

	for (i = 0; i < N / 2 + 1; i++)
		fft_buf1[i] *= std::conj(fft_buf2[i]);

	for (; i < N; i++)
		fft_buf1[i] = std::conj(fft_buf1[N - i]);

	fft.inverse(fft_buf1);

where basically fft_buf1 and fft_buf2 respectively load up with 2W and W samples (roughly) of audio and perform FFT with same number of size 2W. fft_buf2 is zero-padded.

The code then realises and rearranges the obtained coefficients followed by IFFT fft.inverse() to obtain autocorrelation.

Currently I am having

static std::vector<float> unwrapBuffer1(gFftSize, 0.0);
static std::vector<float> unwrapBuffer2(gFftSize, 0.0);

for (int k = 0; k < gFftSize; k++) {
	int index_circular = (k + gWritePointer - gFftSize + gBufferSize) % gBufferSize;
	unwrapBuffer1[k] = gCircularBuffer[index_circular];
}
for (int k = 0; k < gInputSize; k++) {
	unwrapBuffer2[k] = unwrapBuffer1[k];
}
gFft1.fft(unwrapBuffer1);
gFft2.fft(unwrapBuffer2);

std::vector<float> spectrum(gFftSize / 2 + 1);
for (int i = 0; i < gFftSize / 2 + 1; i++) {
	spectrum[i] = gFft1.fda(i) * gFft2.fda(i) + gFft1.fdi(i) * gFft2.fdi(i);
}
for (int i = gFftSize / 2; i < gFftSize; i++) {
	spectrum[i] = spectrum[gFftSize - i];
}
gFft1.ifft(spectrum);  // this line won't work

Because I can't call ifft() this way, I am wondering how to translate the part of IFFT to Bela properly.

giuliomoro

Does your spectrum only have a real part (i.e.: no imaginary part) ?
In that case, you can add the method I posted above - based on the OP's - to the Fft class: https://forum.bela.io/d/1659-a-small-ifft-addition-to-bela-fft-library/4
and pass as a second argument an array of zeros of the same size. Or you could write a method that only takes an array of real values and sets fdi(i) to 0.

lobo

giuliomoro Thanks for the response. I might have a wrong translation of the lines and I think I should have defined two vectors individually for real and imaginary.

for (int i = 0; i < gFftSize / 2 + 1; i++) {
	reIn[i] = gFft1.fdr(i) * gFft2.fdr(i) - gFft1.fdi(i) * gFft2.fdi(i);
	imIn[i] = gFft1.fdr(i) * gFft2.fdi(i) + gFft1.fdi(i) * gFft2.fdr(i);
}
for (int i = gFftSize / 2; i < gFftSize; i++) {
	reIn[i] = reIn[gFftSize - i];
	imIn[i] = -1 * imIn[gFftSize - i];
}

to translate the code

for (i = 0; i < N / 2 + 1; i++)
	fft_buf1[i] *= std::conj(fft_buf2[i]);

for (; i < N; i++)
	fft_buf1[i] = std::conj(fft_buf1[N - i]);

I am not familiar with c++; copy and paste both the header by OP and your function void Fft::ifft raises an error out-of-line definition of 'ifft' does not match any declaration in 'Fft' column. What am I missing?

giuliomoro

well, you'd need to modify the file in Bela/libraries/Fft/Fft.h to add void ifft(const std::vector<float>& reInput, const std::vector<float>& imInput); among the public methods in the class declaration and then modify Bela/libraries/Fft/Fft.cpp to include the code from my post. Editing these is not supported via the Bela IDE, so you'd have to ssh root@bela.local onto the board and navigate to the folder and use your favourite terminal editor. Is this something you could do?
Alternatively, copy the above Fft.* files in your current project by running the following in the console at the bottom of the IDE (replacing MYPROJECT with the name of your current project)

cp /root/Bela/libraries/Fft.* /root/Bela/projects/MYPROJECT

and edit them as above. Then change the #include in render.cpp to #include "Fft.h" so that the libraries are not used and only the local files are used.

lobo

giuliomoro Thanks that was clear! I am modifying Bela/libraries/Fft/Fft.cpp by simply adding the function as you suggested. But I am concerned that the myFft has to be something else defined in Fft.cpp? For Fft.h I replaced ifft() with the header line.

giuliomoro

I amended the code to remove the myFft reference, now it should build.

lobo Fft.h I replaced ifft() with the header line.

You shouldn't replace it, you should add it.

lobo

giuliomoro Yes it works! Thanks. I have a final dumb question: is td(x) equivalent to std:real(x) where x is the result of ifft()? I have extract the real coefficients of the results.

giuliomoro

The code assumes that the frequency domain response has conjugate symmetry and therefore it only sets the first half +1 bins in the frequency domain and only these are processed by the ne10 routine. As you know, a conjugate-symmetrical frequency domain representation corresponds to a real-only time domain signal. The result of the ifft ne10 routine is therefore real-only. So yes, td() gives you a real value, but there is never a complex time-domain value to begin with

giuliomoro

this ifft() with arguments is now on dev.