Поможет БИХ фильтр. Если есть вычислительные ресурсы, то хороший фильтр избавит от гармоник. Единственно, что плохо так это вычисление в несколько этапов.

Фильтр не поможет, т.к. проблема не в гармониках - они и так хорошо разделяются Фурье-преобразованием, а в том, что частота плывет. И в зависимости от того, в какую часть семпла мы смотрим лучше (где наложено окно), тот ответ мы и получаем.
Теперь про результаты анализа файла с 32к оцифровки. Тяжело плавать в серной кислоте с оторванными ногами. У меня матлаб не стоит, так что пришлось потрошить файл вручную. Где-то я ошибся на один семпл и у меня прошел сбой по частоте. Но в целом этот файл гораздо лучше последнего, частота болтается не больше, чем 2 единицы второго знака после запятой, хотя амплитуда маленькая и шумов довольно много. Если сохраните файл в текстовом виде, могу посчитать еще раз поточнее. Первая прикидка частоты - 822.557.

Да конечно пожалуйста. Сконвертировал из *.daq и добавил в архив и текстовый файл, расширение "csv". Нажмите для просмотра прикрепленного файла

И опять засада: у меня 822.5698 Hz

четверти от выборки (по 1024 точек):
1: 822.5570
2: 822.5868
3: 822.5385
4: 822.5531
половинки от выборки (по 2048 точек):
1: 822.5842
2: 822.5649
Полная выборка (4096 точек): 822.5698

Получается, что измеряемая частота сильно плывет? такое теоретически возможно. Но как-то не вяжется с тем что я вижу глазами: я последовательно применяю два метода (подсчет периодов методом пересечение нуля и FFT-based: первый-второй-первый-второй...), и вычисленная каждым из них величина колеблется около своего значения частоты, причем FFT-baised больше. И разность того же порядка как разница между Вашей и моей цифрами.

Вот уточненные результаты: средняя частота 822.5666. Если смотреть по частям (я, правда, делал несколько иначе - оставлял FFT на 4096, но делал узкое окно и сдвигал его относительно данных), то на 1/4 от начала частота 822.5897, на 1/2 - 822.5614, на 3/4 - 822.5633. При делении на 8 получается следующее:
822.5477
822.6066
822.5791
822.5578
822.5175
822.5213
822.5499
Видно, что чем мельче интервал анализа, тем больше разброс частоты. При большем усреднении разброс, естественно, уменьшается. Для данного файла болтанка частоты 9 единиц во втором знаке после запятой.
Я посмотрел на модели - синус плюс белый шум примерно равный тому, что есть в этом файле. Ошибка частоты - примерно 0.03 Гц максимум при делении на 8. Таким образом, здесь кроме шума действительно есть небольшое плавание частоты (примерно на те же 0.03 Гц). На полном интервале данный шум приводит к ошибке меньше 0.01 Гц. Отсюда следуют рекомендации (очевидные) - увеличивайте отношение сигнал/шум и боритесь с плаванием частоты.

Понял, большое спасибо!
Насчет плавания частоты- подозреваю, что оно все-таки имеет чисто физический смысл. Датчик деформации хорошо вибрации берет, и это может быть скажем собственная частота колебаний здания (ну, наверное может мой 7-й этаж так болтаться). То есть такая точность измерений бесполезна с тем типом датчика, что у меня сейчас есть, так он и фазы луны и Камаз в километре от точки измерения будет регистрировать.
В-общем, попрошу померить в плевых условиях/бункерах, а данные оно мне будет на флэшку кидать- соберу статистику по плаванию частоты и искажениям с разными датчиками, может чего умное в голову приползет.........

Ну и опять же про совершенствование моего метода обработки буду думать- так и не врубился, почему у меня значение больше чем у Вас на сотые доли герца. А Вашим результатам я доверяю больше чем своим

Так Вы проверьте свою программу на модели. Сгенерируйте синус программно в файл и напустите на него Вашу программу. Сразу станет понятно, кто врет, т.к. Вы будете знать заранее, что должно получиться. Потом добавьте шума и посмотрите как Ваша программа на него реагирует.

