inverse FFT (DSP-168)

[stephan-ack]

Feature request:

Ability to make an inverse FFT with this library.

Br
Stephan

[dmitry1945]

Hi @StephanAAU,

the inverse FFT in our case is just multiply/divide by N, where N is number of samples.
You can make in with MulC function.

Regards,
Dmitry

[stephan-ack]

Okay, sounds very easily implementable.
I do not understand how the theory behind this works. Could you be so kind to point me towards some literature?

But also, just to make sure I understand correctly.
Example:

I have made an FFT on a signal and have it as an 2point array of the real and complex, with a size of N_SAMPLES.
Then to reconstruct the signal I will take the values at each index and multiply/devide by N_SAMPLES?

(Also, what is an MulC function?)

[stephan-ack]

It just occured to me that we might think of different things when saying "inverse fft".
I want to recontruct a signal to timedomain from frequency domain - I thought this was commonly known as inverseFFT, I might be wrong though.

[dmitry1945]

Hi @StephanAAU,

MulC - please look here:
https://github.com/espressif/esp-dsp/blob/master/modules/math/mulc/include/dsps_mulc.h

inverseFFT = N*FFT (or 1/N * FFT), just try.

But if you want to operate with signals in real-time (make FFT, filtering in frequency domain, then make inverse FFT) it will be more complicated.
You should look to this article "Turning Overlap-Save into a Multiband
Mixing, Downsampling Filter Bank" :
https://ieeexplore.ieee.org/document/1598092?arnumber=1598092&tag=1
I'm not sure is it open or not, but look.

Regards,
Dmitry

[stephan-ack]

Okay, I will give it a try.
Fortunately, what I want to achieve is not in realtime.

Thanks for the links and thanks a lot for your help.

Side note:
I think I got some understanding of how it works after looking at: https://www.embedded.com/dsp-tricks-computing-inverse-ffts-using-the-forward-fft/
It reminds me a bit of laplace transformations.

[stephan-ack]

Hi @dmitry1945

Sorry to bother you again.

I've tried to use the FFT example to create a sine, fft and then inverse fft with MulC, where C=N or C=1/N - but I do not get the same result. I suspect I might be trying to use the wrong signal, but I am abit at loss.
Do you have time to have a look at my code, maybe I am doing something obviously wrong.
Thanks in advance.

Heres my main code:

esp_err_t ret;
    ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);
    if (ret  != ESP_OK)
    {
        ESP_LOGE(TAG, "Not possible to initialize FFT. Error = %i", ret);
        return;
    }

    // Generate hann window
    dsps_wind_hann_f32(wind, N);

    // Generate input signal for x1 A=1 , F=0.1
    dsps_tone_gen_f32(x1, N, 0.9, 0.05,  0);

    // Convert into one complex vector (from two vectors, one with real and one with imag.)
    for (int i=0 ; i< N ; i++)
    {
        //y_cf[i*2 + 0] = y1_audio[i] * wind[i];
        //y_cf[i*2 + 1] = y2_audio[i] * wind[i];
        y_cf[i*2] = x1[i] * wind[i];
        y_cf[i*2 + 1] = 0;
    }
    
    // FFT calculated here
    unsigned int start_b = xthal_get_ccount();
    dsps_fft2r_fc32(y_cf, N);
    unsigned int end_b = xthal_get_ccount();

        
    // Bit reverse 
    dsps_bit_rev_fc32(y_cf, N);
    // Convert one complex vector to two complex vectors
    //dsps_cplx2reC_fc32(y_cf, N);
  
    // y1_cf - Result in log scale
    // y2_cf - Magnitude of your signal in linear scale
    for (int i = 0 ; i < N/2 ; i++) {
        y1_cf[i] = 10 * log10f((y1_cf[i * 2 + 0] * y1_cf[i * 2 + 0] + y1_cf[i * 2 + 1] * y1_cf[i * 2 + 1])/N);
        y2_cf[i] = ((y1_cf[i * 2 + 0] * y1_cf[i * 2 + 0] + y1_cf[i * 2 + 1] * y1_cf[i * 2 + 1])/N);
    }

    // dsps_mulc_f32_ansi(const float *input, float *output, int len, float C, int step_in, int step_out)
    /*
        input: input array
        output: output array
        len: amount of operations for arrays
        C: constant value
        step_in: step over input array (by default should be 1)
        step_out: step over output array (by default should be 1)
    */
    // 1/N = 0.000244140625
    dsps_mulc_f32(y_cf, x1_audio_recreated, N, 0.000244140625, 1, 1);

    // Show power spectrum in 64x10 window from -100 to 0 dB from 0..N/4 samples
    ESP_LOGW(TAG, "Signal 1 - logaritmic frequency domain");
    dsps_view(y1_cf,    N/2,    128,    10,     -20,    80,     '|');

    ESP_LOGW(TAG, "Signal 1 - linear frequency domain");
    dsps_view(y2_cf,    N/2,    128,    10,     0,    100,     '|');

    // ESP_LOGI(TAG, "FFT for %i complex points take %i cycles", N, end_b - start_b);

    ESP_LOGW(TAG, "Signal 1 - time domain");
    dsps_view(x1, 128,   128,    10,     -1,      1,    '.');

    //dsps_view(x1_audio_recreated, 256,   128,    10,     0,      254,    '.');
    ESP_LOGI(TAG, "Original: ");
    for (int i = 0; i < 128; i++) {
        ESP_LOGI(TAG, "%f", x1[i]);
    }

    ESP_LOGI(TAG, "\nRecreated: ");
    for (int i = 0; i < 128; i++) {
        ESP_LOGI(TAG, "%f", x1_audio_recreated[i]);
    }

    ESP_LOGI(TAG,"Done");

