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 в спектре нет других частот, так что отклика по ним фильтр не выдаст.

Сообщение отредактировал Hoodwin - Jun 28 2013, 10:44

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

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


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

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

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

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

Сообщение отредактировал Hoodwin - Jun 28 2013, 14:04

Так петлевой фильтр не нужно на частоте дискретизации делать он может содержать дециматоры, ну сделать его в виде скажем (первое что в голову пришло) С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);