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

[bullit]

Добрый день!

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

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

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

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


А может сигнал быть по уровню с помехой, а то и меньше!

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

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

[yanvasiij]

На мой взгляд в Вашем случае наиболее удобно быстрое преобразование Фурье.

[bullit]

Тогда может Герцеля? Если я ничего не путаю то когда ищется конретная частота, то Герцель предпочтительней?

[yanvasiij]

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

[STAR_IK]

Действительно, тут лучше применить Герцеля, только с окном, чтоб улучшить его АЧХ. Кстати, если нужна информация только о наличии сигнала, то можно не домножать на комплексный коэффициент, тогда алгоритм Герцела превращается в обычный БИХ фильтр 2-го порядка.

[bullit]

А что насчёт частоты дискретизации? Какой её выбрать?

[V_G]

А что насчёт частоты дискретизации? Какой её выбрать?

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

[Guest_TSerg_*]

Начните с изучения помеховой обстановки.
Все остальное - потом.

[bullit]

Сетевые наводки зачастую сравни сигналу. Источник сигнала оч слабый.
Данные АЦП как видите есть в файле, так что найти спектр сигнала не проблема...
Будем посмотреть!

Изучить помеховую обстановку нереально: система мобильная и никто не скажет что будет на новом месте. Возможно что даже всё чисто будет (чистополе), а возможно и мощная энергетическая установка.

Другое дело предварительно замерить помехи и сообщить о недопустимом уровне шума, чтоб не прозевать синал!

[Gyga]

Если частота известна то может использовать корреляцию?

[АНТОН КОЗЛОВ]

цифровой режекторный фильтр на 50Гц или ФНЧ на 30 Гц

[bullit]

ФНЧ? эт чтоб срез на 30 Гц и уже на 50 Гц порядка -40-60 дБ?
Не подскажите какой цифровой фильтр взять?

[V_G]

Их не берут, их рассчитывают. В матлабе.
Получится КИХ примерно 30 порядка при частоте дискретизации 200 Гц. Но с частотой дискретизации сначала все-таки стоит определиться по помеховой обстановке.

[andyp]

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

Вобщем, ключевые слова для гугла - comb notch filter

PS результаты синтеза можно увидеть здесь

Сообщение отредактировал andyp - Apr 13 2013, 22:13

[Lmx2315]

Сетевые наводки зачастую сравни сигналу. Источник сигнала оч слабый.

..если надо бороться с помехами - то используйте фильтрацию, БПФ больших порядков и прочее, если сигнал очень слабый и под шумами - то дополнительно используйте оверсэмплинг - чтобы увеличить цифровое усиление.
Чем больше добьётесь соотношения: частота дискретизации/полоса пропускания фильтра, тем больше вытянете сигнал из шумов.

[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) устройство микропроцессорное, т.е. все фильтры должны быть подсилу мк.
Пример сигнала на рисунке:


А может сигнал быть по уровню с помехой, а то и меньше!

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

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

Я с помощя корреляции идентицировал частоты АОН, 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

[mihalevski]

ФНЧ? эт чтоб срез на 30 Гц и уже на 50 Гц порядка -40-60 дБ?
Не подскажите какой цифровой фильтр взять?

Если нет необходимости в определении точного значения амплитуды сигнала (частота точно неизвестна) можно применить очень простой КИХ фильтр (на основе циклотомических полиномов с единичными коэффициентами) который бедет давить сетевую помеху и ее гармоники на указанные величины. Где-то ранее расчет такого фильтра я здесь приводил. Частота дискретизации 600 Гц.

[bullit]

Возникла такая задача:
Как лучше всего на уровне шума, уровень которого не известен, ловить приход сигнала?
Как подсчитывать уровень шума? при том что о может немного во времени меееедленно меняться!

[_pv]

что значит не известен?
среднеквадратичное значение всегда можно посчитать

[bullit]

В том то и дело, что не знаю где шум, а где сигнал!
Я могу подсчитать хоть средне, хоть среднеквадратичное, но только всего сигнала, ну т.е. какого-то интервала.
Но я не знаю где шум, а где сигнал. тут получается замкнутый круг))

Объясню чуть подробней: вот есть сигнал, как на картинке в вашем посте 28. Визуально я могу выделить сигнал и шум, а как программно - не знаю как подобрать нормальный алгоритм.

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

Да еще чтоб при отсутствии сигнала, не ловить на уровне шума ложные срабатывания.

