Режекторный фильтр на 50 Гц для видеосигнала, как удалить наводку от цепей питания на видеосигнал Опции

[Hoodwin]

Собственно возникла вот какая задача. Имеем видеосигнал АЦП с частотой 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. Возможно ли это?

[Hoodwin]

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дБ к основному значению, то ее можно считать сильно упрощенно, и даже вообще не учитывать.

Как идея?

[_pv]

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