на чистом синусе все шикарно, внутри бина ошибка (разность вычисленной и теоретической частот) соответствует предсказанной для гауссовского окна, именно так я и подбирал коэффициенты для окна (чтобы не более 0.001 Гц максимальная ошибка была).
Более того- и в приборе при подаче чистого синуса все отлично, ну прямо почти единицы наносекунд совпадают с частотомером.
А в реальности что-то не очень. Одни только очень серьезные палки на нечетных гармониках чего стоят. А фильтрануть не могу- полоса прибора позволяет куче гармоник пролазить. Реально хорошо вижу все нечетные гармоники, сколько аппаратный ФНЧ оставляет.

С шумами не смотрел на модели вообще, есть такой грех. Смотрел только разрешающую способность и точность при чистой синусоиде и радовался как слон Посмотрю... Кстати и с гармониками посмотрю....

Есть еще идея, что возбуждение свип-меандром создает массу несобственных колебаний, которые затухают не так быстро как я думал. При низком разрешении/точности это не влияет на результат, а в моем случае вылезает в полный рост. Но это все уже физика процесса, а не обработка данных..... Опять же будем глядеть (да хоть просто попробую от синус-свип-генератора возбудить).
Кстати я приводил картинку, которую производители суперкрутого прибора привели:вот. Ну нет у них гармоник, и сигнал затухает в сотни раз быстрее чем у меня. Может еще и спецдатчики нужны для такого прибора.... Или демпфирование.....

Но это все так, рюшечки для себя любимого и набирание опыта с Матлабом. Поставленная при старте топика задача решена

Еще раз спасибо!

С гармониками бороться довольно просто - нужно, чтобы окно было выбрано так, чтобы на половине расстояния от основного тона до гармоники (по частоте) спектр спадал достаточно сильно. Если взять отношение е-7 - е-8, то достаточно. Хуже с шумом - он-то как раз в полосу попадает и всем мешает.
Производители суперкрутого прибора врут как могут. При соотношении сигнал/шум 21, как у них нарисовано, никакой точности в 3 знака после запятой получить нельзя. Какой там шум на самом деле - не видно, т.к. шкала должна быть логарифмическая.

Попробовал демпфировать датчик резистором- получил увеличение отношения сигнал/гармоника в 2 раза! Это я катушку датчика конкретно так нагружаю.... бум думать, тоже интересно...... Но не работает никто в токовом режиме, это я уже патент какой могу заявлять наверное
Гармоники проверил (на модели). Действительно зря боялся. На определение основной частоты при выбранном мной окне не влияют совсем: при 3-й,5-й 7-й гармониках с амплитудой равной основной гармонике рассчетная частота отличается от поданной примерно на 0.0006 Гц. То есть с этой стороны все в порядке.

А насчет той красивой картинки- так они точность 0.013% от значения обещают. даже для 100Гц это 0.013Гц, для более реальных 500Гц и выше- 0.065 Гц Не лучше чем у меня

Долго читал про потенциальную точность измерения, про то что здесь получается у народа (доли герца), но что-то с моими данными не могу сообразить.


1 вопрос. Имеются данные Нажмите для просмотра прикрепленного файла (после распаковки в матлабе загружаются load dat1). Как для данного конкретного случая (сигнал/шум, объем выборки, частота дискретизации 68кГц) оценить потенциальную точность измерения частоты? Формулу в предыдущих постах видел, но не пойму какое сигнал/шум подставлять, то ли во всей полосе, то ли только в полосе сигнала.

2. Насколько оптимально по точности измерение частоты следующим подходом, какую точность здесь можно ожидать?

Fs=68000;
N=512;
df=Fs/N;

sp=fft(s,N);

fmax=find(abs(sp(1:length(sp)/2))==max(abs(sp(1:length(sp)/2))));

f_ocen=(fmax(1)-1)*df;

f_s=f_ocen-df/1;
f_e=f_ocen+df/1;

m = 128*2;
w = exp(j*2*pi*((f_s-0)-(f_e+0))/(m*Fs));
a = exp(j*2*pi*(f_s-0)/Fs);

L1 = czt(s,m,w,a);
f_izm=find(L1==max(L1))*(2*df/m)+f_s;

