Цифровая обработка аналогового сигнала

[Santik]

...Лично мне картинка от "использовать коррелятор как узкополосный фильтр" очень понравилась. Можно ли её "внедрить" в контроллер?...

Да, конечно, это возможно.
Но надо чётко понимать, что связываться с корреляцией есть смысл, если у Вас есть возможность частотно модулировать исходный сигнал. Если такой возможности нет, то гораздо проще применить цифровой узкополосный фильтр.
Коррелятор реального времени требует от контроллера очень высокой производительности. За время равное шагу дискретизации надо выполнить М операций вида Ах+В, где М - число выборок эталонного сигнала (с которым коррелируется принимаемый сигнал), х - текущая выборка принимаемого сигнала.
В STM32F4 операция вида Ах+В занимает 3 машинных цикла.
В принципе, сделать коррелятор реального времени на таком контроллере, учитывая, что тактовая частота 168 МГц, вполне реально.
Применительно к Вашей задаче - придётся значительно уменьшить частоту дискретизации, что потребует применения качественного аналогового анти-аляйсинг фильтра на входе АЦП.

[Santik]

REC_06.06.2013 23-22-52_signal.txt

Коррелограмма с синусоидальным сигналом 13 Гц 0.5 сек:

[bullit]

Эх не хотелось бы нагружать проц, ибо на нём еще других задачи висит не мало!
А вот по поводу "то гораздо проще применить цифровой узкополосный фильтр" - разве Герцель не делает тоже самое? я считаю его для одной частоты. Хотя результат мне не очень...
Вы можете показать код такого фильтра? и главное как его перестроить на другую частоту 10-30 Гц.
Есть надежда на более чёткое определение границ пучка сигнала с таким фильтром.

[bullit]

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

Такой возможности нет воопче.

А вот с цифровым фильтром мне разобраться не удалось...проект делаю в фоне...

[Santik]

... А вот с цифровым фильтром мне разобраться не удалось...

Это плохо!

...проект делаю в фоне...

Это хорошо!

Чтобы разобраться с цифровыми фильтрами, начнём издалека...

Режекторный фильтр на 50 Гц:

y(k) = G·[b0·x(k) +b1·x(k-1)+b2·x(k-2)] – a1·y(k-1) – a2·y(k-2).

х(k) - исходные выборки
y(k) - выборки отфильтрованного сигнала
R=1.01 ; R2=R*R
G=0.99

Для Fd=10416 Гц (FN=5208) и частоте режекции Fs=50 Гц можно вычислить коэффициенты b1,a1, a2
zn = cos (3.1415926*Fs/Fn)+j·sin(3.1415926*Fs/Fn)
zp=zn*R
b0=1
b1 = -2·Re{zn} =-1.999090369781252
b2 =1

a1 = - (2·Re{zp})/R2 =-1.9786996480094832296
a2 = 1/R2 = 0.98.
Окончательно:
Y(k)=0.99*(X(k) -1.999090369781252 *X(k-1) + X(k-2)) +1.9786996480094832296* Y(k-1)-0.98* Y(k-2)

Таким образом мы 50 Гц вырезали.

Аналогично можно построить фильтр пропускающий только сигнал 13 Гц (или 22 Гц - только коэффициенты будут другие)

[bullit]

Спасибо большое за фильтр!

Буду мучать в матлабе, слава богу Матлаб освил достаточно))

... А вот с цифровым фильтром мне разобраться не удалось...

Это плохо!

...проект делаю в фоне...

Это хорошо!

Честно говоря жаль что с ЦОСом имею мало опыта, щас вот тяхонько "закрываю белые пятна" знаниях. матлаб здорово помогает закрепить теорию практикой.

[Santik]

Забыл выражение для коэффициента G:
G=(1.0+(1.0+2.0*Real(Zp))/(R**2))/(2.0+2.0*Real(Zn))
И в значении кажется ошибся...

[Santik]

Селекторный фильтр.
Передаточная функция режекторного фильтра:



Передаточная функция селекторного фильтра:

или
где
Уравнение фильтра:
Y(k)=(1.0-G)*X(k) +(a1-G*b1)*X(k-1) +(a2-G)* X(k-2) -a1*Y(k-1)-a2* Y(k-2)
На рисунке данные из файла REC_06.06.2013 23-26-58_signal.txt после режекторного фильтра 50 Гц и селекторного фильтра 13 Гц

Фрагмент программы (Фортран):

Код: Выделить всё

!Комментарий
! Режекторный фильтр 50 Гц. Алгоритм фильтра: y(k) = G[x(k) + b1 x(k-1) + x(k-2)] - a1 y(k-1) -a2 y(k-2).
! X(k)  - выборки исходного сигнала

R=1.01 ! Параметр отвечающий за "устойчивость" фильтра. Должен быть > 1
Zn= CMPLX(cos(Pi*50.0/5208.), sin(Pi*50.0/5208.)) !Мнимую часть Zn можно не вычислять
Zp=R*Zn
G=(1.0+(1.0+2.0*Real(Zp))/(R**2))/(2.0+2.0*Real(Zn))
a1=-2.0*Real(Zp)/R**2
a2=1.0/(R**2)
b1=-2.0*Real(Zn)
b2=1.0

X1(1)=0.0
X1(2)=0.0
do k=3, N2 ! N2-общее число выборок сигнала
X1(k)=G*(X(k)  +b1*X(k-1) +b2*X(k-2)) -a1* X1(k-1)-a2* X1(k-2)
end do

do k=1, N2
X(k)=X1(k) 
end do

! Селекторный фильтр 13 Гц 
R=1.001    !параметр уменьшен для большего подавления частот не равных 13 Гц
Zn= CMPLX(cos(Pi*13.0/5208.), sin(Pi*13.0/5208.))
Zp=R*Zn
G=(1.0+(1.0+2.0*Real(Zp))/(R**2))/(2.0+2.0*Real(Zn))
a1=-2.0*Real(Zp)/R**2
a2=1.0/(R**2)
b1=-2.0*Real(Zn)
b2=1.0
  
X1(1)=0.0
X1(2)=0.0
do k=3, N2
X1(k)=(1.0-G)*X(k) +(a1-G*b1)*X(k-1) +(a2-G)* X(k-2)  -a1*X1(k-1)-a2* X1(k-2)
end do

do k=1, N2
X(k)=X1(k) 
end do

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

На картинках исходный сигнал -> режекторный фильтр -> 3x и 5х селекторный

Ну вот где-то так...

[bullit]

Будем пробывать, Спасибо!

PS Фортран? я думал этот язык уже исчез... мой первый язык программирования, изучал по советской книжке в красной обложке! гугль тока чёт не выдает по поиску...

[Santik]

Матлабом никогда не пользовался... Но думаю, что переписать с Фортрана труда не составит.
Тем более что: (Материал из Википедии )
"MATLAB как язык программирования был разработан Кливом Моулером (англ. Cleve Moler) в конце 1970-х годов... Целью разработки служила задача дать студентам факультета возможность использования программных библиотек Linpack и EISPACK без необходимости изучения Фортрана... До сих пор в Интернете можно найти версию 1982 года, написанную на Фортране, распространяемую с открытым исходным кодом.