Jump to content

    
p_v

Подскажите с алгоритмом определения фундаментальной частоты

Recommended Posts

Задача: хочется определять обороты коллекторного двигателя по пульсациям тока, которые вызывает переключение обмоток. И чтобы математика на Cortex M0+ вписывалась в 0.1 сек за итерацию.

  • Искомая частота - 600...6000 герц (450...45000 rpm для 8-полюсного коллекторника бормашинки).
  • Желаемая точность измерений - 1% (относительно измеряемой частоты)

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

 

Что имеем:

  • https://github.com/speedcontrols/dc_sc_grinder/tree/scripts/scripts - исходники, с которыми экспериментирую (в бранчах есть полные дампы).
  • Написал честный FIR-дециматор, и перегнал дампы на 16 килогерц (примерно как в реальном девайсе)
  • Дампы и прочие выклопы нормализовал в целочисленные.

Дамп на высоких оборотах

high_16khz.thumb.png.e7453b409129f5fce770d15c9f747f86.png

 

Дамп на низких оборотах

low_16khz.thumb.png.8176ca34584fcf5f3d98f263adf8d2e5.png

 

Спектр на высоких оборотах, 512 точек, обнулена постоянная составляющая (чтобы не портила масштаб). ~ 32 герца на точку. Нужен второй большой пик (который "троится"). ~ 4000 герц, похоже на правду (30000 rpm).

high_spectre.thumb.png.d122613a896c6c06a0cd6c44ded9d79c.png

 

Спектр на низких оборотах, 512 точек, обнулена постоянная составляющая. Тоже второй пик, ~ 600 герц (4500 rpm), сходится.

low_spectre.thumb.png.33b618772bce08896f0944d74a9cd941.png

 

С точностью на низах примерно понятно. Если видим что пик внизу - включаем ресамплер, и делаем FFT еще раз. Сетевые пульсации и гармоники в стороне от нужного диапазона - можно спокойно игнорировать. Правда, критично вырастает время выборки. Получается 0.8 сек при точности 1.3%. Терпимо. В крайнем случае похерю точность до 2%, приподняв частоту самплинга.

 

Остались вопросы:

  • Что все-таки смотреть после FFT, реальную часть или корень из суммы квадратов?
  • Просто искать максимум или как-то проверять соседние отсчеты?
  • Как определять "отсутствие" оборотов (или когда ниже допустимых, что актуально для старта). Ведь если PID/ADRC подхватит ложный пик - может не дать разогнаться. Мерить при калибровке амплитуду частоты на низах и ставить порог в половину?
  • Есть ли смысл накладывать перед FFT оконную функцию?
  • Может еще чего-то не знаю и пропустил?
Edited by p_v

Share this post


Link to post
Share on other sites
33 minutes ago, p_v said:

Остались вопросы:

  • Просто искать максимум или как-то проверять соседние отсчеты?
  • Есть ли смысл накладывать перед FFT оконную функцию?
  • Может еще чего-то не знаю и пропустил?

 

И тд..

Share this post


Link to post
Share on other sites

Может быть вам поможет автокорреляция?  Попробуйте измерить время от пика до первого перехода через 0, как четверть периода частоты

Share this post


Link to post
Share on other sites
10.10.2021 в 08:37, p_v сказал:

Что все-таки смотреть после FFT, реальную часть или корень из суммы квадратов?

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

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

Share this post


Link to post
Share on other sites
On 10/10/2021 at 8:37 AM, p_v said:

И чтобы математика на Cortex M0+ вписывалась в 0.1 сек за итерацию.

А сейчас сколько времени эти все вычисления выполняются?

И зачем вам, собственно, БПФ? Попробуйте дискретные косинусные преобразования. Правда, флеш сильно раздуется таблицами косинусов…

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

P.S. А еще мне непонятно, зачем мучиться с определением частоты оборотов по пульсациям, если у вас УЖЕ есть устройство, которое эту самую частоту задает? Ну так и берите данные о скорости вращения со своего привода!..

Edited by Eddy_Em

Share this post


Link to post
Share on other sites
10.10.2021 в 08:37, p_v сказал:

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

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

в качестве примера можно попробовать такой алгоритм на сигнале из Вашего файла hilda_15625Hz_zero_to_middle.

я вычислил fft для 512 отсчетов на протяжении всего файла. для каждого результата убрал первые несколько отсчетов - в них мы ожидаем сильную низкочастотную помеху.

получилась такая спектрограмма.

1580207989_.thumb.png.f70c68946d30d54aa6c927843ea609e6.png356960761_.thumb.png.a9aa2b174ec1696571e3bdac07cd7b4a.png

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

в начале файла шум не превышает порог

1931207851_.thumb.png.0d0ae62d4ab920ba63d4ea52dc34dc0d.png

во время включения уже виден нужный пик, но частота слишком быстро изменяется на интервале одного fft, поэтому порог не превышен

609641312_.thumb.png.70bc36724e4e598312fb39bd00f1ff1f.png

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

1786396732_.thumb.png.98c8b8ab32419bbd9c095e5d62643910.png

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

 

вот как выглядит спектрограмма с порогом