Я, как (к сожалению) далекий от теории человек, делал проверку метода экспериментально на модели: задавал чистую синусоиду и вычислял разницу между этой подаваемой частотой и той, что мне метод насчитывал. В цикле для разных частот, потом на график. Если зависимость немонотонная- то детально обследовал внутри этой немонотонности (у меня, например, достаточно между двумя бинами ДПФ исследовать). Далее то же, но с использованием именно тех типов данных и размерностей, которые у Вас реально применимы в железе, не все и не везде double имеют. Если на этом этапе не нравится- то надо что-то в консерватории менять.
Ну а далее- добавляете мусору в исходный сигнал по вкусу- хоть гармоники, хоть белый шум, хоть еще что-нибудь, присущее вашей физике исходного сигнала.

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

Прежде чем экспериментировать и бросаться в пучину расчётов, рекомендую рассчитать два массива данных для интересующей длительности наблюдения и сравнить их.
1. Допустим, речь о частоте Fи
2. Желаемая точность dFи
3. Частота дискретизации Fд
4. Разрядность АЦП - N
5. Период наблюдения - T.
Таким образом, имеем T*Fд выборок, квантованных (в лучшем случае) по 2^(N-1) уровням...
Вычисляем 2 идеальных массива для частоты Fи и Fи + dFи и сравниваем их. В ряде случаев эти массивы (при чрезмерно оптимистичном подходе ) будут просто идентичными, тогда и копья ломать не стоит.
Если идти чуть дальше, то можно соотнести невязку этих двух массивов и одним из них с величиной, например, сигнал/шум реального АЦП. Это тоже покажет вам примерные идеальные возможности вашей системы.

Ruslan1, подскажите пожалуйста адреса ресурсов, где вы нашли функции по вычислению окна Гаусса и гаусовской интерполяции.
К примеру, на сайте http://www.dsplib.ru/content/winadd/win.html указана немного другая формула по вычислению коэффициентов окна Гаусса.

Просто я попробовал применить ваш алгоритм на С++ на реальном сигнале от гитары, так он выдает всегда результат интерполяции на 0,5 больше чем исходная частота.
В математике я слаб, но ваш метод интерполяции мне показался уж очень простым

http://www.nicholson.com/rhn/dsp.html - там есть подборка ресурсов на тему уточнения FFT
http://www.mathworks.com/help/toolbox/signal/gausswin.html


Да каких только формул я не видел, пока искал
Я уже писал, что не уверен в том, что мой вариант верный. Но матлаб для меня это авторитет, а у него внутренняя функция так же работает как я считаю.


0.5 Гц? странно, о такой разнице я еще не слышал. Можно исходники и окна и интерполяции в студию? И файл с исходным оцифрованным сигналом, очень интересно.....
Если на 0.05 Гц- это уже ближе к телу....

Сначала я в одном источнике прочитал что точнее будет, если логарифмическую кривую аппроксимировать параболой, а потом совсем в другом месте прочитал что это и есть гауссовская интерполяция.

Правильное Гауссовское окно считается так:
sqrt(sqrt(PI)) / sqrt(sigma) * sqrt(FSIZE) * exp (-(x-FSIZE/2.0)*(x-FSIZE/2.0)/2.0/sigma/sigma);
PI - это число пи, FSIZE - размер Фурье в точках, sigma - параметр, задающий ширину окна.
При таком нормировочном коэффициенте накладывание окна не приводит к дополнительному коэффициенту в амплитудах.
По поводу гауссовской интерполяции - все правильно, но полезно вводить весовые коэффициенты, пропорциональные квадрату амплитуды сигнала - улучшает точность, т.к. шум при малом сигнале меньше влияет на определение положения вершины.

Извиняюсь за очередной глупый вопрос: а как выбрать ширину окна (sigma)?
Например, хочу задавить третью гармонику, которая в худшем случае находится на расстоянии 105 бинов, если рассматривать частотную область. Это значит, что мне нужно взять окно, которое в временной области задавит находящийся на расстоянии 105 семплов от ценрального отчета (FSIZE/2) ? Или все совсем не так?
У меня получилось что при сигма=43 имею подавление 60 dB во временной области.

Или достаточно, чтобы к границам окна наблюдения (отчеты 0 и FSIZE) упало до нуля ( ну или -60 dB хотя бы) ?

