Собственно возникла вот какая задача. Имеем видеосигнал АЦП с частотой 66.(6) МГц, но вообще частота может программно меняться, от 25 до 80 МГц. Но некоторые приемники выдают сигнал с замешиванием в него частоты 50 Гц. И вот требуется эту помеху из сигнала убрать, однако DC нужно оставить. Обработка делается на ПЛИС и основная частота обработки в ней - 100 МГц.

Пока что в голову приходит такой способ.
1) ФНЧ, оставить полосу порядка 200-300 Гц.
2) децимация сигнала на частоту работы полосового фильтра.
3) Полосовой фильтр вокруг частоты 50 Гц.
4) ресемплинг на исходную частоту оцифровки
5) Вычитание из входного сигнала выделенной составляющей

Однако вот попробовал синтезировать какие-то фильтры в MATLAB, смущает, что на шагах 1-3 составляющая приобретает сдвиг по фазе, и становится непригодна для простого вычитания из исходного сигнала. Причем если для одной частоты еще можно как-то подравнять ФЧХ на 0 именно для частоты фильтрации, то для полосы 50-60 Гц фаза довольно сильно изменяется в полосе пропускания.

Пока задачка стоит подавить составляющую на 50 Гц где-нибудь на 35-45 дБ. На данный момент амплитуда помехи оценена в 10% динамического диапазона АЦП. При реализации фильтра на ПЛИС вполне можно успеть делать КИХ на 128-256 порядков на одном блоке памяти. Но вот все же хотелось бы давить не только одну частоту в 50 Гц, но диапазон от 50 до 60. Возможно ли это?

А теперь представьте себе что по телевизору показывают зебру крупным планом. Причём на всей левой половине экрана полоска чёрная, а на правой - белая. А потом камера перемещается к хвосту и теперь верхняя половина экрана чёрная, а нижняя белая.
И как тут отличить полезный сигнал от сетевой наводки?
А если зебра испугается и побежит?

С помехой лучше бороться там, где она возникает. Пинайте конструкторов, чтобы нормальное заземление и экранировку обеспечили.
А в сетевой наводке кроме 50 Гц есть еще и 100, и 150 и т.д. - все не отфильтруешь.

зебры на куда более высоких частотах водятся.
50Гц, это полэкрана и больше.
плюс наверное можно и фазу 50Гц отловить и выпилить куда более качественно нежели просто полосовым фильтром на 50Гц.

Ну вообще в моем случае видеосигнал - это сигнал с приемника РЛС, и он зебру никогда не показывает, особенно бегающую крупным крупом. И проблема в том, что приходится сопрягаться с множеством разных РЛС, а помеха в 50 Гц есть не у всех, а лишь у некоторых. Пинать чужих конструкторов дело неблагодарное, весь предыдущий опыт показывает, что с чужой конторой (и даже несколькими сразу) где сядешь, там и слезешь.

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

Особенности сигнала таковы, что ему нельзя удалять DC, иначе потом трудно восстанавливать уровень сигнала для выделения целей. Кроме этого, после приведения DC к нулю еще производится временная автоматическая регулировка усиления (ВАРУ), и поэтому наличие такой вот помехи перед ВАРУ приводит к ее усложнению после ВАРУ. И еще проблема в том, что разные РЛС имеют видеосигналы разной полярности и мы обычно выделяем сигнал только одной полярности, заменяя значения от другой полярности на нули. Ну и вот если сигнал начинает из-за помехи "проваливаться" в другую полярность, то приходится либо ему ложный сдвиг выдумывать, чтобы вытащить целиком в свою полярность, либо он просто обнуляется и портит выделение целей. И то, и другое плохо, лучше всего убрать 50 Гц до того, как отсекать нужную полярность. Фактически сейчас и так есть коррекция ошибки смещения АЦП (конвейерного типа), но она статическая, и если к ней добавить медленную динамическую составляющую, то в общем почти ничего сильно менять даже не надо. Но вот нужно выделить эту составляющую из видеосигнала.

А какая у вас частота следования зондирующих сигналов? Откуда взялась проблема с 50гц?

Частота в зависимости от станции и ее режима может быть от 750 Гц до 3 КГц. Обычно у наших станций частота в зависимости от шкалы меняется дискретно: 750 -> 1500 -> 3000 КГц. В любом случае на период частоты 50 Гц, приходится несколько десятков зондирований.

Проблема 50Гц взялась, я думаю, от особенностей блоков питания собственно станции и нашей обработки. И особенностей организации цепей, уравнивающих нули в аппаратуре РЛС и обработки. Трудно сказать однозначно, я не видел общей схемы всего этого хозяйства, этим занимались интеграторы. Например, на японских небольших станциях с питанием от +24В этой проблемы нет У наших станций она тоже не всегда есть, правда.

