Jump to content

    
p_v

Проверьте пожалуйста алгоритмы IIR фильтров

Recommended Posts

https://github.com/speedcontrols/dc_sc_grinder/tree/scripts/scripts/lib

 

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

 

Обычный батерворт 2 порядка, из которого лепится потом каскад 8-12 порядка. Исходник пока на яваскрипте, чтобы по-быстрому алгоритм поиска фундаментальной частоты отладить.

Edited by p_v

Share this post


Link to post
Share on other sites
17 минут назад, p_v сказал:

https://github.com/speedcontrols/dc_sc_grinder/tree/scripts/scripts/lib

 

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

 

Обычный батерворт 2 порядка, из которого лепится потом каскад 8-12 порядка. Исходник пока на яваскрипте, чтобы по-быстрому алгоритм поиска фундаментальной частоты отладить.

 

2-й порядок. Частота среза  -3дб. 4-й -6дб. 8-й -12дб итд. Так и задумано? Если хочется -3дб на частоте среза при любом порядке нужно синтезировать фильтр для заданного порядка, а не городить их из 1-й биквадратной секции.

Share this post


Link to post
Share on other sites
47 minutes ago, thermit said:

2-й порядок. Частота среза  -3дб. 4-й -6дб. 8-й -12дб итд. Так и задумано? Если хочется -3дб на частоте среза при любом порядке нужно синтезировать фильтр для заданного порядка, а не городить их из 1-й биквадратной секции.

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

 

Фильтрами надо просто найти нужную полосу, где от сигнала останется синус, и дальше уже точную частоту определить через zero cross. Как там заваливаются края полосы и какую точку считать частотой среза - не принципиально.

 

На глаз, для исходного сигнала надо 8-12 порядок батерворта, чтобы лишнее убрать.

Edited by p_v

Share this post


Link to post
Share on other sites

Частоту среза по -3дБ считаете? При последовательном включении фильтров она будет изменяться. Т.к. в точке где один фильтр давал -2дБ, два дадут уже -6дБ. Т.е. частота среза (точка -3дБ) сместится.

Share this post


Link to post
Share on other sites
1 minute ago, arhiv6 said:

Частоту среза по -3дБ считаете?

Если речь об исходниках, которые я прошу проверить, то да, там стандартная терминология, -3дБ.

 

22 minutes ago, arhiv6 said:

Т.к. в точке где один фильтр давал -2дБ, два дадут уже -6дБ. Т.е. частота среза (точка -3дБ) сместится.

Ну вроде это не должно стать большой проблемой. Мне фильтры из каскадов второго порядка нужены только для выделения чистого синуса, а точное значение я все равно по zero-cross буду определять.

 

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

Share this post


Link to post
Share on other sites

Формулы странные.

Должно быть так:

w=tan(pi*f0/sf);

norm=1.0/(w*(w+sqrt(2))+1.0);

a1=2.0*norm*(w*w-1.0);

a2=norm*(w*(w-sqrt(2))+1.0);

lpf

b1=2*norm*w*w

b2=b0=0.5*b1

hpf

b1=-2*norm

b0=b2=-0.5*b1

 

Разностное уравнение

u(n)=x(n)-a1*u(n-1)-a2*u(n-2);

y(n)=b0*u(n)+b1*u(n-1)+b2*u(n-2);

 

Edited by thermit

Share this post


Link to post
Share on other sites
7 hours ago, thermit said:

Формулы странные.

Напарник откуда-то сам выводил, поэтому немного не классическая запись получилась. У нас разделение труда, я не лезу в математику/физику пока совсем не припрет.

 

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

 

 

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

 

На первый взгляд - ничего такого, что нельзя было бы почистить простенькими фильтрами. Частота сигнала в диапазоне 500-15000 герц. Дискретизация сейчас 250 килогерц, чтобы пока не возиться с пред-фильтрацией, в реальном девайсе в итоге будет ~50 килогерц.

Edited by p_v

Share this post


Link to post
Share on other sites

Переделал фильтры на "канонические". Принципиально ничего не поменялось. Но теперь хотя бы нет сомнений, что фильтры правильные. @thermit, еще раз спасибо.

 

Это после ФНЧ 15 кГц (батерворт, 2 порядок, 4 раза)

lp_15k.thumb.png.569a0bfa8ea8d45386e68857b7c246e4.png

 

Это после ФНЧ 15 кГц и ФВЧ 300 Гц.

lp_15k_hp_300.thumb.png.a41a7b125e7128916284930d1bb29a8b.png

 

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

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

Дальше строим несколько сеток полосовых фильтров, уменьшая шаг сетки, и "половинным делением" подбираем полосу, когда сигнал станет больше какого-то порога. Ширина полосы должна быть ~0.7 в каждую сторону от середины, чтобы гармоники не прихватить. В фурье перегонять не хочется - дорого.

 

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

 

Осталось понять, как вычислять порог полезного сигнала (когда с выхода полосового фильтра фундаментальная частота поперла). Из того что приходит в голову - на исходном сигнале взять положительную часть, посчитать среднее арифметическое (что-то типа средней амплитуды), и порогом принять 25-50%. Я в привильную сторону думаю?

 

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

Edited by p_v

Share this post


Link to post
Share on other sites

1-й вопрос, что в голову приходит: а с какой точностью нужно измерять?

Проверенные исходники? fftw. Правда с плавающей точкой. В фиксированной точке сходу могу предложить исходники от ti.  Это все вычисление дпф. На счет непосредственно вычисления сонограммы - такого не встречал. Придется допиливать самостоятельно. Да там ничего сложного. Кстати, что за контроллер?

Share this post


