This short bit of code is a full fledged radix-2 fast fourier transform.
It is kept short and easy to read to make it easy for students to learn and understand the fast fourier algorithm.
As this algorithm is recursive, it does not require the butterfly reordering of input or output data.
Length of the input must be a power of 2.
y = mfft(1:2^10)
n = 2^10;
data = sind((1:n)*10)'+rand(n,1);
y = mfft(data);
Friedrich (2019). mfft (https://www.mathworks.com/matlabcentral/fileexchange/67515-mfft), MATLAB Central File Exchange. Retrieved .
Yep. Works now. I like it. Didn't know that is is possible to make it that short.
Thank you for your findings. I actually introduced the bug by my last update, it worked before, but you are right about the explicit definition of the output.
The bug slipped through my test as the old version must have somehow been cached. It will work now.
PS: old code for reference (without explicit definition of the output)
function x = mfft(x)
% X = MFFT(X) calculates the fast discrete fourier transform of vector x.
% Length of x must be a power of two. This implementation is
% intended to help understand the fast fourier algorithm. It is not
% optimized for speed.
% Copyright 2018
% Friedrich Welck, 28 May 2018
len = length(x);
if len >= 2
% calculate fft of odd elements
odd = mfft(x(1:2:len));
% calculate fft of even elements
even = mfft(x(2:2:len));
% rotate even elements to prepare for recombination
even = exp( (0:len/2-1)*(-2i*pi/len) ) .* even;
% recombine fft of odd and even elements
x = [odd + even , odd - even ];
A simple implementation of the FFT algorithm is very welcome, as this allows to study how the recursive computation works. Unfortunately, this does not run in R2017a. The example above results - as expected - in the error:
Output argument "y" (and maybe others) not assigned during call to "mfft"
Even if a different release of Matlab accepts this code, it should be changed, because the intended audience - students - should learn that the case len==1 is handled explicitly by an assignment to y and not by any black-box magic one release might accept. This would also help to make the code portable. It is very desirable to have code "Compatible with any release" if no special features are necessary to solve the problem.
Will be very useful once the bug is corrected.
bugfix as found by Heinrich Acker