Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

BaRaGoZ
Сообщения: 4
Зарегистрирован: 07 апр 2014, 18:32

Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

Сообщение BaRaGoZ »

Есть задача: МК обрабатывает данные с акселерометра, сбор данных происходит каждые 2мс, отсчеты собираются в течении 1с, те 500(можно 512, если использовать БПФ) за 1 секунду. Далее необходимо найти уровни сигналов с частотами от 2-5 Гц, после чего принимается решение о превышении уровня сигналов допустимого значения. Делаю так: набираю массив данных в течение 1с, затем отсчеты скидываются по КОМ-порту на ПК в Матлаб. Затем полученные отсчеты прогоняю через FFT и на выходе должен увидеть спектр частот. Все прекрасно работает если частота воздействия на акселерометр от 100 Гц и выше(более 250Гц не проверял). Но если воздействую на акселерометр с частотой примерно 2-6 Гц, то в спектре эти частоты практически «не выделяются», вот картинка(1 график-исходный сигнал,2график-после FFT ф-ии): http://sb.uploads.ru/qsm8r.png
Причем воздействие видно даже визуально+амплитуда сигнала воздейтсвия много выше уровня шума! Попробовал просто «вырезать» шум, те в коде программы просто делаю (если noise_min<=сигнал<=noise_max, то сигнал = 0), тогда ситуация улучшилась, и сигналы нужных частот стали видны в спектре, сделал я это как то интуитивно, из-за того что уровень амлитуды определенной частоты - это энергия сигнала данной частоты на 1ом графике, а там визуально видно, что воздействие хоть и отчетливое, но уж очень кратковременное и по сравнению с энергией шума, его энергия мала. Поэтому не знаю насколько такое решение правомерно и как оно будет работать в общем случае, те в разных условиях и не понятно какие границы шума обрезать и тд. В связи с этим появились вопросы: 1)Возможно ли при таком исходном сигнале выделить низкие частоты, если да то как? 2)Надо ли обрабатывать дополнительно исходный сигнал, до преобразований Фурье, если да, то как грамотно это сделать?

BaRaGoZ
Сообщения: 4
Зарегистрирован: 07 апр 2014, 18:32

Re: Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

Сообщение BaRaGoZ »

Уточню задачу и алгоритм: Надо сделать устройтсво, которое улавливает шаги на расстоянии N метров на основе акселерометра. Был предложен алгоритм: собираем информацию от дачика в течение 6 секунд(на секунду 500 отсчетов), задем делаем для каждой секунды ДПФ(БПФ), находим среднюю ампл-ду каждой секунды в полосе от 2 до 5 Герц, далее находим коэффициент К=Аср6/( (Аср1+Аср2+Аср3+Аср4+Аср5)/5 ), если К превышает порог, значит шаги есть. Затем делаем А1=А2;А2=А3;А3=А4;А4=А5;А5=А6; а А6 считаем снова( т.е. собираем 500 отсчетов за секунду и тд.). Алгоритм и математика крутятся в МК, просто скидываю все отсчеты и результаты по интерфейсу в Windows(написал программку на WinApi и использую также Матлаб) для того, чтобы проверить правильность работы алгоритма МК(сравниваю что получаем на МК и в матлабе например)+для того,чтобы все результаты видеть наглядно. Т.к. алгоритм изначально был предложен, я его и придерживаюсь. По поводу RC фильтра аппаратного - этот вариант не подходит, дабы в перспективе могут понадобиться и другие частоты, например ловить буду не только шаги и тд. А вот по поводу программного фильтра, это интересней. Я сегодня как то интуитивно понял, что пики сигналов приходящих надо сгладить,сделал так - разбил 500 отсчетов на N отрезков(пробовал 100 и 50 отрезков), каждый отрезок усреднял, на выходе получается набор «столбиков», потом все это хозяйство подсовывал в Фурье, результат намного лучше!Теперь осталось осознать какой фильтр выбрать, и как его правильно рассчитать,просто под рукой нету интернета, поэтому есть некоторая ограниченность при экспериментах. Поэтому теперь такие вопросы появились: 1)Какой фильтр лучше выбрать в данной ситуации? 2)Повлияет ли фильтр на «точность» измерений, просто на больших расстояниях всплески сигнала будут не такие большие... ? И если можно пример расчета фильтра или советы по его расчетам! Спасибо, заранее всем благодарен!

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

