Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: Цифровая обработка низкочастотного аналогового сигнала
Форум разработчиков электроники ELECTRONIX.ru > Цифровая обработка сигналов - ЦОС (DSP) > Алгоритмы ЦОС (DSP)
bullit
Добрый день!

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

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

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

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

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

Заранее огромное спасибо!
yanvasiij
На мой взгляд в Вашем случае наиболее удобно быстрое преобразование Фурье.
bullit
Тогда может Герцеля? Если я ничего не путаю то когда ищется конретная частота, то Герцель предпочтительней?
yanvasiij
Да, согласен. Вам же нужна только одна частота. Можно еще попытаться воспользоваться обычным цифровым полосовым фильтром, вот только частота сети и искомого сигнала слишком близко друг другу, потребуется большой порядок фильтра.
STAR_IK
Действительно, тут лучше применить Герцеля, только с окном, чтоб улучшить его АЧХ. Кстати, если нужна информация только о наличии сигнала, то можно не домножать на комплексный коэффициент, тогда алгоритм Герцела превращается в обычный БИХ фильтр 2-го порядка.
bullit
А что насчёт частоты дискретизации? Какой её выбрать?
V_G
Цитата(bullit @ Apr 11 2013, 22:39) *
А что насчёт частоты дискретизации? Какой её выбрать?

Если сетевые наводки серьезные, и нуждаются именно в цифровой фильтрации, то стоит посмотреть спектр (хотя бы через звуковую карту компьютера), и выбрать частоту, где гармоники сетевой наводки уже будут ниже требуемого уровня шумов. На эту частоту строим аналоговый антиалиасинговый ФНЧ, дискретизируем с частотой по Котельникову.
TSerg
Начните с изучения помеховой обстановки.
Все остальное - потом.
bullit
Сетевые наводки зачастую сравни сигналу. Источник сигнала оч слабый.
Данные АЦП как видите есть в файле, так что найти спектр сигнала не проблема...
Будем посмотреть!

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

Другое дело предварительно замерить помехи и сообщить о недопустимом уровне шума, чтоб не прозевать синал!
Gyga
Если частота известна то может использовать корреляцию?
АНТОН КОЗЛОВ
цифровой режекторный фильтр на 50Гц или ФНЧ на 30 Гц
bullit
ФНЧ? эт чтоб срез на 30 Гц и уже на 50 Гц порядка -40-60 дБ?
Не подскажите какой цифровой фильтр взять?
V_G
Их не берут, их рассчитывают. В матлабе.
Получится КИХ примерно 30 порядка при частоте дискретизации 200 Гц. Но с частотой дискретизации сначала все-таки стоит определиться по помеховой обстановке.
andyp
Хороший режекторный фильтр можно сделать используя comb filter и правильно подгадав частоту дисретизации, чтобы минимум комба попал как раз на частоту, которую нужно давить (ну т.е. она должна быть кратной тому, что давите). Получается дешево и сердито. Он же успешно задавит гармоники помехи.

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

PS результаты синтеза можно увидеть здесь
Lmx2315
QUOTE (bullit @ Apr 12 2013, 07:43) *
Сетевые наводки зачастую сравни сигналу. Источник сигнала оч слабый.

..если надо бороться с помехами - то используйте фильтрацию, БПФ больших порядков и прочее, если сигнал очень слабый и под шумами - то дополнительно используйте оверсэмплинг - чтобы увеличить цифровое усиление.
Чем больше добьётесь соотношения: частота дискретизации/полоса пропускания фильтра, тем больше вытянете сигнал из шумов.
tmtlib
Цитата(bullit @ Apr 11 2013, 11:47) *
Ранее в проекте закладывали 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
Цитата(bullit @ Apr 11 2013, 12:47) *
Добрый день!