Link to post
Share on other sites
14 часов назад, p_v сказал:

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

 

Здравствуйте!

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

я скачал файл hilda_50kHz_rpm_high.txt и попробывал применить медианный фильтр.

оригинальный сигнал:

1154958688_.thumb.png.9442dd3f01df5fe04cca106ccbf9656f.png

результат обработки медианным фильтром третьего порядка:

1960551823_.thumb.png.c88b869528b03390c46fd32cf103427c.png

результат обработки медианным фильтром пятого порядка:

1517149767_.thumb.png.cfff8bbbbbf79c567a81e380461f27f5.png

по результатам обработки видно, что тут можно рассмотреть применение медианного фильтра до фильтрации ФНЧ и ФВЧ.

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

 

Share this post


Link to post
Share on other sites
11 minutes ago, thermit said:

1-й вопрос, что в голову приходит: а с какой точностью нужно измерять?

Ну если быть реалистом, то для юзеру отклонение оборотов бормашинки на 10% глубоко пофик. Если учесть запас для работы  ADRC  т.п. - погрешность меньше 1% точно не нужна.

 

12 minutes ago, thermit said:

Проверенные исходники? fftw. Правда с плавающей точкой. В фиксированной точке сходу могу предложить исходники от ti.  Это все вычисление дпф. На счет непосредственно вычисления сонограммы - такого не встречал. Придется допиливать самостоятельно. Да там ничего сложного.

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

 

В идеале ведь сигнал содержит частоту, пропорционалдьную оборотам, и 100 герц. Поэтому теоретически достаточно через ФВЧ полностью завалить 100 герц, а дальше сделать сетку ФНЧ в диапазоне 500-15000 герц, и подобрать ближайший к нулю, который начнет пропускать сигнал (это и будет наша искомая частота).

 

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

 

Можно наверное то же самое сделать через DFT, забив несколько таблиц с шагом 2, 4, 8, ..., и "уточняя" только тот участок где найден максимальный пик спектра. Но по-моему с фильтрами проще.

 

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

 

28 minutes ago, thermit said:

Кстати, что за контроллер?

https://oshwlab.com/speed/dc-speed-control

 

Cortex-M0+. В наличии только умножение 32*32. То есть, вычисления с фиксированной точкой 16.16 (или может быть 8.24). Отсюда и батерворты 2 порядка, объединенные последовательно. Другие конфигурации от нехватки точности расколбасит.

5 minutes ago, stealthisname said:

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

Спасибо что не пожалели времени.

 

IMHO в медианном фильтре не будет особого смысла, т.к. с теми пиками, которые он способен вырезать, ФНЧ высокого порядка справится не хуже. А он намного проще. Ну и как я писал выше, если процессить не все полуволны целиком, а у каждой похерить 25% с начала и 25% с конца, то прострелы станут не особо критичными.

 

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

Share this post


Link to post
Share on other sites

Прогнал сэмплы через эту библиотеку https://github.com/indutny/fft.js, 16834 отсчетов, частота 250 килогерц.
 

const FFT = require('fft.js')
const fft_size = 16384

const fft_input = Array.from(data.subarray(4000, 4000 + fft_size))
const fft_out = new Array(fft_size)

const f = new FFT(fft_size)
f.realTransform(fft_out, fft_input);
f.completeSpectrum(fft_out)

Не уверен что правильно использовал, на выходе получается 32К точек, на графике первые 2000 (дальше либо нули либо мусор):

spectre1.thumb.png.bd6d72f13750c165624173075abbd933.png

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

Есть идеи что подкручивать?

Share this post


Link to post
Share on other sites
05.10.2021 в 10:22, p_v сказал:

Фильтрами надо просто найти нужную полосу, где от сигнала останется синус, и дальше уже точную частоту определить через zero cross

 

18 часов назад, p_v сказал:

На первый взгляд - ничего такого, что нельзя было бы почистить простенькими фильтрами. Частота сигнала в диапазоне 500-15000 герц.

 

6 часов назад, p_v сказал:

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

 

4 часа назад, p_v сказал:

В идеале ведь сигнал содержит частоту, пропорционалдьную оборотам, и 100 герц. Поэтому теоретически достаточно через ФВЧ полностью завалить 100 герц, а дальше сделать сетку ФНЧ в диапазоне 500-15000 герц, и подобрать ближайший к нулю, который начнет пропускать сигнал (это и будет наша искомая частота).

 

Вот часть спектра для сигнала hilda_50kHz_rpm_high.txt. Просто тупо FFT по всему сигналу, без всяких предварительных обработок. Ну и где здесь искомая частота, пропорциональная оборотам?:) Пусть даже без учета леса гармоник от 100Гц...

High_0_4400.thumb.jpg.60eb41e4410a0b023486d726d9c8f6e6.jpg

Share this post


Link to post
Share on other sites

Какая-то у вас ерунда со спектром. Наверное надо все-таки 100 герц и их гармоники сначала обкоцать. Еще в бренче "raw_data" лежат семплы на 250 килогерц. Лучше пока их юзать, чтобы не было всяких левых отражений. 50 килогерц я сделал простым прореживанием, это не совсем корректно.

 

Вроде разобрался с FFT. Выборка 250 килогерц, 16К точек, на графиках первые 1024 точки

 

 

abs(real)

spectre_re.thumb.png.f60be27d6e0041a8a9fed5f03b6a4d83.png

 

 

sqrt(real^2 + imaginery^2)

spectre_re_img.thumb.png.82ea15c8b08349f8a750920d24ac1d28.png

 

Не знаю, по какому правильнее смотреть, вроде по первому?

 

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

Edited by p_v

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.