[Zuse 0]

Всем привет!

Пытаюсь сделать на МК 1986ВЕ92 (близок к STM32F103) устройство, которое бы непрерывно оцифровывало поступающий на вход АЦП в диапазоне нескольких кГц сигнал, выделяло в сигнале гармоническую составляющую с наибольшей амплитудой и формировало с помощью ЦАП синус с той же частотой и фазой. 

Для этого каждые 5 мкс АЦП непрерывно производит АЦ-преобразования, получаются массивы по 1024 отсчета, которые фильтруются цифровым ФНЧ, затем делается БПФ по 256 отсчетам (используются каждый четвертый отсчет из 1024). Далее методами интерполяции, описанными по ссылке 

https://www.dsprelated.com/freebooks/sasp/Quadratic_Interpolation_Spectral_Peaks.html   

определяется частота и фаза спектрального пика, рассчитывается массив из 1024 значений выходного синуса, который через DMA скармливается ЦАП.

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

Коллеги, подскажите, плиз, где можно почитать поподробнее про определение фазы спектрального пика?    

       

[_pv 105]

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

ну и вместо интерполяции фазы, можно интергалы (*sin(w0*t), *cos(w0*t)) на найденной частоте с максимальной амплитудой пересчитать честно, их арктангенс == фаза.

[petrov 18]

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

Вот пример использования скользящего ДПФ для выделения гармоники несущей частоты, чтобы на приёме убрать вращение созвездия с разностной частотой.

https://electronix.ru/forum/topic/23652-model-8psk-modema/page/13/#comment-1755048

[Zuse 0]

13 hours ago, _pv said:

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

А это не одно и тоже ?

[vervs 87]

36 минут назад, Zuse сказал:

А это не одно и тоже ?

гармоники компоненты сигнала кратные основной частоте, видимо _pv не понял, что Вы имели в виду бины спектра, возле интересующей частоты
кстати учитывайте фазу ФНЧ, а при интерполяции еще и окно

[blackfin 76]

On 1/14/2026 at 5:45 PM, Zuse said:

Для этого каждые 5 мкс АЦП непрерывно производит АЦ-преобразования, получаются массивы по 1024 отсчета...
В результате вижу, что частота определяется достаточно точно, ... а вот с фазой проблемы.

Частота, это производная фазы по времени. Поэтому ошибка измерения фазы "Δφ" равна интегралу по времени от ошибки измерения частоты "ΔF".
На интервале измерения равном "T" это дает ошибку: Δφ = 2*pi*ΔF*T.
Поэтому, если допустимая ошибка измерения фазы равна "Δφ", то необходимо измерять частоту "F" гармонического сигнала с точностью: ΔF ≤ Δφ/(2*pi*T).
Это условие выполняется?

[Zuse 0]

11 minutes ago, blackfin said:

Частота это производная фазы по времени. Поэтому ошибка измерения фазы "Δφ" равна интегралу по времени от ошибки измерения частоты "ΔF".
На интервале измерения равном "T" это дает ошибку: Δφ = 2*pi*ΔF*T.
Поэтому, если допустимая ошибка измерения фазы равна "Δφ", то необходимо измерять частоту "F" гармонического сигнала с точностью: ΔF ≤ Δφ/(2*pi*T).
Это условие выполняется?

Чему равен Т ?

[blackfin 76]

On 1/15/2026 at 9:20 AM, Zuse said:

Чему равен Т ?

Из первого поста понятно, что T равен:

On 1/14/2026 at 5:45 PM, Zuse said:

1024*5 мкс.

 

[Zuse 0]

1 hour ago, blackfin said:

Очевидно, что T равен...

По приведенному отношению получается, что при ошибке измерения частоты 5 Гц (у меня что-то типа того) ошибка по фазе будет составлять 9 градусов. Думаю, похоже на то, что вижу на осциллографе.

Но меня не точности сейчас интересуют, а интересует правильность метода определения фазы, который я описал.