977232755_.thumb.png.6aa997591994781a422f2a1692a02f22.png1495844636_.thumb.png.868aa84fada294cb76eb71285da40a31.png

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

результат работы алгоритма - оценка частоты

1368709200_.thumb.png.19440211771dcb8d2c55929d14b3334a.png198641435_.thumb.png.baf205553ff7b5346e5b7bfa9c3795cd.png

аналогичные результаты для сигнала из файла hilda_15625Hz_zero_to_low  с тем же размером fft и с тем же множителем при вычислении порога

143583384_.thumb.png.18c578ae2674ee1bab858125aeee4448.png2037838990_.thumb.png.a9748a6a60273be8257f44d8b7be9682.png

 

228477947_.thumb.png.76fdfa9e7448364100e7db67659998ef.png163997626_.thumb.png.0efed0f06dc9ff6673f459ece8586b6f.png

Edited by stealthisname

Share this post


Link to post
Share on other sites

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

image.thumb.png.0f8328422ddf5dc198860b98815bd14d.png

просто пропустив через ФВЧ чтобы отрезать 50Гц, а от остального посчитать отношение среднеквадратичного значения производной к среднеквадратичному значению сигнала, и среднеквадратичное значение сигнала в качестве порога и использовать.

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

Share this post


Link to post
Share on other sites

ФВЧ можно и аппаратный поставить перед АЦП, сильно упростив себе работу (да и МК меньше придется телодвижений совершать). Про использование более шустрых ДКП вместо БПФ я уже выше говорил. Но сдается мне, какой-нибудь простой вариант периодограммы будет еще быстрей.

Share this post


Link to post
Share on other sites
On 10/12/2021 at 12:28 PM, _pv said:

просто пропустив через ФВЧ чтобы отрезать 50Гц, а от остального посчитать отношение среднеквадратичного значения производной к среднеквадратичному значению сигнала, и среднеквадратичное значение сигнала в качестве порога и использовать.

Я не понял, что это за магия, и почему она вдруг стала пропорциональна оборотам, а не погоде на марсе например.

 

On 10/12/2021 at 7:21 AM, stealthisname said:

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

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

 

Share this post


Link to post
Share on other sites
1 hour ago, p_v said:

Я не понял, что это за магия, и почему она вдруг стала пропорциональна оборотам, а не погоде на марсе например.

потому что синус отличается от производной синуса по амплитуде ровно в w=2*Pi*f раз. то есть их отношение поделённое на 2Pi и есть частота.

ну а через среднеквадратичное значение просто удобно амплитуды синусов считать, независимо от частоты.

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

Share this post


Link to post
Share on other sites

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

 

Скажем так, скорость самого БПФ на 1024 точки не самый критичный кусок. Там больше прижимает время набора отсчетов (1000/16000 => 0.0625 сек). Я пока попробую на БПФ, т.к. там предсказуемый результат, но если что-то пойдет не так - буду думать над вашим.

 

PS. Ну и что скрывать, просто интересно БПФ поюзать.

Edited by p_v

Share this post


Link to post
Share on other sites
12.10.2021 в 14:46, Eddy_Em сказал:

ФВЧ можно и аппаратный поставить перед АЦП

А что, разве можно его не ставить? Наложение спектров (aliasing) придумали слабаки и трУсы?

Share this post


Link to post
Share on other sites
9 hours ago, p_v said:

Плюс подсчет девиации - не самая дешевая штука.

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

БПФ ~100 тактов на отсчёт.

Share this post


Link to post
Share on other sites
5 hours ago, Сергей Борщ said:

А что, разве можно его не ставить? Наложение спектров (aliasing) придумали слабаки и трУсы?

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

On 10/12/2021 at 1:46 PM, Eddy_Em said:

ФВЧ можно и аппаратный поставить перед АЦП, сильно упростив себе работу

Аналоговый фильтр перед АЦП не сильно поможет, так как фильтр там на самом деле довольно злой должен быть (а не скользящее среднее как в примере на картинке с частотой среза 1.5кГц), у ТС частоты от 600Гц начинаются и он 1% хочет, при этом Гц до 200-300 надо бы всё задавить довольно хорошо. Ну то есть можно конечно попробовать такой фильтр сгородить снаружи, но это совсем не упростит, а скорее наоборот.

Share this post


Link to post
Share on other sites
2 hours ago, _pv said:

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

Если быть точнее, там выпрямленное напряжение, поэтому 100 либо 120 герц, плюс конская вторая гармоника, и возможно приличные третяя-четвертая.

Поэтому надо сурово резать от нуля до 400 - 480 герц. И учитывать от 600 герц. Чтобы такое разделить - это от десятого порядка фильтры. Без плавучки - либо каскады IIR либо очень длинный FIR. Можно, но пипец как грустно.

 

7 hours ago, _pv said:

БПФ ~100 тактов на отсчёт.

Часто встречал такую оценку на этом форуме. А для каких она условий? Случайно не для аппаратной плавучки?

 

https://community.st.com/s/question/0D50X00009XkhK9/fft-in-stm32f0 - тут в четветом ответе пишут для M0+ 600К циклов на 1024 точки для Q15 (мне возможно понадобится Q31, т.к. постоянная составляющая очень большая).

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.