Мне кажется нужно использовать то, что отличает шум от сигнала: 1) это уровень 2) гармоничность 3) длительность

Сейчас вот есть такой алгоритм:
1) провожу подсчет среднего за N точек (синий на графике) (1/N кратен (частота сэмплов/частота искомого сигнала)) - по-моему это децимация называется? поправьте если ошибаюсь.
2) отрицательные значения умножаем на -1
2) считаем среднее за интервал длительности сигнала
Получаем что-то подобное:


Критерий: сигнал выше среднего 60-80 % времени - но критерий какой-то странный))

[bullit]

Вот кстати различия между шумом и сигналом:

[_pv]

В том то и дело, что не знаю где шум, а где сигнал!
Я могу подсчитать хоть средне, хоть среднеквадратичное, но только всего сигнала, ну т.е. какого-то интервала.
Но я не знаю где шум, а где сигнал. тут получается замкнутый круг))

когда уровень вдруг становится больше чем 3 (или больше) сигма, скорее всего (99.7%) это уже не шум а сигнал, и считать среднеквадратичное значение шума надо прекратить.

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

что-то мне подсказывает что интеграл Фурье на заданной частоте это оптимальный полосовой фильтр,

Сейчас вот есть такой алгоритм:
1) провожу подсчет среднего за N точек (синий на графике) (1/N кратен (частота сэмплов/частота искомого сигнала)) - по-моему это децимация называется? поправьте если ошибаюсь.
2) отрицательные значения умножаем на -1
2) считаем среднее за интервал длительности сигнала
Получаем что-то подобное:


Критерий: сигнал выше среднего 60-80 % времени - но критерий какой-то странный))

если уровень меньше чем k * sigma (k=3 например) считаем среднеквадратичное значение дальше, если нет, то это сигнал и не обновляем среднеквадратичное значение.

причём чтобы не хранить кольцевой буфер для честного вычисления среднего и среднеквадратичного значений, можно считать немного проще:
Y0 += (x - Y0) / N;
Y2 += (x - Y0)^2 / N;
stDev = sqrt(Y2);

[bullit]

Если Вы посмотрите последний рисунок, то увидите, что там отношение более чем 3.
Может лучше считать rms на интервале половины или кратной искомой частоты?

За формулы rms спасибо!

Ну кстати!
Сделал расчёт на интервалах пол периода искомого сигнала в течении секунды. И получился достаточно не плохой результат. Даже трехкратного запаса не надо. Единственное как первоначально понять что мы не на сигнале, а в паузе. Либо первоначально смотреть на весь период следования сигнала (2 сек = 1 сек сигнал + 1 сек пауза) и уже исходя из этого смотреть. Ну вроде очень даже не плохой результат!
Спасибо. Будем дальше испытывать!

Что-то у меня значение RMS по выше приведенной формуле просто постоянно растёт - никак понять не могу что я делаю не так.

[_pv]

Что-то у меня значение RMS по выше приведенной формуле просто постоянно растёт - никак понять не могу что я делаю не так.

потому что она неправильная.
Y0 += (x - Y0) / N;
Y2 += ((x-Y0)^2 - Y2) / N

[bullit]

Растёт, только с большей скоростью)
Я так понимаю N это всё таки константа?
И судя по всему это нечто иное чем RMS, не? Тем более, если N не константа.
Поправьте если я ошибаюсь!

[_pv]

И судя по всему это нечто иное чем RMS, не?

N - константа времени.
и да, по ссылке правильно.
а это просто замена честного усреднения БИХ фильтром.
разница лишь в том что если считать честно не с начала до конца, а только за последние n отсчётов, то придётся хранить все эти n отсчётов в неком кольцевом буфере, а так не надо.

[Corner]

Для точного определения где шум, а где сигнал наиболее удобно FFT. Белый шум равномерно размазан по спектру, а детерминированный сигнал будет превышать эту величину. Причем при узком спектре полезного сигнала соотношение С/Ш в рамках одной полоски будет тем выше, чем больше точек у FFT.

[bullit]

Огромное спасибо за подсказку - это прекрасная идея. Только я не стал брать БПФ. А просто из отфильтрованного сигнала "вырезал" режкторным фильтром полезный сигнал - получил уровень шума. есть конечно нюансы по масштабу сигнала - но преодолимо!
И работает прекрасно - но нужно конечно подпиливать!