Может имеет смысл попробовать реализовать фильтр в виде IIR 32х битной точности ?

Тогда зачем отфильтровывать 50гц? Выкинте все от 0 до частоты зондирования, как минимум.

Похоже, нужен классический адаптивный компенсатор помех.
Читайте Уидроу Нажмите для просмотра прикрепленного файла
Но обеспечить 40 дБ в Ваших условиях будет непросто.

Methane

Насколько я понимаю, любой фильтр вносит задержку в сигнал. Если сигнал состоит из одной частоты, то можно так подобрать параметры фильтра, что задержка станет равна одному периоду, и сдвиг по фазе обнулится, так что сигнал станет пригоден для прямого вычитания. Если же фильтр выделяет совокупность частот, то период повторения совокупности уходит в бесконечность и подобрать параметрами фильтра нулевую задержку уже невозможно. Остается довольствоваться задержкой, какая есть. Так вот, даже если этот сдвиг будет, скажем, всего 5 градусов для 50 Гц, то это эквивалентно задержке на 15-18 тысяч отсчетов на основной частоте дискретизации видео. Мне вот жалко в ПЛИС столько памяти тратить на тупую задержку сигнала. То есть, вырезать спектр от 0 до 700 Гц не проще, чем от 50 до 60. Проще всего давить ровно одну частоту - 50 Гц. И так по одной и давить все неугодные помехи.

polyakovav
Спасибо, попробую полистать на буднях. Хотя мне для начала нужно оценить, какие подходы есть и какова трудоемкость каждого. В общем то сейчас уже есть фильтры, которыми эта помеха давится, но они работают уже после ВАРУ и после выбора полярности сигнала, а это уже поздновато.

Еще вот такой вопрос. Если мне нужно выделить всего одну частоту 50 Гц из сигнала, оцифрованного на 50-60 МГц, то мне явно нужно где-то 5 порядков частоты сбросить. Во сколько ступеней лучше делать децимацию такого рода? Насколько я понимаю, чем более резкие перепады в АЧХ, тем больше нужен порядок фильтра и больше разрядность коэффициентов. А в данном вопросе интересно сократить частоту дискретизации во много раз, угробив на это минимум ресурсов ПЛИС.

Если помеха действительно гармоническая, то адаптивный компенсатор не является фильтром в прямом смысле, и понижение частоты не требуется.
Можно формировать цифровым синтезатором в ПЛИС копию помехи с подстройкой по амплитуде и частоте (фазе), и вычитать из входной смеси сигнала с помехой.
Адаптация амплитуды и частоты (фазы) - по минимуму выходной энергии. Методы адаптации описаны в книжке.
В молодости применял для фильтрации сетевой помехи в кардиографии. При стабильных характеристиках помехи работает хорошо. Если помеха меняется во времени, важно правильно выбрать коэффициент адаптации (скорость сходимости).

Поддерживаю такую идею. Была как-то такая задача и решена была именно таким способом. Вполне успешно

Судя по ващим ответам, эквалайзер, невозможен в принципе. И зачем делать в ПЛИС? Просто из результатов по дальности вычтите задержку. Одна операция.


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

Вполне жизненный способ, только пункты 1 и 3 реализуются в пункте 2, и называться этот пункт будет просто «децимация сигнала».




Используйте КИХ фильтры с линейной фазой и счастье Вам гарантировано.

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

Нажмите для просмотра прикрепленного файла



Задержка все равно будет нужна, как минимум на половину порядка эквивалентного фильтра. Но ее легко реализовать на блочной памяти.
И таки да, вырезать спектр от 0 до 700 Гц проще, чем от 50 до 60.



Наименее затратно (по ресурсам) делать децимацию (и интерполяцию), используя максимально возможное количество ступеней. Но зачастую можно пойти на компромиссы.

1-3 вроде просто Герцелем (интегралом Фурье) можно сделать чтобы понять амплитуду и фазу 50Гц.

а то и просто завести цифрвой осциллятор на ~50Гц
y[n+1] = 2*cos(w) * y[n] - y[n-1]
сделать фапч, плавно подкручивая ему частоту и амплитуду подстраиваясь под помеху, и вычитать потом из сигнала.

Так вот в том то и дело, что на какой частоте дискретизации делать этот осциллятор? На частоте сигнала неудобно:
1) Как я писал выше, частота сигнала у меня может дискретно меняться от 25 до 80 МГц. Синтезатор хочется попроще, на одной частоте.
2) Для частоты сигнала, ну даже в 25 МГц, имеем cos(w) = cos( 6.28 * 50 / 25E6) = 1 - 8E-11. Ну то есть, это нужно считать в какой-то очень безумной разрядности.