Ширина выбирается в зависимости от задачи. Мне нужно было, чтобы вклад второй гармоники в амплитуду основного сигнала не превышал 1е-8, ну я и подбирал сигму так, чтобы посередине между основным тоном и второй гармоникой был провал не хуже 1е-8. Для определения только частот можно, думаю, ограничиться разделением на уровне 1е-6.
По поводу соотношения окна в частотной и временной областях - что-то там похожее должно быть, но я не исследовал точно. Может какой-то коэффициент вылезти. Лучше смотреть сразу в спектре.

Ruslan1, извиняюсь, при копи-пасте потерял еденичку в формуле интерполяции. Пока дебаггером не прошел, не замечал.
За ссылки - большое спасибо, немного подразобрался с теорией.

да пожалуйста.



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

Шкалу вертикальную сделайте линейную, а не логарифмическую - тоже ничего не увидите. У них же не децибелы, а амплитуда.

Строго говоря, приводит, т.к. коэффициент у вас рассчитывается для всей прямой, а не для ограниченного носителя.

Согласен, теоретически - да, но поскольку сигма выбирается так, чтобы значение на краях было очень мало, то ошибка тоже крайне невелика.

Вопросы о разрешающей способности измерения частоты за доли секунды с разрешением меньше Герца можно посмотреть здесь, очень невероятные вещи узнаете!
http://electronix.ru/forum/index.php?showt...mp;#entry703059

Для струны (струнные датчики растяжения) есть еще две проблемы, мы на них наезжали.

1. Частота колебаний зависит от их амплитуды. Грубо говоря, если струна колеблется с большей амплитудой - она больше растянута и ее жесткость выше. Единственный выход - минимизировать амплитуду ((или стабилизировать или выбирать участок с определенным ее значением).

2. В струне, особенно не очень корректно закрепленной, присутствуют колебания на нескольких очень близких частотах. Это связано, например с тем, что струна имеет чуть разные длины и жесткости для вертикального и горизонтального колебания, например на гитаре. А горизонтальная и вертикальная моды переходят одна в другую, в результате появляется медленная модуляция амплитуды на фоне общего затухания сигнала.

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

А наиболее правильный способ - или сделать на этой струне автогенератор, или принудительно возбуждать ее колебания внешним воздействием сканируя разные частоты. Последний способ - наиболее долгий (тут много и подводных камней) но точность возрастает на порядок. Не знаю применимо ли это к вашим датчикам. Основное преимущество - колебание становится непрерывным и амплитуда колебаний строго фиксирована. При этом спектр близких частот (разные моды колебаний - не путать с гармониками!) четко вычисляется по амплитудной модуляции за неограниченное время измерения.

Интересно, спасибо. А еще, мне всегда казалось, что струна выдает не чистую синусоиду, а все-таки с гармониками, как минимум с четными, то есть у струны отдельно есть колебания ее половины, четверти...

Ессно гармоники есть, но при определении частоты по фурье они не мешают. А вот близкие частоты мод возбуждения - мешают сильно, поскольку сбивают форму распределения (оно перестает быть гауссовским) и снижают точность определения частоты.

Т.е. если есть моды с частотами 1 кГц и 1.001 кГц, то погрешность определения частоты может достигать 1 Гц. Даже формально в этом случае непонятно, что называть резонансной частотой датчика. А подобное наблюдалось - неправильное крепление струны приводит к появлению сильно разнесенных по частоте мод. Спасает только селективное возбуждение только одной из них - а это уже не ударное возбуждение, а долгая (более 1 сек) накачка соответствующей частотой.

Если кому ещё интересна эта тема, вот хороший материальчик:
http://stackoverflow.com/questions/4633203...-between-frames

Метод уточнения частоты с использованием данных о фазе (нужно брать 2 фурье, лучше с нахлёстом):


лучше будет работать на искусственных линейно изменяющихся сигналах, но и для реальных есть смысл.

Хороший алгоритм по типу штангенциркуля. Только Вы забыли сказать, что такое k. Из-за этого, если не читать матерьяльчик, то непонятно, что здесь написано. Здесь k - это номер бина в котором максимум энергии. В матерьяльчике автор показывает, что правильную частоту можно получить не только по этому "максимальному" бину но и по FFT_N/FFT_STEP соседних. Поэтому не возникает проблем с частотами на границе между бинами при перекрытии блоков FFT. Без нахлёста с такими пограничными частотами алгоритм работать не будет, шумы не дадут. Систематическая ошибка у алгоритма отсутствует. Интересно посмотреть как алгоритм ведёт себя по отношению к шуму - судя по всему выйдет на точность оценки максимального правдоподобия CRLB, особенно если ещё и научиться делать взвешенную оценку по нескольким значениям возле максимума (остальные FFT_N/FFT_STEP оценок из-за малой энергии бинов реально утонут под шумом)

А можно этот код расписать иначе? Так как в матлабе не силён. Если можно то блок схемой или просто словами.

Я планирую выложить примерчик с интересными методами здесь. Но думаю это будет не скоро. А пока попробую на словах, например, для фурье на 1024 отсчёта:
1) Берём FFT и получаем массивы амплитуд мнимой и действительной части по 512 элементов в каждом Im[0..511], Re[0..511].
2) Для каждого из 512 элемента считаем фазу. Если по простому, то фаза[i]:=ArcTan(Im[i] / Re[i]); Можно сделать получше, расписав разные случаи (если Re[i]=0 и т.д.)
3) храним фазы за предыдущий проход FFT и за текущий. m_fftLastPhase[i] - это старая фаза, а phase - это на самом деле phase[i], т.е. фаза[i] от текущего FFT. Заметил, что метод лучше работает если применить сглаживающее окно.