Сообщение Бахурин Сергей »

Ну то что вы разбили на куски и усредняли это не совсем корректно. Лучше если будете сдвигать по одному отсчету и пересчитывать выход. Т.е. первый шаг


Второй шаг


и так далее сдвигать по одному отсчету. Ваши столбики еще более сгладятся и должно стать еще лучше.

BaRaGoZ
Сообщения: 4
Зарегистрирован: 07 апр 2014, 18:32

Re: Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

Сообщение BaRaGoZ »

Спасибо,попробую, только формула немного не ясна,мне кажется
1) в сумме нижний предел должен увеличиваться. а верхний уменьшаться, а то последний последний шаг получиться: от к=N-1 до 2N, а это как то не понятно...
2) надо нормировать y(i), например делить на количество членов суммы.
Или я не прав?
Пробовал еще фильтр НЧ, который в матлабе идет функцией fir1(), а затем результат закидывать в Фурье, частоты действительно низкие отображаются в спектре, но как то не точно, т.е при 4х всплесках максимальная частота все равно 11ая или 17ая, а нужные 2-5Гц высокие, но не максимальные,да и как то стабильности нет(график строить ся 1 раз в секунду) при примерно схожих воздействиях, т.е. стук происходит примерно 4 раза за секунду, а графики меняются, причем меняются доминирующие частоты... пробовал разные параметры ф-ии fir1 (M и FC-частоту среза), также пробовал сам считать весовые коэффициенты, результат один, подозрение что я что то не понимаю, или неправильно подбираю параметры М и FC... Как их правильно рассчитать?

Еще забыл спросить, значения с акселерометра бывают и отрицательные, частотный фильтр никак "не страдает" от этого?

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

Сообщение Бахурин Сергей »

1) разумеется вы не должны допускать выхода индекса за пределы
2)Да надо нормировать
3)Как использовали фильтр который посчитался фукцией fir1?
Фильтр не страдает от отрицательных значений

BaRaGoZ
Сообщения: 4
Зарегистрирован: 07 апр 2014, 18:32

Re: Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

Сообщение BaRaGoZ »

1)пробовал предложенный вами вариант(только смещал лишь нижний предел суммы и усреднял), график почти в линию превратился - больно гладкий+ последний элемент всегда получается максимальный(т.к. он не усредняется, просто равен сам себе)
2)если коэффициенты считал сам то делал так(это пример из книги Стивена Смита - "Цифровая обработка сигналов", только тут язык бэйсик наверное, я переделывал в матлабе, ну это сути не меняет):
100 'НИЗКОЧАСТОТНЫЙ ОКОННЫЙ ФИЛЬТР
110 'Данная программа предназначена для обработки 500 отсчётов
130 '
140 DIM X[500] 'X[ ] – массив отсчётов входного сигнала
150 DIM Y[500] 'Y[ ] – массив отсчётов выходного сигнала
160 DIM H[100] 'H[ ] – массив весовых коэффициентов НЧ-фильтра
170 '
180 PI = 3.14159265
190 FC = .012 'Определение частоты среза (в диапазоне 0…0.5)
200 M% = 65 'Определение порядка фильтра (101-й порядок)
210 '
230 '
240 ' 'Вычисление весовых коэффициентов НЧ-фильтра (16.4)
250 FOR I% = 0 TO M
260 IF (I%-M%/2) = 0 THEN H[I%] = 2*PI*FC
270 IF (I%-M%/2) <> 0 THEN H[I%] = SIN(2*PI*FC * (I%-M%/2)) / (I%-M%/2)
280 H[I%] = H[I%] * (0.54 - 0.46*COS(2*PI*I%/M%) )
290 NEXT I%
300 '
310 SUM = 0 'Нормирование коэффициентов НЧ-фильтра для получения
320 FOR I% = 0 TO M 'единичного коэффициента усиления на нулевой частоте
330 SUM = SUM + H[I%]
340 NEXT I%
350 '
360 FOR I% = 0 TO M
370 H[I%] = H[I%] / SUM
380 NEXT I%
390 '
400 FOR J% = M TO 499 'Операция свёртки с входным сигналом
410 Y[J%] = 0
420 FOR I% = 0 TO M
430 Y[J%] = Y[J%] + X[J%-I%] * H[I%]
440 NEXT I%
450 NEXT J%
460 '
470 END
Если коэффициенты считал через fir1(), то делал только свертку. Пробовал еще:
100 'ОДНОРОДНЫЙ НЕРЕКУРСИВНЫЙ ФИЛЬТР (ФИЛЬТР СКОЛЬЗЯЩЕГО СРЕДНЕГО)
110 'Данная программа предназначена для обработки 500 отсчётов
120 'входного сигнала однородным фильтром 101-го порядка.
130 '
140 DIM X[500] 'X[ ] – массив отсчётов входного сигнала
150 DIM Y[500] 'Y[ ] – массив отсчётов выходного сигнала
160 '

