Как рассчитать БПФ большого размера, используя БПФ меньшего размера?

Вот алгоритм вычисления БПФ размера P с использованием двух меньших функций БПФ размеров M и N (в исходном вопросе называются размеры M и k).

Входные данные:
P — размер большого БПФ, который вы хотите вычислить.
M, N выбираются таким образом, что MN=P.
x[0...P-1] — входные данные.

Настройка:
U — двумерный массив с M строк и N столбцов.
y — это вектор длины P, который будет содержать БПФ x.

Алгоритм:
шаг 1. Заполнить U от x по столбцам, чтобы U выглядело так:
x(0) x(M) ... x(P-M)
x(1) x(M+1) ... x(P-M+1)
x(2) x(M+2) ... x(P-M+2)
... ... ... ...
x(M-1) x(2M-1) ... x(P-1)

шаг 2. Замените каждую строку U собственным БПФ (длины N).
шаг 3. Умножьте каждый элемент U(m,n) на exp(-2*pi*j*m* н/п).
шаг 4. Замените каждый столбец U своим собственным БПФ (длиной M).
шаг 5. Считайте элементы U по строкам в y, например:

y(0) y(1) ... y(N-1)
y(N) y(N+1) ... y(2N-1)
y(2N) y(2N+1) ... y(3N-1)
... ... ... ...
y(P-N) y(P-N-1) ... y(P-1)

Здесь это код MATLAB, который реализует этот алгоритм. Вы можете проверить это, набрав fft_decomposition(randn(256,1), 8);

function y = fft_decomposition(x, M)
% y = fft_decomposition(x, M)
% Computes FFT by decomposing into smaller FFTs.
%
% Inputs:
% x is a 1D array of the input data.
% M is the size of one of the FFTs to use.
%
% Outputs:
% y is the FFT of x.  It has been computed using FFTs of size M and
% length(x)/M.
%
% Note that this implementation doesn't explicitly use the 2D array U; it
% works on samples of x in-place.

q = 1;   % Offset because MATLAB starts at one.  Set to 0 for C code.
x_original = x;
P = length(x);
if mod(P,M)~=0, error('Invalid block size.'); end;
N = P/M;

% step 2: FFT-N on rows of U.
for m = 0 : M-1
    x(q+(m:M:P-1)) = fft(x(q+(m:M:P-1)));
end;

% step 3: Twiddle factors.
for m = 0 : M-1
    for n = 0 : N-1
        x(m+n*M+q) = x(m+n*M+q) * exp(-2*pi*j*m*n/P);
    end;
end;

% step 4:  FFT-M on columns of U.
for n = 0 : N-1
    x(q+n*M+(0:M-1)) = fft(x(q+n*M+(0:M-1)));
end;

% step 5:  Re-arrange samples for output.
y = zeros(size(x));
for m = 0 : M-1
    for n = 0 : N-1
        y(m*N+n+q) = x(m+n*M+q);
    end;
end;

err = max(abs(y-fft(x_original)));
fprintf( 1, 'The largest error amplitude is %g\n', err);
return;
% End of fft_decomposition().

Вы можете просто использовать последние log2 (k) проходы БПФ по основанию 2, предполагая, что предыдущие результаты БПФ получены из подмножеств данных с соответствующим чередованием.

Ну, БПФ - это в основном рекурсивный тип преобразования Фурье. Это основано на том факте, что, как говорит Википедия:

Самые известные алгоритмы БПФ зависят от факторизации N, но существуют БПФ со сложностью O(N log N) для >всех N, даже для простых N. Многие алгоритмы БПФ зависят только от того факта, что e^(-2pi* i/N) является N-м первообразным корнем из единицы и, таким образом, может применяться к аналогичным преобразованиям над любым конечным полем, таким как теоретико-числовые преобразования. Поскольку обратное ДПФ — это то же самое, что и ДПФ, но с обратным знаком экспоненты и коэффициентом 1/N, любой алгоритм БПФ можно легко адаптировать для него.

Так что это в значительной степени уже было сделано в БПФ. Если вы говорите о получении сигналов с более длительным периодом из вашего преобразования, вам лучше выполнить ДПФ для наборов данных с ограниченными частотами. Возможно, есть способ сделать это из частотной области, но IDK, если кто-то действительно это сделал. Ты мог быть первым!!!! :)

Ответ kevin_o сработал довольно хорошо. Я взял его код и устранил циклы, используя некоторые базовые приемы Matlab. Функционально идентичен его версии

function y = fft_decomposition(x, M)
% y = fft_decomposition(x, M)
% Computes FFT by decomposing into smaller FFTs.
%
% Inputs:
% x is a 1D array of the input data.
% M is the size of one of the FFTs to use.
%
% Outputs:
% y is the FFT of x.  It has been computed using FFTs of size M and
% length(x)/M.
%
% Note that this implementation doesn't explicitly use the 2D array U; it
% works on samples of x in-place.

q = 1;   % Offset because MATLAB starts at one.  Set to 0 for C code.
x_original = x;
P = length(x);
if mod(P,M)~=0, error('Invalid block size.'); end;
N = P/M;

% step 2: FFT-N on rows of U.
X=fft(reshape(x,M,N),[],2);

% step 3: Twiddle factors.
X=X.*exp(-j*2*pi*(0:M-1)'*(0:N-1)/P);

% step 4:  FFT-M on columns of U.
X=fft(X);

% step 5:  Re-arrange samples for output.
x_twiddle=bsxfun(@plus,M*(0:N-1)',(0:M-1))+q;
y=X(x_twiddle(:));

% err = max(abs(y-fft(x_original)));
% fprintf( 1, 'The largest error amplitude is %g\n', err);
return;
% End of fft_decomposition()