persistent_fft.m

来自「一个关于FFT的例子」· M 代码 · 共 42 行

function f=persistent_fft(x)
% F = persistent_fft( X )
%   Calculates the fft of a column vector x
%   results are, theoretically, same as matlab's fft
%   but you can expect discrepancies ~ 1e8 in normalized L_2 error.
%
%   Performs one step of the radix-2 algorithm, then calls matlab's fft
%
%   Vectorized, and uses persistent variables to save time
%   ... but alas, still slower than matlab's fft
%
%   For the fft of a fixed length, this code runs several orders
%   of magnitude faster on *subsequent* calls, due to the persistent
%   variables (note: matlab's fft also runs faster on subsequent calls, due
%   to learned "wisdom" in the FFTW algorithm)
%
%   Why you want to use this code over matlab's own fft()?  I do not have
%   an answer to that.
%
%   Stephen Becker, 2/22/2008  [e-mail]

persistent weights N

n = length(x);
if size(x,2) ~= 1 || n < 4 || round( n/2 ) ~= n/2
    error('Input must be a column vector of even length >= 4')
end

% Decide whether to recalculate the weights or not:
if isempty(N) || n ~= N 
    N = n;
    w = exp(-2*pi*i / N);
    helper = w .^ [1:(N/4-1)].' ;
    weights = [1;helper; -i ;flipud( -real(helper)+i*imag(helper) ) ];
%     weights = w .^ [0:(N/2-1)].' ;    % same thing, but slower
end

% A simple radix-2 DIT:
E = fft( x(1:2:N) );
O = fft( x(2:2:N) );
O = O .* weights;
f = [ E + O ; E - O ];