Если теперь рассмотреть, что ФАПЧ и все такое все же лучше делать на какой-то промежуточной частоте, скажем, 1 КГц, то в полный рост вылезают снова проблемы с децимацией и апсемплингом.
И опять остается открытым вопрос, во сколько ступеней проще скинуть частоту с 25-80МГц до 1Кгц, чтобы фильтр был низкого порядка и разрядность коэффициентов не больше 12-16.

Да осциллятор можно проще сделать - ввиде таблички синуса и аккумулятора, старшие N разрядов аккумулятора == индекс элемента в таблице. Число на входе аккумулятора это частота. Как в DDS делают.

Можно использовать режекторный фильтр на основе всепропускающего фильтра
H(z)= 1+(a – bz^-1 + z^-2)/(1 - bz^-1 + az^-2),
который имеет плавную настройку частоты в пределах 0,02 - 0,45 частоты дискретизации и добротности.
Причем добротность может быть очень высока.
Если имеем частоту дискретизации 50 МГц,
то чтобы отрезать частоту 50 Гц надо в нем все регистровые задержки увеличить в 20000 раз -
и получится гребенчатый режекторный фильтр с 20000 провалами.
Поскольку входной сигнал более-менее случайный, то лишние провалы не повлияют на результат.
Похожие аналоговые гребенчатые фильтры используют в радиолокации очень давно.

Не всегда удобно вырезать из сигнала какую-то компоненту. Бывает, что проще подавить всё остальное, оставив необходимую частоту. Посмотрите, возможно, имеет смысл задавить весь сигнал на высоких, убрать НЧ и оставить только 50 Гц. Разумеется, при этом 50 Гц после данного фильтра приобретёт какую-то задержку по фазе. Параллельно сдвигаем исходный сигнал по фазе (на ОЗУ с заданной длиной кольцевого буфера, например) и на выходе вычитаем выделенную частоту 50Гц.

А вот я еще спрошу.

Мне вот понравилась идея, вытекающая из постов pv + alex_os. А именно, не делать никаких ФНЧ и ресемплингов, а просто взять и посчитать в квадратурах коэффициенты для sin и cos составляющих фильтруемой частоты. И вот я написал пару интегралов:
1) S(N) = Int([0..N*T0], sin(w*t)*sin(w0*t))
2) C(N) = Int([0..N*T0], cos(w*t)*sin(w0*t))

В идеале для очень большого N результат должен стремиться к интегралу Фурье. Но я пока построил графики для N=1, N=4, N=32. Вот они.
На графиках зеленым цветом показан интеграл 1, а синим цветом - интеграл 2. По горизонтали отложена безразмерная частота w/w0.
Видно, что
1) для w=w0 выполняется условие ортогональности sin(w0*t) и cos(w0*t).
2) для w=0 и для целых w/w0, больших 2, оба интеграла дают 0, что вполне законно.

Но вот при промежуточных значениях w, близких к w0, возникают паразитные отклики, которые выдают довольно большие амплитуды синуса и косинуса. Будет ли в этом случае наблюдаться эффект компенсации помехи, или, наоборот, фильтр начнет сам добавлять w0.

Можно ли для расчета амплитуд для синуса и косинуса делать измерения для нескольких разных N? Например, (S(8)*3 + S(4))/4 подавляет часть боковых лепестков на приведенном графике. Вначале я было подумал, что в этом есть смысл, но вот когда подал в качестве эталонного входа cos(wt) вместо sin(wt), то амплитуда sin-составляющей для некоторых стала w даже больше, чем при использовании только S(4) для оценки амплитуды sin-составляющей.

В принципе такая схема, получается, работает, если известно, что вблизи w0 в спектре нет других частот, так что отклика по ним фильтр не выдаст.

А зачем Вам сдались эти интегралы? Вам предлагают сделать ФАПЧ которая будет цеплятся за помеху в районе 50 Гц. ФАПЧ измеряет фазу и амплитуду помехи, далее полученную оценку помехи вычитаете из сигнала. Полоса пропускания ФАПЧ примерно равна удвоенной частоте среза разомкнутой петли. Соответственно режектор получится как бы полосовой. Сам генератор фапч работает на частоте дискретизации, но подстраивать фазу амплитуду можно и нужно с гораздо меньшей частотой, соответственно не будет никаких страшных фильтров для которых требуется не ординарная точность.