Забыл написать: например, константа SAMPLE_RATE=44100, FFT_N =1024, FFT_STEP = 100 (количество точек, на которые смещаем исходные данные. 100 намного меньше 1024, поэтому результат лучше)

И всё таки не много не понятно.
Что делает эта функция?

вычисляет остаток от деления дельты на два пи, вероятно...

Фазу приводят к диапазону -3.14...+3.14, пример в градусах: это что-то наподобие того, что 350 градусов привратится в -10 и т.д. Т.е. вектор (Im,Re) всегда в правом полукруге от -90 до 90 градусов.

Именно на этом этапе алгоритм перестанет работать, если нет перекрытия. Как показывает автор "матерьяльчика" в последней табличке (Pass 4) для отсутствия перекрытия.
http://www.dspdimension.com/admin/pitch-sh...g-using-the-ft/
Для частот на границе бинов получаем 2 значения частоты для примерно одинаковых величин спектральной мощности, и только одно из них правильное. Поскольку максимум двух равных величин при реальных вычислениях определяется шумом, то в этом случае алгоритм будет выбирать правильное значение с вероятностью 0.5. Впрочем можно попытаться определить частоту другим методом по каждому из блоков отдельно и потом взять то значение которое ближе к этой оценке, но это уже как бы гибрид



Для определения частоты вам не нужны все FFT_N значений фазы, поскольку только FFT_N/FFT_STEP значений фазы около максимума спектральной мощности правильные, а остальные искажены врэпингом на 2*pi (Pass 5). Более того в реальной ситуации присутствуют шумы и большинство из этих FFT_N/FFT_STEP правильных оценок (если их много) будут реально сильно зашумлены из-за очень быстрого падения уровня сигнала в бине при отклонении от центрального. Мы помним, что отклик бинов на синусоиду падает в соответствии с функцией спектрального окна, и в лучшем случае он будет спадать как синк или быстрее. В реальной ситуации только 2-4 центральных бина будут давать значимую оценку частоты, причем наименее искажены шумом будут только одна (если частота близко к центру бина) или две (если частота на границе бинов) оценки для центральных бинов (где максимум спектральной мощности). Остальные "утонут" под шумом
Все 512 значений могут представлять интерес только для визуализации "сонограм", поскольку они, в большинстве, шумоподобны.
Поэтому приведеный алгоритм "недоделаный" - мало того что он неэффективен, но он не содержит последнего шага, поскольку он генерирует FFT_N значений частоты из которых в лучшем случае (отсутствие шума) только FFT_N/FFT_STEP правильны и не сказано, как правильные выбрать из массива оценок, большинство из которых неправильны, см. таблички у автора
http://www.dspdimension.com/admin/pitch-sh...g-using-the-ft/
По смыслу понятно, что выбирать нужно где-то возле максимума спектральной мощности

