FFTWN: Round up to efficient FFT integer

FFTWN Round up to appropriate integer for efficient FFT computation

It is well known that the FFT algorithm is more efficient for some
data lengths. Often, it is recommended that users zero-pad data
to the next power of 2 to maximize FFT efficiency and minimize
computation time. This strategy may significantly increase the
data length and actually increase the FFT calculation time due
to the burden of calculating the extra points (See example 2
below).

The MATLAB FFT implementation uses the FFTW library. With default
settings, FFTW is highly efficient for data lengths of the form
N = 2^a*3^b*5^c*7^d*11^e*13^f, where e + f = 0 or 1, and the rest
of the exponents are arbitrary.

This function accepts a non-negative integer as its argument, and
outputs an integer greater or equal to the input that conforms
to the above equation. The algorithm simply increments the value
and re-checks until it arrives at an acceptable integer.

FFTW Reference:
http://www.fftw.org/fftw3_doc/Real_002ddata-DFTs.html#Real_002ddata-DFTs

Usage:
m = fftwn(n)

n: non-negative integer to be checked
m: next highest FFTW integer >= n

Example 1: Typical usage example
N = 1546; T = 100; % no. of data points, total time
dt = T/N % time step
t = dt*(0:N-1)'; % time vector
x = 3*sin(2*pi*3.2*t) + 0.05*randn(N,1); % signal + noise
NFFT = fftwn(N) % NFFT = 1560, optimum fft length
X = fft(x,NFFT)/N; % compute normalized DFT
X_pos = 2*abs(X(1:NFFT/2+2)); % single-sided magnitude
df = 1/(dt*NFFT); % frequency step
f = df*(0:NFFT/2+1)'; % frequency vector
stem(f,X_pos); % stem plot to check result

Example 2: Demonstrate potential time advantage
timeit is available from the file exchange, File ID: #18798
Fs = 1024; % arbitrary
N = 17^6; % N = 24,137,569; chosen to be tough for FFTW
Nfftwn = fftwn(N) % Nfftwn = 24,147,200; increase of 9,631
Nnp2 = 2^nextpow2(N) % Nnp2 = 33,554,432; increase of 9,416,863
t = 1/Fs*(0:N-1)';
x = 3*sin(2*pi*14*t) + 0.05*randn(N,1);
f = @() fft(x,N);
g = @() fft(x,Nfftwn);
h = @() fft(x,Nnp2);
timeit(f) % ans = 6.43; time for unmodified N
timeit(g) % ans = 6.05; FFTWN(N) faster in this case
timeit(h) % ans = 7.40; 2^np2(N) slower in this case

See also: FFT, FFTW