Дело в том, что когда я подаю на вход моего устройства синус частотой допустим около 4 кГц, частоту которого можно немного варьировать, и пытаюсь интерполировать фазу спектрального пика по двум наибольшим бинам, то фазовый сдвиг между входным синусом и синусом на выходе ЦАП стабилен только на частотах более-менее кратных основной частоте 1/(1024*5 мкс). При отклонении от кратных частот фазовый сдвиг начинает дрожать/скакать тем сильнее, чем больше отклонение. 

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

Чтобы себя перепроверить я попробовал повторить алгоритм интерполяции фазы в Excel. Сформировал массив значений синуса с частотой примерно в 10 раз большей, чем основная частота, но не кратной основной частоте. Затем отдельно для каждой четверти массива сделал БПФ, потом для каждой четверти с помощью интерполяции по двум наибольшим бинам определил фазу спектрального пика и в конце попробовал восстановить исходный сигнал, заново рассчитав значения синуса внутри четвертей, используя полученные с помощью интерполяции четыре значения фаз. В результате первая и вторая четверти друг с другом по фазе не состыковавались.

Короче,  возникли у меня сомнения по части правильности алгоритма в результате.

 

 

[dOb 10]

БПФ используем для первого приближения, а дальше методом наименьших квадратов подбираем синус.

Правда вопрос: хватит ли быстродействия. 

[blackfin 76]

On 1/15/2026 at 10:28 AM, Zuse said:

Но меня не точности сейчас интересуют, а интересует правильность метода определения фазы, который я описал.

Метод, очевидно, неправильный..

1. Децимация по времени не нужна, поскольку ухудшает точность измерения частоты.
2. Цифровой ФНЧ не нужен, поскольку вносит ошибку в фазу измеряемого сигнала даже для КИХ-фильтра с линейной ФЧХ.
3. Перед вычислением БПФ входной сигнал нужно умножить на оконную функцию, напр. на: gausswin. Это уменьшит влияние боковых лепестков sinc-функции прямоугольного окна.
4. Для увеличения точности измерения частоты можно сделать интерполяцию в частотной области добавлением нулей и вычислением БПФ от массива размером 4096 точек или больше.
5. После этого фазу гармонического сигнала можно посчитать по формуле: \(\varphi_s=\arg(C_k)-(\omega_s-\omega_k)\cdot\frac{T}{2}\)

 

Где:
Ck - комплексный бин БПФ с максимальной амплитудой;
ωs - частота гармонического сигнала найденная интерполяцией по трем точкам;
ωk - частота бина БПФ с максимальной амплитудой;

T - интервал измерения.

 

Подробнее можно почитать по ссылкам: Частота сигнала, Измерение частоты основной гармоники.

[Zuse 0]

9 minutes ago, blackfin said:

Метод, очевидно, неправильный

Неконкретно сформулировал вопрос. Переформулирую: линейная интерполяция фазы двух максимальных бинов даст фазу спектрального пика?

PS Все прочие детали (ФНЧ, децимация итд) приведены для формирования общей картины

[_pv 105]

32 minutes ago, Zuse said:

Неконкретно сформулировал вопрос. Переформулирую: линейная интерполяция фазы двух максимальных бинов даст фазу спектрального пика?

фаза там между бинами настолько же линейна, насколько максимум похож на параболу.

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

при непопадании длины преобразования в период, хвосты от несовпадения на концах этот макисмум размажут и сместят (пусть и меньше с оконной функцией), заодно и фазу.

добавте оконную функцию и уберите фнч, возможно этого хватит.

[Alex11 13]

Вот Вам пример как точно считать частоту и фазу.

Fourier.zip

Здесь считается для 50 Гц и окрестностей с БПФ на 8192 отсчета, возможно, при коротком БПФ нужно будет что-то поменять.

И при гауссовском окне пик не парабола. Парабола - логарифм от пика.

[petrov 18]

Zuse

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

И не надо ничего интерполировать, ФЧХ полосовых фильтров в ДПФ не выровнены по фазе, при переключении между бинами происходит скачёк фазы, нужно выход каждого фильтра умножить на соответствующий поворачивающий коэффициент и скачков не будет.