Как насчёт следующего алгоритма:
1) Берём FFT 128 точек, грубо находим максимум
2) В окрестности максимума запускаем DFT с шагом, как у FFT на 4096 точки, но только для маленького диапазона частот. Т.е. получаем кусок спектра из 50-60 точек с повышенной точностью.

Есть смысл? Как это будет соотноситься с поиском центра масс только по FFT-128?

Не уверен, но по моему так ничего не выйдет.
Попробуйте посмотреть на алгоритм герцеля.

P.S.: Не спешите, может еще кто чего дельного подскажет

1.А как вы его возьмёте? У вас обычно нет данных для этого... Добавить нулей? Можно, уже много раз обсуждалось...
Или Вы просто не хотите делать большой БПФ4096, а хотите поиметь только часть частот около максимума?? А данные добавлять всё равно нулевые...Тоже возможно, сделать на нужных частотах Герцеля, но не БПФ, поскольку это одно и то же по результату. В любом случае потом по точкам вблизи максимума подгоняют параболу (или даже функцию спектрального окна, но это не для реального времени)


2.Вот Вам простой и самый точный из известных мне способ интерполяции вблизи максимума, если без дополнительных преобразований. Значительно точнее чем по центру масс. Но как и в случае с добавлением нулей будет содержать небольшие систематические ошибки. Т.е. при очень высоком SNR будет работать теоретически не оптимально

Y is a three-element complex vector with the
% DFT output magnitude maximizer as the center element Y(kmax-1), Y(kmax), Y(kmax+1)).
%
% Returns -0.5 < x < 0.5, which is the fraction of the sample
% spacing (i.e., bin width) about the center element where the
% peak is estimated to be.

% Ref: Macleod, M.D., "Fast Nearly ML Estimation of the Parameters
% of Real or Complex Single Tones or Resolved Multiple Tones,"
% IEEE Trans. Sig. Proc. Vol 46 No 1,
% January 1998, pp141-148.

ref = Y(2); % Isolate phase reference.
R = real(Y.*conj(ref)); % Generate phase corrected coefficient vector.
gamma = (R(1)-R(3))/((2*R(2))+R(1)+R(3)); % Calculate offset.
delta = (sqrt(1 + 8*gamma*gamma)-1)/(4*gamma); % Final estimate.
x = delta;

Можете найти полжение максимума по этой формуле, потом прогнать Герцель в окрестности. Герцель это фильтр, который в точности даёт результат БПФ. но только в одной частотной точке. Можно найти в первом приближении центральную частоту по Macleod, прогнать Герцель на этой частоте f0 и двух соседних f0-0.25df, f0+0.25df, df шаг БПФ на 128. Потом провести параболу и посмотреть на какой частоте её вершина. Называют это авторы ML-extension

Да, я хотел делать именно так:
БПФ128 находит максимум X,
потом идёт цикл ДПФ4096, но i пробегает значения не 0..4095, а (X*32-10)..(x*32+10)
В обоих случаях ничего нулей не добавляю.

пс Вот смотрю приведённый вами интересный метод.

Это не ДПФ. По определению ДПФ содержит только N ортогональных частот, где N - это количество исходных отсчетов
Формально Вы можете вычислить такие же суммы для более мелкой равномерной сетки, но это будет математически всё равно что добавить нулей до 4096 точек и взять ДПФ на 4096. Конечно, если Вы считаете в лоб ДПФ (или Герцеля), а не по БПФ, при вычислениях нулей не добавляете, зачем, их вклад в сумму нулевой. Результат во всех случаях будет один

С "грубо берём максимум" есть ещё одна проблема, если потом захочется брать очень тонко. Можно иногда промахнуться.
Кроме обычных ошибок, вызванных шумом, возможны ещё и ошибки другой природы при частоте на границе бинов. А именно максимум ложный. Если значения мощности в двух бинах примерно равны, то шум может перекинуть максимум к ложному бину. После этого мелкая сетка Вам уже не поможет... Поэтому при последовательном уточнении нужно быть аккуратным с максимумом - это он максимум, когда есть явно доминирующее значение, а если они примерно равны лучше взять вместо максимума что-то более похожее на центр масс

По точности http://dsplib.ru/forum/viewtopic.php?f=7&t=142
По методам http://forum.vingrad.ru/forum/tech-algorit...que-method.html

Можно.
Но БПФ128 надо повторить 32 раза для 32 входных массивов, затем умножить на вращающие коэффициенты и сделать
ДПФ32 с бинами, которые интересуют.
Будет укороченный БПФ4096 по основанию 128 и 32 с прореживанием по частоте и с 128;32 инверсией на выходе.

Так вот эти N ортогональных частот можно не все считать, а только одну, пять, шесть и т.п.. Это как раз-таки ДПФ, потому как

но это не по "горизонтали" так сказать, а по "вертикали".
То есть I - это не номер отсчёта, а частота. Отстчёты же идут по J от 0 до 4095.
Соответственно никаких нулей никто никуда не добавляет.

Может я неправильно донёс свою мысль. Попробую по-другому.
Беру БПФ - получаю спектр
Беру ДПФ - получаю спектр, приём тютелька в тютельку амплитуды все те же самые. Идея-то в чём, у ДПФ амплитуда i-той частоты amp[i] получается как сумма перемножений j-тых отсчётов на соответствующие коэффициенты. Можно и i=12345-той частоты посчитать, а остальные просто не трогать даже. Это в БПФ фиг что выбросишь, всё уже завязано в клубок.

Финальная идея такова:
1. БПФ на 128 используя отсчёты I=0..127 , интерполяция спектра до 2048 точек и поиск максимума N.
2. ДПФ по частотам J=N-10...N+10, используя отсчёты I=0..4095
3. После захвата максимума делаем уже только скользящее ДПФ, при потере - переход на пункт 1.

Скользящее ДПФ: для частоты J вычитаем из суммы произведений самое "старое" произведение и прибавляем "новый" отсчёт умноженный на соответственный коэффициент. Для 20 частотных бинов 20*2 вычитаний + 20*2 умножений + 20*2 сложений на каждый пришедший отсчёт с АЦП (*2, т.к. спектр комплексный).