What is numpy.fft.rfft and numpy.fft.irfft and its equivalent code in MATLAB

I am converting a python code into MATLAB and one of the code uses numpy rfft. In the documentation of numpy, it says real input.

Compute the one-dimensional discrete Fourier Transform for real input.

So what I did in MATLAB is using abs but the results are different.

Python code

ffta = np.fft.rfft(a) 


ffta = abs(fft(a));

What have I misunderstood?


  • The real FFT in numpy uses the fact that the fourier transform of a real valued function is so to say "skew-symmetric", that is the value at frequency k is the complex conjugate of the value at frequency N-k for k=1..N-1 (the correct term is Hermitian). Therefore rfft returns only the part of the result that corresponds to nonpositive frequences.

    For an input of size N the rfft function returns the part of the FFT output corresponding to frequences at or below N/2. Therefore the output of rfft is of size N/2+1 if N is even (all frequences from 0 to N/2), or (N+1)/2 if N is odd (all frequences from 0 to (N-1)/2). Observe that the function floor(n/2+1) returns the correct output size for both even and odd input sizes.

    So to reproduce rfft in matlab

    function rfft = rfft(a)
         ffta = fft(a);
         rfft = ffta(1:(floor(length(ffta)/2)+1));

    For example

    a = [1,1,1,1,-1,-1,-1,-1];
    rffta = rfft(a)

    would produce

    rffta =
     Columns 1 through 3:
       0.00000 + 0.00000i   2.00000 - 4.82843i   0.00000 + 0.00000i   
     Columns 4 through 5:
       2.00000 - 0.82843i   0.00000 + 0.00000i

    Now compare that with python

    >>> np.fft.rfft(a)
    array([ 0.+0.j        ,  2.-4.82842712j,  0.-0.j        ,  
            2.-0.82842712j,  0.+0.j        ])

    Reproducing irfft

    To reproduce basic functionality of irfft you need to recover the missing frequences from rfft output. If the desired output length is even, the output length can be computed from the input length as 2 (m - 1). Otherwise it should be 2 (m - 1) + 1.

    The following code would work.

    function out = irfft(x, even)
        if nargin == 1
            even = true;
        % `n`: the output length
        % `s`: the variable that will hold the index of the highest
        %      frequency below N/2, s = floor((n+1)/2)
        N = numel(x);  % function assumes 1D input
        if even
            n = 2 * (N - 1);
            s = N - 1;
            n = 2 * (N - 1) + 1;
            s = N;
        outf = zeros(1, n);
        outf(1:N) = x;
        outf(N+1:n) = conj(x(s:-1:2));
        out = ifft(outf);

    Now you should have

    >> irfft(rfft(a))
    ans =
       1.00000   1.00000   1.00000   1.00000  -1.00000  -1.00000  -1.00000  -1.00000

    and also

    abs( irfft(rfft(a)) - a ) < 1e-15

    For odd output length you get

    >> irfft(rfft(a(1:7)),even=false)
    ans =
       1.0000   1.0000   1.0000   1.0000  -1.0000  -1.0000  -1.0000