С учетом сказанного простейший алгоритм будет выглядеть так:
1) Берём FFT и получаем массивы амплитуд мнимой и действительной части по 512 элементов в каждом Im[0..511], Re[0..511].
2) Вычисляем энергию для каждого k E[k]=Re*Re+Im*Im
3) Находим максимум E[k] - пусть будет E() для kmax. Здесь вообще-то нужно проверять что kmax не прыгает. Если изменение kmax от фрейма к фрейму больше 1 по модулю - то это не синусоида, выход по ошибке
4) Для элементов kmax, kmax-1, kmax+1 считаем фазу. фаза[]:=ATAN2(Im[], Re[]); другие элементы на самом деле не нужны, в них отсутствует полезная информация
3) храним посчитаные фазы за предыдущий проход FFT и за текущий. m_fftLastPhase[i] - это старая фаза, а phase - это на самом деле phase[kmax], т.е. фаза[kmax] от текущего FFT. Поскольку kmax гарантировано не меняется больше чем на 1, то в любом случае значение фазы для kmax посчитаны и для предыдущего фрейма и есть в массиве

Если перекрытие FFT_N/FFT_STEP >2 то оценку можно немного улучшить если использовать не одно значение максимума, а два первых максимума и построить взвешенную оценку для частоты по этим двум оценкам.
Если главный максимальный бин k0=kmax, а следующий по энергии k1=arg max(E(kmax-1), E(kmax+1)), то можно построить оценку типа
(f(k1)*E(k1)+f(k0)*E(k0))/(E(k0)+E(k1)) с большим иммунитетом к шуму, главным образом для частот на границе бинов (поскольку в этом случае будет две равноценных оценки)

А так да, идея алгоритма на вскидку смотрится очень хорошо

f(k1) это измеренная частота для бина k1 ? А что если приращение фазы измерять по нескольким бинам в окрестности
максимума? Взвешивание получится автоматически.



X_1 - предидущий выход FFT (вектор столбец),
X - текущий выход FFT,
kmax - индекс бина с максимальной энергией,
n - число 0...3

Можно бы и по всем. только не все бины одинаково полезны. Поэтому это будет не оптимально.
Усреднять надо в соответствии с уровнем сигнала каждого бина. Вычисление фазы в бинах с низким уровнем сигнала по отношению к шуму бесмысленны - это уже фаза компоненты шума а не сигнала. А уровень сигнала убывает при отклонении от центральных бинов достаточно быстро (быстрее чем 1/n), в то время как шум остается постоянным. Поэтому при низком уровне сигнал/шум лучше усреднять по 2-3-м, причем с учетом энергии (можно даже расписать точные формулы, выписывая плотность вероятности ошибки фазы в бине как функцию отношения сигнал/шум - это можно найти в учебниках по спектральному анализу или статрадиофизике - и максимизмизируя произведение вероятностей).

При очень высоком уровне сигнал/шум, хотя бы 20дб, действительно усреднять вроде можно по всем четырём (FFT_N/FFT_STEP в общем случае перекрытия, но тоже с учетом энергии бина), поскольку ошибки измерения фазы статистически независимы, но различны для разных бинов, в зависимости от энергии

да, f(k) - измеренная частота для бина k. Я не максимизировал функцию вероятности, а просто написал на вскидку формулу, которая будет вести себя оптимально для частот и посередине бина (f(k0)) и на границе бинов (f(k0)+f(k1))/2,
если этого не сделать, то точность упадёт на границе

ЗЫ. Вообще-то вычисление фазы через квадратуры всегда усиливает шум и помехи. Поэтому при вычислении частоты через фазу не стоит пренебрегать использованием спектральных окон, как в примере - чтобы подавить хвосты от возможных синусовых помех. Они и всегда мешают, а этому методу будут особенно

Так предлагаемое скалярное произведение и учитывает энергию бинов (может быть даже оптимальным образом ).
Бины с большей энергией вносят больший вклад в результирующую фазу, чем бины с меньшей энергией.

Виноват. Про angle (в Матлабе?) не знал, думал среднее арифметическое фазы.

Сорри, не написал что матлаб. angle -аргумент комплексного числа.

Более правильно - по МНК подбирается параметры гаусса (огибающей) (при соответствующей весовой функции) с весами, пропорциональными мощности. При этом точность возрастает на порядок.

Так предложили уже взвешивать, если говорить о алгоритме с разностью фаз ДПФ перекрывающихся фреймов.
Или Вы о чём то своём?