Fourier Transform for real input optimisation

We have real signal with N points. esp-fftr.zip

#include "dsps_fft2r.h"

void main()
{
    enum {N=8};
    float x1[N] = {0,1,2,3,4,5,6,7};
    float y_cf[N*2];

    dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);

    for (int i = 0; i < N; i++)
    {
        y_cf[i*2 + 0] = x1[i];
        y_cf[i*2 + 1] = 0.0; // only real part
    }
    // FFT
    dsps_fft2r_fc32_ansi(y_cf, N);
    // Bit reverse
    dsps_bit_rev_fc32_ansi(y_cf, N);

    for (int i = 0 ; i < N ; i++) {
        printf("%f,%f\n", y_cf[i * 2 + 0], y_cf[i * 2 + 1]);
    }
}

Output:

28.000000,0.000000
-4.000000,9.656854
-4.000000,4.000000
-4.000000,1.656854
-4.000000,0.000000
-4.000000,-1.656854
-4.000000,-4.000000
-4.000000,-9.656855

Result is same as numpy.fft.fft(range(8)). As we can see output is symmetric. There is a special optimized function in numpy for this case called rfft. From numpy doc:

When the DFT is computed for purely real input, the output is Hermitian-symmetric, i.e. the negative frequency terms are just the complex conjugates of the corresponding positive-frequency terms, and the negative-frequency terms are therefore redundant. This function does not compute the negative frequency terms, and the length of the transformed axis of the output is therefore n//2 + 1.

It it would be great to have such rfft in esp dsp.

Maybe it's easier to show on python. This is straightforward dft:

def compute_dft_complex(input):
  n = len(input)
  output = []
  for k in range(n):  # For each output element
    s = complex(0)
    for t in range(n):  # For each input element
      angle = 2j * cmath.pi * t * k / n
      s += input[t] * cmath.exp(-angle)
    output.append(s)
  return output

This is optimised rdft:

def compute_rdft_complex(input):
  n = len(input)
  output = []
  for k in range(n//2+1):  # !!! almost 2x time faster !!!
    s = complex(0)
    for t in range(n):  # For each input element
      angle = 2j * cmath.pi * t * k / n
      s += input[t] * cmath.exp(-angle)
    output.append(s)
  return output

Test:

import cmath
import numpy as np
np.allclose(np.fft.fft(x), compute_dft_complex(x))
np.allclose(np.fft.rfft(x), compute_rdft_complex(x))

espressif / esp-dsp

Fourier Transform for real input optimisation #22