Перейти к содержанию
    

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

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

 

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

 

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

Изменено пользователем p_v

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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-й биквадратной секции.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

47 minutes ago, thermit said:

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

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

 

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

 

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

Изменено пользователем p_v

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

1 minute ago, arhiv6 said:

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

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

 

22 minutes ago, arhiv6 said:

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

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

 

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

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

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);

 

Изменено пользователем thermit

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

7 hours ago, thermit said:

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

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

 

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

 

 

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

 

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

Изменено пользователем p_v

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Переделал фильтры на "канонические". Принципиально ничего не поменялось. Но теперь хотя бы нет сомнений, что фильтры правильные. @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%. Я в привильную сторону думаю?

 

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

Изменено пользователем p_v

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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кГц. возможно, что применение медианного фильтра до понижения частоты лучше уберет рандомные "прострелы" и при этом меньше исказит полезный сигнал.

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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% с конца, то прострелы станут не особо критичными.

 

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Прогнал сэмплы через эту библиотеку 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

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Какая-то у вас ерунда со спектром. Наверное надо все-таки 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 фильтрами надежно обкоцать.

Изменено пользователем p_v

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

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

Гость
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...