Есть задача: необходимо определить наличие сигнала (синусоида) конкретной частоты (заранее известной, скажем 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 Гц.
Нажмите для просмотра прикрепленного файла
И файл маткада для быстрого старта...
Нажмите для просмотра прикрепленного файла
_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
Цитата(bullit @ Apr 22 2013, 12:38) *
А правильно ли я понимаю что строка "s = s + ((data + coeff * sprev - sprev2) - s) * Kt" это окно?
А имеет ли смысл "ставить" цифровой фильтр предварительно? Ведь сам алгоритм и есть фильтр!

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

Kt здесь это просто аналог дифференцирующей RC цепочки с характерным временем 1/(1-Kt) сэмплов.
иначе мощность посчитанная Герцелем на заданной частоте будет просто монотонно растущая функция, так как это будет мощность от самого начала и до текущего момента, данный БИХ ФВЧ помогает Герцелю "забыть" отчсёты которые были давно, позже чем 1/(1-Kt) сэмплов назад.
bullit
В приложенном файле 3 набора: сильный сигнал, чуть сильнее помехи, и слабый сигнал.
Нажмите для просмотра прикрепленного файла
_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
Потрясающий результат, огромное спасибо!
mihalevski
Цитата(bullit @ Apr 13 2013, 16:23) *
ФНЧ? эт чтоб срез на 30 Гц и уже на 50 Гц порядка -40-60 дБ?
Не подскажите какой цифровой фильтр взять?


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

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

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

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

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

Сейчас вот есть такой алгоритм:
1) провожу подсчет среднего за N точек (синий на графике) (1/N кратен (частота сэмплов/частота искомого сигнала)) - по-моему это децимация называется? поправьте если ошибаюсь.
2) отрицательные значения умножаем на -1
2) считаем среднее за интервал длительности сигнала
Получаем что-то подобное:
Нажмите для просмотра прикрепленного файла
Критерий: сигнал выше среднего 60-80 % времени - но критерий какой-то странный))
bullit
Вот кстати различия между шумом и сигналом:
Нажмите для просмотра прикрепленного файла
_pv
Цитата(bullit @ Aug 4 2013, 15:46) *
В том то и дело, что не знаю где шум, а где сигнал!
Я могу подсчитать хоть средне, хоть среднеквадратичное, но только всего сигнала, ну т.е. какого-то интервала.
Но я не знаю где шум, а где сигнал. тут получается замкнутый круг))

когда уровень вдруг становится больше чем 3 (или больше) сигма, скорее всего (99.7%) это уже не шум а сигнал, и считать среднеквадратичное значение шума надо прекратить.
Цитата(bullit @ Aug 4 2013, 15:46) *
К сожалению от алгоритма герцеля пришлось отказаться, так как ооочень много шума, пришлось ставить ПФ порядком выше.
Хотя может после ПФ герцеля поставить?

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

Цитата(bullit @ Aug 4 2013, 15:46) *
Сейчас вот есть такой алгоритм:
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
Цитата(bullit @ Aug 4 2013, 20:53) *
Что-то у меня значение RMS по выше приведенной формуле просто постоянно растёт - никак понять не могу что я делаю не так.

потому что она неправильная.
Y0 += (x - Y0) / N;
Y2 += ((x-Y0)^2 - Y2) / N
bullit
Растёт, только с большей скоростью)
Я так понимаю N это всё таки константа?
И судя по всему это нечто иное чем RMS, не? Тем более, если N не константа.
Поправьте если я ошибаюсь!
_pv
Цитата(bullit @ Aug 5 2013, 01:40) *
И судя по всему это нечто иное чем RMS, не?

N - константа времени.
и да, по ссылке правильно.
а это просто замена честного усреднения БИХ фильтром.
разница лишь в том что если считать честно не с начала до конца, а только за последние n отсчётов, то придётся хранить все эти n отсчётов в неком кольцевом буфере, а так не надо.
Corner
Для точного определения где шум, а где сигнал наиболее удобно FFT. Белый шум равномерно размазан по спектру, а детерминированный сигнал будет превышать эту величину. Причем при узком спектре полезного сигнала соотношение С/Ш в рамках одной полоски будет тем выше, чем больше точек у FFT.
bullit
Огромное спасибо за подсказку - это прекрасная идея. Только я не стал брать БПФ. А просто из отфильтрованного сигнала "вырезал" режкторным фильтром полезный сигнал - получил уровень шума. есть конечно нюансы по масштабу сигнала - но преодолимо!
И работает прекрасно - но нужно конечно подпиливать!
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Invision Power Board © 2001-2025 Invision Power Services, Inc.