ну или не фапч, а просто посчитать интегралы Фурье только для w=50Гц чтобы понять амплитуду и фазу 50Гц помехи и затем вычесть.


да оно вроде более менее устойчиово ко всяким округлениям, ну а сделать всего пару умножений хоть и в 40+ разрядов имхо всё равно проще чем фильтры городить на 50Гц при 50МГц сэмплинга.

alex_os
Так если делать ФАПЧ, то проблема с частотой дискретизации не пропадает совсем, так как ФАПЧ содержит петлевой фильтр нижних частот, который должен интегрировать ошибку по фазе. Причем этот фильтр будет работать на частоте дискретизации сигнала, то есть 25-80 МГц, а выделять должен полосу от 0 до 2 Гц где-нибудь. Вот это место непонятно, как его сделать простым фильтром?

Когда я задумался об интегралах Фурье, то я понял так, что мы считаем параллельно два интеграла, получаем сразу два коэффициента перед синусом и косинусом, что эквивалентно амплитуде и фазе. Но то что мы считаем их на частоте дискретизации - это просто, скажем, 15-18 лишних бит к сумматору, на котором идет интегрирование, можно даже ничего не фильтровать и спектр не переносить. Вот. И это даже будет работать, если считать, что помеха ровно на одной частое w0.

С автоподстройкой, конечно, лучше было бы, но надо придумать, как петлевой фильтр реализовать на частоте дискретизации.

Так петлевой фильтр не нужно на частоте дискретизации делать он может содержать дециматоры, ну сделать его в виде скажем (первое что в голову пришло) СIC дециматора. Генератор ФАПЧ можно подстраивать один раз на тысячи входных отсчетов. Напрямую типа ДПФ можно подсчитать фазу и амплитуду, только что Вы будете делать дальше с этими данными?. Что бы воспользоваться измеренной фазой и амплитудой нужен будет здоровенный буфер размером время_измерения * частоту_дискретизации.

зачем буфер? интеграл Фурье считается на лету, а 50Гц более менее постоянные, то есть не будет особых проблем если, грубо говоря, секунду считать амплитуду и фазу 50Гц, а потом следующую секунду вычитать синус с этими параметрами, никуда он за секунду не денется.
плюс посчитанные интегралы Фурье можно наверное через какой-нибудь совсем простенький БИХ фильтр пропустить y+= (x-y) / K, чтобы сделать нечто вроде окна.
и тогда в любой момент времени есть фаза и амплитуда помехи. осталось только вычесть.

Как-то у вас все сложно с вашими фильтрами. А нельзя проще, примерно так.

Берем счетчик, который циклично работает и задает нужную нам частоту, равную, например 50 Гц * 32 точки=1600 Гц. Для 80 МГц это будет 50000 - всего 16 разрядов.
Берем сумматор, который суммирует все отсчеты сигнала в пределах периода счетчика. Если разрядность АЦП 10 - то разрядность сумматора будет порядка 36 (51200000) 26. В каждом периоде запихиваем старшие 10 бит сумматора в кольцевой буфер из 32 10-разрядных регистров или в 32 ячейки памяти.

В результате в буфере образуется 32-х точечный образ 50-герцовой помехи произвольной формы рассчитанный для предыдущих 32 точек. Его можно использовать для вычитания помехи в следующих 32 точках. Т.е.там будет не только синус помехи 50 Гц, но и ее гармоники, но с меньшей "точечностью", т.е. гармоника 100 Гц будет представлена 16 точками, зато дважды, и т.д. Небольшие отклонения частоты будут мешать, но немного. Можно не просто "вычислять" образ помехи, но и делать это с усреднением - новое значение в точке равно старому + разность нового и старого помноженная на некий коэффициент (или сдвинутая на несколько бит) (экспоненциальное усреднение). При сдвиге нам вообще не нужен перемножитель.

Далее вычитаем из сигнала в следующем периоде нашей тактовой частоты сигнал из регистра буфера, с адресом равным (текущему+1) - это сосчитанное в предыдущем "большом" периоде (32 точки) значение помехи 50 Гц для следующего такта.

При реализации на ПЛИС это гораздо проще, чем куча перемножителей или Фурье. Да и при реализации на процессоре, в общем, тоже. Особенно, если в основной задаче фурье не нужно.

P.S. пардон, по невнимательности наврал в разрядности сумматора - она будет 80 Мгц * 20 мсек / 32 точки *1024 = 51200000, т.е. всего 26 разрядов.