180 '
190 FOR I% = 50 TO 449 'Цикл по числу отсчётов выходного сигнала
200 Y[I%] = 0 'Предварительное обнуление
210 FOR J% = -50 TO 50 'Цикл для вычисления очередного отсчёта
220 Y[I%] = Y[I%] + X[I%+J%]
230 NEXT J%
240 Y[I%] = Y[I%]/101 'Завершение усреднения – операция деления
250 NEXT I%
260 '
270 END
Результаты похожи, причем и графики отфильтрованного сигнала похожи на синусоиду, нежели на некий шум, что радует, 1 фильтр поглаже немножко. Еще не понял такой момент, при использовании фильтров некоторые значения выходного сигнала = 0(тк забиваем выходной массив в цикле с значения индекса M). Насколько это правильно, и в чем смысл таких 0? При в Фурье скидываю значения выходного массива , но без 0ей., тогда кстати спектр становится похож на спектр без фильтрации. Вообщем худо-бедно получается выделить например усредненный сигнал 1-5Гц, т.е К = Аср(средняя сумма ампл. 1-5Гц)/ШУМсредний при воздействии достигает от 2-10. Хотелось бы только максимально возможно увеличить это К, убрав шум и не нужные частоты.Для этого, как мне кажется, надо правильно выбрать фильтр и рассчитать коэффициенты М и Fср, пока пользуюсь для этого опцией fdatool (М получается 65), что можете сказать по этому поводу?
Еще есть интересная статья, на тему улавливания шагов, только там сейсмограф, но смысл один: http://sem-proceedings.com/26i/sem.org- ... raffic.pdf
Там они фильтруют, затем берут модуль сигнала, затем пропускают все это через блок "Decimate", а затем уже Фурье... Вот не могу понять зачем брать модуль и что это за блок такой "Decimate"???
Заранее спасибо!
P.S.: к сожалению не смог сегодня выложить графики, в будущем постараюсь!

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: Выделение низких частот(2-6Гц) сигнала при помощи ДПФ

Сообщение Бахурин Сергей »

BaRaGoZ писал(а):Еще не понял такой момент, при использовании фильтров некоторые значения выходного сигнала = 0(тк забиваем выходной массив в цикле с значения индекса M). Насколько это правильно, и в чем смысл таких 0? При в Фурье скидываю значения выходного массива , но без 0ей.,
Да это нормально поскольку фильтр имеет задержку выхода относительно входа. Просто надо не брать в расчет начальные нули.
BaRaGoZ писал(а): Для этого, как мне кажется, надо правильно выбрать фильтр и рассчитать коэффициенты М и Fср, пока пользуюсь для этого опцией fdatool (М получается 65), что можете сказать по этому поводу?
Ну 65 так 65. Значит вы задали такие требования на фильтр что он должен иметь 65 коэффициентов чтобы ваши требования обеспечить. Если вы сможете его реализовать используя ваши ресурсы то все ок.
BaRaGoZ писал(а): Там они фильтруют, затем берут модуль сигнала, затем пропускают все это через блок "Decimate", а затем уже Фурье... Вот не могу понять зачем брать модуль и что это за блок такой "Decimate"???
Decimate значит прореживание. Они оставляют от сигнала с частотой дискретизации 1000 Гц каждый 30 отсчет только в результате понижают частоту дискретизации в 30 раз до 33.3 Гц. Теперь на более низкой частоте дискретизации легче вести обработку. Например ваш фильтр будет гораздо лучше работать.

Ответить