[tmtlib 0]

Как насчёт следующего алгоритма:

1) Берём FFT 128 точек, грубо находим максимум

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

 

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

[eugen_pcad_ru 0]

Не уверен, но по моему так ничего не выйдет.

Попробуйте посмотреть на алгоритм герцеля.

 

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

[fontp 0]

Как насчёт следующего алгоритма:

1) Берём FFT 128 точек, грубо находим максимум

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

 

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

 

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

[tmtlib 0]

1.А как вы его возьмёте? У вас обычно нет данных для этого... Добавить нулей? Можно, уже много раз обсуждалось...

Или Вы просто не хотите делать большой БПФ4096, а хотите поиметь только часть частот около максимума??

Да, я хотел делать именно так:

БПФ128 находит максимум X,

потом идёт цикл ДПФ4096, но i пробегает значения не 0..4095, а (X*32-10)..(x*32+10)

В обоих случаях ничего нулей не добавляю.

 

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

[fontp 0]

Да, я хотел делать именно так:

БПФ128 находит максимум X,

потом идёт цикл ДПФ4096, но i пробегает значения не 0..4095, а (X*32-10)..(x*32+10)

В обоих случаях ничего нулей не добавляю.

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

 

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

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

 

С "грубо берём максимум" есть ещё одна проблема, если потом захочется брать очень тонко. Можно иногда промахнуться.

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

[Santik 0]

По точности http://dsplib.ru/forum/viewtopic.php?f=7&t=142

По методам http://forum.vingrad.ru/forum/tech-algorit...que-method.html

 

[aser 0]

Можно.

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

ДПФ32 с бинами, которые интересуют.

Будет укороченный БПФ4096 по основанию 128 и 32 с прореживанием по частоте и с 128;32 инверсией на выходе.

[tmtlib 0]

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

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

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

потом идёт цикл ДПФ4096, но i пробегает значения не 0..4095, а (X*32-10)..(x*32+10)

но это не по "горизонтали" так сказать, а по "вертикали".

То есть I - это не номер отсчёта, а частота. Отстчёты же идут по J от 0 до 4095.

Соответственно никаких нулей никто никуда не добавляет.

 

Может я неправильно донёс свою мысль. Попробую по-другому.

Беру БПФ - получаю спектр

Беру ДПФ - получаю спектр, приём тютелька в тютельку амплитуды все те же самые. Идея-то в чём, у ДПФ амплитуда i-той частоты amp получается как сумма перемножений 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, т.к. спектр комплексный).