[dmitry1945]

Hi @StephanAAU,
I will look like this:

    // Convert into one complex vector (from two vectors, one with real and one with imag.)
    for (int i=0 ; i< N ; i++)
    {
        y_cf[i*2] = x1[i] * wind[i];
        y_cf[i*2 + 1] = 0;
    }
    // Here you have your complex signal. 

    // FFT calculated here
    dsps_fft2r_fc32(y_cf, N);
    // Bit reverse 
    dsps_bit_rev_fc32(y_cf, N);
    // Here you have your complex FFT spectrum

    // Here you can make your operations with complex spectrum at y_cf

    // FFT calculated here
    dsps_fft2r_fc32(y_cf, N);
    // Bit reverse 
    dsps_bit_rev_fc32(y_cf, N);
    // Here you have your original complex signal in time domain. 
    // y_cf[i*2 + 0] == x1[i]*wind(i)*N;
    // y_cf[i*2 + 1] == 0;
  
    // You can male multiplication with constant here:
    dsps_mulc_f32(y_cf, y_cf, 2*N, 0.000244140625, 1, 1);

[stephan-ack]

Thanks again @dmitry1945 !

[Aavu]

Hi @dmitry1945,
I was looking at computing ifft for real signals too. With the implementation you shared above, since you feed in the freq domain values without taking the complex conjugate, won't the output indices be reversed? Shouldn't we reverse the order of y_cf except the 1st sample?
For the order to not mess up, I think we should take the complex conjugate before performing the fft again (for ifft)? For complex signals we also need to take the complex conjugate of the output to get the original signal back.

So, inverseFFT = *(1/N)(*FFT), where * denotes complex conjugate operation.

Please correct me if I am missing something here.

[michikite]

Hi @dmitry1945, I was looking at computing ifft for real signals too. With the implementation you shared above, since you feed in the freq domain values without taking the complex conjugate, won't the output indices be reversed? Shouldn't we reverse the order of y_cf except the 1st sample? For the order to not mess up, I think we should take the complex conjugate before performing the fft again (for ifft)? For complex signals we also need to take the complex conjugate of the output to get the original signal back.

So, inverseFFT = *(1/N)(*FFT), where * denotes complex conjugate operation.

Please correct me if I am missing something here.

https://www.dsprelated.com/showarticle/800.php

[ayavilevich]

I was looking at computing ifft for real signals too.

@Aavu were you looking to do an ifft with a real output or an ifft with a real input? Did you manage to make an efficient code to calculate either of those? Thanks.

[dmitry1945]

Hi @ayavilevich,

You can use just normal FFT that exist, and change the sign of the output for imag part.

Regards,
Dmitry

[ayavilevich]

Hi @ayavilevich,

You can use just normal FFT that exist, and change the sign of the output for imag part.

Regards, Dmitry

@dmitry1945 yes, I understand that and I can do that after "mirroring" the results into N complex numbers and doing a full N fft. I was looking for a more optimal solution for the real case. I created another issue about it so not to hijack this thread.
Please take a look at #60 . There is more info there.

[coderforlife]

I have a working solution for doing the real IFFT. A slight correction needs to be made the Fourier-space data to make it inverse correctly. The below code implements forward and inverse FFTs using the ESP-DSP functions to do most of the work. It does require 2 additional passes through the data to perform complex conjugation. These passes could be avoided by making another version of dsps_cplx2real_fc32() and dsps_fft4r_fc32() with several of the signs flipped.

void rfft_init(int n_reals) {
    dsps_fft2r_deinit_fc32();
    dsps_fft4r_deinit_fc32();
    dsps_fft2r_init_fc32(NULL, n_reals >> 1);
    dsps_fft4r_init_fc32(NULL, n_reals >> 1);  // 4r is used for BOTH radix-2 and radix-4 FFTs!
}

void _rfft_core(float *x, int n) {
    if (n & 0xAAAAAAAA) {  // n is a power of 4
        dsps_fft4r_fc32(x, n);
        dsps_bit_rev4r_fc32(x, n);
    } else {
        dsps_fft2r_fc32(x, n);
        dsps_bit_rev2r_fc32(x, n);
    }
}

// complex conjugate of n complex numbers as pairs of floats (real, imag); operates in-place
void _conj(float *x, int n) { for (int i = 0; i < n; i++) { x[2*i+1] = -x[2*i+1]; } }

void rfft(float *x, int n) {
    // x has n real values or n/2 complex values, n must be a power of 2
    _rfft_core(x, n >> 1);
    dsps_cplx2real_fc32(x, n >> 1);
}

void irfft(float *x, int n) {
    // x has n complex values or n*2 real values, n must be a power of 2
    float re = x[0], im = x[1];

    _conj(x, n);               // conjugate the input for IFFT
    dsps_cplx2real_fc32(x, n); // reverse the post-processing (although not quite right)

    // Fix issues with `dsps_cplx2real_fc32()`
    x[0] = (re + im) * 0.5;  // fixes even values (and helps odd values) [DC component]
    x[1] = (im - re) * 0.5;  // fixes odd values [Nyquist component]

    _rfft_core(x, n);  // perform the FFT

    _conj(x, n); // fix the sign of every other real value

    // Scale the output by 1/n
    float scale = 1.0f / n;
    for (int i = 0; i < n * 2; i++) { x[i] *= scale; }
}