Python rfft algorithm

Question

Numpy has a rfft function to compute FFT or Real input. I would like to code my own code for that. Here's my actual code of DFT and FFT [ref].


def DFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array x"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

def FFT(x):
    """A recursive implementation of the 1D Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]

    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    elif N <= 32:  # this cutoff should be optimized
        return DFT_slow(x)
    else:
        X_even = FFT(x[::2])
        X_odd = FFT(x[1::2])
        factor = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + factor[:N / 2] * X_odd,
                               X_even + factor[N / 2:] * X_odd])

I don't know much about Fourier transform, any idea how I would compute the rfft from these codes?

score 0 · Accepted Answer · answered Feb 03 '20 at 23:49

0

The rfft is just the positive-frequency components of the full fft of the input, since the fft of a real-valued signal is conjugate-symmetric in frequency.

def rDFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array of real data x"""
    x = np.real(np.asarray(x, dtype=float))
    N = x.shape[0]
    n = np.arange(N)
    k = np.arange(N//2+1).reshape((N//2+1,1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

There are some trickier factors that go into implementing the recursive version when you're only calculating half of the coefficients, but the slow rDFT should get you started.

answered Feb 03 '20 at 23:49

Mark Snyder

1,635
3
12

How can I compute the inverse of such a transform? I tried to use a pseudo inverse but the result is not very good... – Feb 04 '20 at 18:12
1

@Dave If you want to do it without using `np.fft.irfft`, you can reconstruct the signal using `N//2+1` sinusoids with amplitudes and phases proportional to the transform components. You can see my answer [here](https://stackoverflow.com/a/59726152/12482432) for more details - the second implementation only uses `rfft` coefficients. – Mark Snyder Feb 04 '20 at 18:36
Thank you, very helpfull. Thanks to your code I made a good `ifft` function. But I have trouble to compute the `irfft`. I tried [this](https://pastebin.com/raw/ZYDP7R9e) but it returns a sine twice as big as the original. Can you help me? – Feb 05 '20 at 13:36
@Dave When you calculate the individual sinusoids, you should divide by `len(x)` (the length of the reconstructed signal), not `N` (the length of the `rfft`). – Mark Snyder Feb 05 '20 at 17:53
1

@Dave Also, you should be careful: if the original `len(x)` is odd, `2*(N-1) != len(x)`. Instead, `2*(N-1) == len(x)-1`. You can fix this by checking if `X[-1].real == X[-1]`. If it does, the original `len(x)` was `2*(N-1)` because `len(x)` was even. If it doesn't, the original `len(x)` was `2*N-1` because `len(x)` was odd. – Mark Snyder Feb 05 '20 at 18:01

Python rfft algorithm

1 Answers1