pythonmatlabnumpyfft

# 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)
``````

MATLAB code

``````ffta = abs(fft(a));
``````

What have I misunderstood?

Solution

• 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));
end
``````

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;
end
% `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;
else
n = 2 * (N - 1) + 1;
s = N;
end
outf = zeros(1, n);
outf(1:N) = x;
outf(N+1:n) = conj(x(s:-1:2));
out = ifft(outf);
end
``````

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
``````