Цифровая обработка низкочастотного аналогового сигнала, Полезный сигал ниже 50 Гц Опции

[tmtlib]

Ранее в проекте закладывали oversampling. Но практической пользы я пока не нашел.

видимо неправильно закладывали. Должно помогать.
Смотрю у вас частота дискретизации 10500Гц, а хотите померять 10-30Гц, вот вам и оверсэмплинг.
Осталось сделать downsampling и вытащить низкие частоты вверх.
Советую записать ваш сигнал в WAV-файл, открыть звуковой редактор Audacity, после загрузки файла выбрать режим отображения Spectrogram вместо Waveform. Затем Effect -> Change Speed -> 200%, потом ещё раз и ещё и смотрите на спектре, появляются ли ваши частоты. Появиться они дожны соответсвтенно выше, чем на 10-30Гц, в зависимости от того, насколько вы ускорили запись в Change Speed. Если что-то увидите, есть смысл делать алгоритм. Вот нетребовательный к памяти метод:

Code :
//Filtres décimateurs
// T.Rochebois
// Based on
//Traitement numérique du signal, 5eme edition, M Bellanger, Masson pp. 339-346
class Decimateur5
{
private:
float R1,R2,R3,R4,R5;
const float h0;
const float h1;
const float h3;
const float h5;
public:
Decimateur5::Decimateur5():h0(346/692.0f),h1(208/692.0f),h3(-44/692.0f),h5(9/692.0f)
{
R1=R2=R3=R4=R5=0.0f;
}
float Calc(const float x0,const float x1)
{
float h5x0=h5*x0;
float h3x0=h3*x0;
float h1x0=h1*x0;
float R6=R5+h5x0;
R5=R4+h3x0;
R4=R3+h1x0;
R3=R2+h1x0+h0*x1;
R2=R1+h3x0;
R1=h5x0;
return R6;
}
};

отсюда, есть на девять точек http://www.musicdsp.org/archive.php?classid=3#231
Легко переделывается для многокаскадности, в 2-4-8-16 и т.д. раз, место нужно только под переменные h,R каждого из каскадов. Я цеплял 8-каратную децимацию прямо на прерывание АЦП, работало очень хорошо и быстро, без всяких простоев в фоновом режиме.

Если ничего не получается, выложите WAV-файл, может появятся идеи.

[bullit]

Помеховую обстановку не оценить - система не статична.

Записать на звукавуху не получится. Сигнал оочень слабый.

Пока с сигналом играюсь в матлабе! Попробую выше предложенные алгоритмы в нём прогнать!

[vgo1]

Добрый день!

Есть задача: необходимо определить наличие сигнала (синусоида) конкретной частоты (заранее известной, скажем 25 Гц) на предмет её наличия. Диапазон от 10-30 Гц. Сигнал слабый, и при этом "перемешен" с 50 Гц (наводки/помехи). Сигнал имеет паузы: 1 секунду сигнал есть, 2 сек. нет.

Ранее в проекте закладывали oversampling. Но практической пользы я пока не нашел.

Может кто нить подсказать, что даёт этот метод, и как им пользоваться.

Исходные данные:
1) необходимо детектировать (сказать что сигнал такой то частоты присутствует!) сигнал частотой скажем 25Гц
2) частота дискретизации порядка 10500 Гц.
3) устройство микропроцессорное, т.е. все фильтры должны быть подсилу мк.
Пример сигнала на рисунке:
[attachment=76335:first.JPG]
А может сигнал быть по уровню с помехой, а то и меньше!

Что необходимо сделать?

Заранее огромное спасибо!

Я с помощя корреляции идентицировал частоты АОН, DTMF и ответ станции. Атмеловский мк ATmega8 работаящий с кварцем 16МГц успевал. Чтобы ускорить расчеты пользовался таблией заранее расчитанных значений синусов(256 значений).
Но от аналогового (или на переклячаемых конденсаторах)фильтра на входе АЦП вам наверное не уйти. После фильтрации можно взять частоту дискретизации намного раз меньше(зависит от частоты среза) и ваш oversampling станет ненужным.

[bullit]

Ну аналоговый фильтр стоять будет, да только частоту его среза буду отодвигать подальше от 30, ну и конечно 50 Гц. Так как очень широкий диапазон температур, в том числе с "быстрой" сменой температурного диапазона (скажем улица - машина/здание). Так что он больше от высокочастотных помех будет, делать сложный аналоговый при широкой цифровой возможности...зачем?
Тут еще момент такой! Спектр сигнала широкий(21-25 Гц), это не пик скажем как у 50 Гц (см. рисунок), а более широкое распределение (оч похоже на Гауссовское распределение). Поэтому разница "сигнал - шум" не такая большая (это уровне сигнала как в топике). Может стоить этот момент как-то учесть? Хорошо бы площадь под ней подсчитать... но стоит ли? да и как?