Вопрос - а от кадра к кадру у данного сигнала есть чистые от полезной составляющей поля? Гашение обратного хода луча или что-то из этой оперы? Может быть, "модель" помехи брать из этих интервалов и вычитать до следующего "обновления" модели?
А что если решить задачу по-пионерски железно: т.е. сгородить схемку оцифровки сетевого напряжения и для каждого канала установить необходимые фазу и амплитудный коэффициент компенсации? Отслеживать частоту наведёнки не придётся, а другие два параметра (фазовое смещении и коэффициент проникновения в сигнал) особо не плывут.

alex_os
Что-то меня идея с дециматором и классическим ФАПЧем отпугивает, почему-то кажется слишком навороченной.

Вот я тут еще на досуге подумал и вот какую штуку сочинил.
1) Идея о том, что синус формируем по табличке - правильная. Если табличка есть, то можем формировать синусоиду и косинусоиду по значениям фазы сигнала. Допустим, в ПЛИС у нас есть 8Кбит блочной памяти (у меня Cyclone III), это массив на 512 16-разрядных коэффициентов. Без особых трудностей можно считать, что на один период приходится 1000 точек фазы (не 1024 только потому, что удобнее делать 50Гц из N*5 МГц без дробной части; 500 точек на половину периода, вторая половина - это копия первой с обратным знаком). Итак, при частоте работы фильтра в 25МГц имеем интервал времени на приращение фазы такой: 25E6 / 50 / 1000 = 500 тактов. То есть, раз в 500 тактов мы можем увеличивать индекс и обновлять значения синусов и косинусов. В дальнейшем мы будем управлять периодом инкремента фазы, изменяя его в диапазоне от 480 до 520, что будет примерно соответствовать частоте от 48 до 52 ГЦ.

2) Вместо сложного ФАПЧ будем следить за тем, чтобы за N оборотов фазы на текущей частоте интеграл C(N) = Int([0..N*T(w)], a(t)*cos(w*t+p)) был равен нулю. Идея в том, что для любой конкретной частоты синус ортогонален косинусу, и если мы подберем фазу так, чтобы интеграл от сигнала, умноженного на косинус фазы был равен нулю, то это будет означать, что мы как раз синус фазы совпадает с сигналом (с точностью до знака). А значит, по значению интеграла от синуса можно проводить вычисление амплитуды синфазной составляющей. И теперь остается только решить вопрос с подстройкой фазы.

3) Подстройку фазы можно проводить по значению интеграла C(N). Как видно из графиков, которые я постил выше (синие кривые), интеграл sin(w*t)*cos(w0*t) дает значения разных знаков в некоторой окрестности w0. Значит, по знаку интеграла C(N) можно производить изменение интервала приращения фазы, при этом, однако, удерживая его в диапазоне регулирования от 480 до 520. Можно, наверное, ускорить сходимость фазы к фазе сигнала, но для фиксированной стационарной частоты помехи это, наверное, особо и не требуется.

4) Таким образом, подгоняя интеграл C(N) близким к нулю, можно обеспечить синхронность фазы и тогда интеграл S(N) = Int([0..N*T(w)], a(t)*cos(w*t+p)) будет давать значение A*N*T(w)/2, где A - амплитуда помехи на частоте w. И вот измерив A можно будет компенсировать помеху дальше.

Значение N можно выбрать побольше. Таким, чтобы полупериод на графике C(N) покрывал диапазон регулирования. Для диапазона 50+-2Гц, это достигается примерно для N=16.

Сложность с фильтрами и децимацией уходит вообще, остается только нормировка для определения A, так как при изменении интервала приращения фазы изменяется число слагаемых в интеграле - N*T(w). Так вот, если считать, что T(w) ~ (500+-20), то A ~ 1/(500+-20). Так как 20 сильно меньше 500, по можно считать по формуле 1/500*(1-+(20/500)). А учитывая, что даже на краю диапазона поправка в нормировке дает изменение всего в -28дБ к основному значению, то ее можно считать сильно упрощенно, и даже вообще не учитывать.

Как идея?

значение фазы прямо вычисляется из интегралов Фурье, впрочем как и величина амплитуды.
for(i...){
re += y[i] * cos(2_pi_w * t[i]);
im += y[i] * sin(2_pi_w * t[i]);
}
phase = atan(im/re);

pv
Интегралы Фурье на конечном отрезке времени работают только для периодических функций. Я же выше графики приводил, аж три штуки. Из них видно, что нельзя разложить близкие к w0 частоты на суперпозицию синуса и косинуса частоты w0. Поэтому бессмысленно считать фазу, если не уметь подстраивать частоту в некотором диапазоне. А если частота помехи строго 50 Гц, то фазу как таковую нет смысла считать, - можно просто посчитать амплитуды для синуса и косинуса, и потом эту суперпозицию и вычесть.