0
$\begingroup$

I am working on a fast autocorrelation function using the FFT (similar to scipy.signal.fftconvolve). The algorithm is pretty straight forward but I feel like I'm leaving an optimization on the table and I can't figure out how to make it work. Here is an outline of the (usual?) algorithm:

  1. Run a real signal $x$ through the discrete Fourier transform: $X = \mathfrak{F}\{x\}$.
  2. Replace the entries in $X$ with their norm: $Y = X \cdot X^*$. (Don't care about DC for now.)
  3. Run $Y$ through the inverse discrete Fourier transform the get the autocorrelated real signal: $y = \mathfrak{F}^{-1}\{Y\}$

With the forward FFT I can make use of the fact that I have a real signal, so the computation is faster than for complex signals of the same length. So far so good. Here is where I'm stuck: $Y$ is also a real signal, since the product of a complex number with its conjugate is always real. Right now I am using a backward real FFT function that expects a complex signal and returns a real signal. Is there a trick similar to the forward real FFT? Something like a Decimation In Time (DIT) algorithm for the input of the IFFT.

Here is an example in Python:

import numpy as np
x = [2, 1, 3, 1, 3, 1, 2, 1]
X = np.fft.rfft(x)
Y = X * np.conj(X)
y = np.fft.irfft(Y)  # <- feel like this could be improved, Y is a real signal.
print(x, X, Y, y, sep='\n')

# x: [2 1 3 1 3 1 2 1]
# X: [14.+0.j -1.-1.j  0.-0.j -1.+1.j  6.+0.j]
# Y: [196.+0.j   2.+0.j   0.+0.j   2.+0.j  36.+0.j]
# y: [30. 20. 29. 20. 28. 20. 29. 20.]
$\endgroup$

1 Answer 1

0
$\begingroup$

There is such a function, though numpy might not expose it: it's still a real-valued-input FFT! So, it's, aside from an exponent's sign, the same as your rfft. Matter of fact, FFT libraries like FFTw do expose a backwards r2c FFT variant.

@Richard Lyons has a nice article on how to compute the inverse FFT just with the help of forward FFTs and basic transformations. What you'll like is this method #4:

Method #4, description below

You need to

  1. take the complex conjugate of your input, then
  2. calculate the forward FFT of the result, then
  3. complex conjugate of the output.

And fun fact: your input to the rfft is real, so step 1. is a NOP (no-operation). Step 2 is your rfft, and step 3 is still relatively trivial to compute.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.