ЗЫ а дайте почить про корреляцию. Чёт ничего путного не нашел!

[tmtlib]

А в каком формате данные? Если есть возможность выложите кусочек

[bullit]

Могу дать в текстовый файле в архиве!
Частота оцифровки 10416 Гц.
 test26.rar ( 124.94 килобайт ) Кол-во скачиваний: 75

И файл маткада для быстрого старта...
 openfileandgoerthel.rar ( 705 байт ) Кол-во скачиваний: 60

[_pv]

так чем Герцель-то не угодил?

test.lua:

Код

file = io.open("test26.txt")
f = 22
fs = 10416
coeff = 2 * math.cos(2 * 3.1415926535 * f / fs)
s = 0
sprev = 0
sprev2 = 0
Kt = 0.995

for data in file:lines("*n") do
  s = s + ((data + coeff * sprev - sprev2) - s) * Kt
  sprev2 = sprev
  sprev = s
  power = sprev2 * sprev2 + sprev * sprev - coeff * sprev * sprev2
  print(power * 1e-9)                
end
file:close()

lua52 test.lua > out
gnuplot -e "plot 'out' with lines, 'test26.txt' with lines" -

[bullit]

Неплохие результаты!
Особенно поигравшись с коэффициентом.
В принципе Герцеля я и предполагал поставить изначально!

А правильно ли я понимаю что строка "s = s + ((data + coeff * sprev - sprev2) - s) * Kt" это окно?

А имеет ли смысл "ставить" цифровой фильтр предварительно? Ведь сам алгоритм и есть фильтр!

PS с ЦОСом я на вы, вот белые дыры тяхонько закрываю...

[_pv]

А правильно ли я понимаю что строка "s = s + ((data + coeff * sprev - sprev2) - s) * Kt" это окно?
А имеет ли смысл "ставить" цифровой фильтр предварительно? Ведь сам алгоритм и есть фильтр!

это что-то вроде окна, сделанное простым БИХ фильтром.
то есть по идее, наверное, надо было просто посчитать производную от power и опять отфильтровать. что примерно данный фильтр и делает.
предварительно фильтр не особо поможет, т.к. герцель это просто преобразование Фурье для заданной частоты т.е. лучше уже не будет.
разве что для децимации, чтобы не тащить такой поток данных, но так как данный Герцель будет делаться быстрее чем нормальный ФНЧ для децимации, то пожалуй это бесполезно, а то и хуже сделает.

[bullit]

К сожалению при слабом сигнале (уровень сигнала = шуму) необходимо менять коэффициент Kt. Иначе либо хорошо видно при мощном сигнале, либо при слабом... менять в процессе имхо не дело!
И как бы итоговый массив после Герцеля обработать? т.е. проверка на то что сигнал длится столько то периодов и т.д.?

[_pv]

образцы сигналов в студию.

Kt здесь это просто аналог дифференцирующей RC цепочки с характерным временем 1/(1-Kt) сэмплов.
иначе мощность посчитанная Герцелем на заданной частоте будет просто монотонно растущая функция, так как это будет мощность от самого начала и до текущего момента, данный БИХ ФВЧ помогает Герцелю "забыть" отчсёты которые были давно, позже чем 1/(1-Kt) сэмплов назад.

[bullit]

В приложенном файле 3 набора: сильный сигнал, чуть сильнее помехи, и слабый сигнал.
 tests.rar ( 690.05 килобайт ) Кол-во скачиваний: 80

[_pv]

Код

power_der = 0
power_filt = 0
power_prev = 0
Kt = 1 - 1e-3
Kt2 = 1e-3
Kt3 = 1-1e-5
peak = 0
output = 0
for data in file:lines("*n") do
  s = data + coeff * sprev - sprev2
  sprev2 = sprev
  sprev = s
  sprev =  sprev * Kt
  sprev2 =  sprev2 * Kt
  power = sprev2 * sprev2 + sprev * sprev - coeff * sprev * sprev2
  power_der = (power - power_prev)
  power_prev = power
  if (power_der < 0) then power_der = 0 end
  power_filt = power_filt + (power_der - power_filt) * Kt2
  if power_filt > peak then peak = power_filt end
  peak = peak * Kt3
  if power_filt > peak * 0.5 then output = 1 else output = 0 end
end

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

[bullit]

Потрясающий результат, огромное спасибо!

[bullit]

Тесты
 tests.rar ( 1.37 мегабайт ) Кол-во скачиваний: 74