09-1 Глава 5 Случайные сигналы (Лекционный курс)

5.1Частотные характеристики цифровых фильтров

До сих пор мы не рассматривали такую важную характеристику ЦФ, как частотная. В какой-то мере роль частотной характеристики играет системная функция. Однако использовать только системную функцию не всегда удобно, особенно при анализе цифровых систем обработки сигналов, аналогичных непрерывным фильтрам.

Частотная характеристика (частотный коэффициент передачи) аналогового фильтра может быть определена как отношение спектра выходного сигнала к спектру входного сигнала (для -функции спектральная плотность S₁()=1):

K() = S₂()/S₁()

В качестве входного сигнала удобно использовать обобщенный гармонический сигнал e^jt, тогда выходной сигнал будет иметь вид:

K() e^jt

Для определения ЧХ ЦФ на его вход (точнее, на вход дискретного фильтра, т.к. эффекты квантования не учитываются) подают дискретизированный сигнал x(nT)= e^jnT. Тогда выходной сигнал будет иметь вид: y(nT)=K() e^jnT.

Поскольку входной и выходной сигналы дискретны во времени, а дискретные сигналы обладают периодическими спектрами с периодом, равным частоте дискретизации, то отношение их спектров даст также периодическую частотную характеристику с тем же периодом 2/T (рад/с).

Подставим обобщенный дискретизированный гармонический сигнал x(nT)=e^jnT в разностное уравнение (4.33) для цифрового фильтра общего вида:

K()e^jnT(1-b₁e^-jT-a₂e^-j2T...-b_Ne^-jNT) =

= e^jnT(a₀+a₁e^-jT+a₂e^-j2T+...+a_Me^-jMT)

Откуда найдем частотную характеристику:

(4.40)

Это выражение совпадает с выражением для системной функции цифрового фильтра, если
e^-jT заменить на z^-1. Таким образом, установлена простая связь между частотной характеристикой цифрового фильтра и его системной функцией:

K() = H(e^jT) (4.41)

Из (4.41) следуют многие свойства частотной характеристики цифровых фильтров, в частности, периодичность.

Соотношение (4.41) позволяет определить положение нулей и полюсов частотной характеристики K() по известному положению особых точек системной функции H(z).

(4.41) позволяет записать связь между частотной и импульсной характеристиками цифрового фильтра. Заменяя z на e^jT в формуле (4.26), получим:

(4.42)

Видно, что частотная характеристика связана с системной функцией соотношением, подобным Дискретному Преобразованию Фурье.

В качестве примера рассмотрим цифровой фильтр 1^го порядка с импульсной характеристикой:

g(nT)=e^-nT/,

для которого системная функция была найдена ранее в виде:

Заменяя z^-1 на e^-jT, получим выражение для частотной характеристики:

Рис. 4.1. АЧХ ЦФ: 1-при T/=1, 2-при T/=1/4. Для |f|<1/(2T) график ЧХ хорошо совпадает с ЧХ фильтра-прототипа (3)

Теперь найдем амплитудно-частотную и фазочастотную характеристики:

Определим положение особых точек передаточной функции K(p) рассмотренного цифрового фильтра. Единственный полюс системной функции этого фильтра расположен в точке z=e^-T/. Заменяя z на e^pT, получим e^pT= e^-T/, откуда:

Рис. 4.2. Расположение полюсов передаточной функции

Единственному полюсу системной функции соответствует бесконечное число полюсов передаточной функции, расположенных вдоль прямой j с интервалом 2/T (см. Рис. 4 .2).

5.2 Цифровая фильтрация в задачах восстановления пропущенных данных и интерполяции

По различным причинам в длинной записи данных могут быть пропущены одно или несколько изолированных значений. Например:

- измерение вообще не производилось;

- измерение сделано неправильно или запись измерения искажена и артефакты впоследствии удалены;

- при передаче цифровых сигналов по радиоканалу зафиксированы ошибки в изолированных отсчетах сигнала;

- при обработке сигналов была получена неопределенность типа sin(x)/x при x=0 и ЭВМ отказалась делить на нуль.

Обычным является способ восстановления пропущенных значений, который основан на предположении, что отрезок сигнала, содержащий пропуск, является многочленом (полиномом) P_n(x) некоторой степени n (обычно выбирают нечетные степени).

Для получения интерполяционной формулы полезна теорема, которая утверждает, что разностный оператор (n+1)-ой степени ⁽ⁿ⁺¹⁾ аннулирует полином P_n(x) степени n, т.е. разность (n+1)-го порядка для полинома P_n(x) тождественно равна нулю. Не доказывая эту теорему, исследуем ее справедливость для полиномов 0,1 и 2 степеней.

Запишем разностный оператор в виде:

[u_m]= u_m+1 - u_m

[P_n(mT)]= P_n(mT+T) - P_n(mT)

Полезно отметить, что конечно-разностная аппроксимация производной для дискретных систем с равномерным шагом дискретизации T определяется через разностный оператор:

u'|_m= (u_m+1 - u_m)/T = [u_m]/T

А) Полином нулевой степени - это прямая, параллельная оси абсцисс. Поскольку значение полинома неизменно, то первая же разность равна нулю:

[P₀(x)] 0

Б) Для полинома первой степени характерна постоянная скорость изменения значения (неизменный угол наклона), следовательно, для любых x разности при равномерном шаге дискретизации T окажутся одинаковыми. Значит, любой полином первого порядка после применения к нему разностного оператора переходит в полином нулевого порядка:

[P₁(x)] P₀(x)

Следующая разность оказывается нулевой (смотри п.А):

[[P₁(x)]]  ²P₁(x)  0

В) Для полинома второй степени P₂(x) определим первую разность в общем виде и получим полином первого порядка P₁(x):

Пусть P₂(x) = Ax² + Bx + C

Тогда: P₂(x+T) = A(x+T)² + B(x+T) + C

[P₂(x)] = P₂(x+T) - P₂(x) = A[(x+T)²-x²]+BT

[P₂(x)] = A(2Tx+T)+BT = A₁x+B₂ = P₁(x)

Следующие две разности "уничтожают" остаток (смотри пп.Б,А), т.е.:

³P₂(x)  0

Задание для самостоятельной работы: доказать теорему для полиномов 3-го порядка.

Предположим, что пропущенное значение принадлежит полиному третьего порядка. Следовательно, разность четвертого порядка для этого полинома должна тождественно равняться нулю: ⁴P₃(x)  0.

Запишем формулы вычисления разностей 2-го, 3-го и 4-го порядков для полинома P(x):

[P(x)]=P(x+T) - P(x)

[P(x+T)]=P(x+2T) - P(x+T)

²[P(x)]=P(x+2T) - 2P(x+T) + P(x)

²[P(x+T)]=P(x+3T) - 2P(x+2T) + P(x+T)

³[P(x)]=P(x+3T) - 3P(x+2T) + 3P(x+T) - P(x)

³[P(x+T)]=P(x+4T) - 3P(x+3T) + 3P(x+2T)-

- P(x+T)

⁴[P(x)]=P(x+4T) - 4P(x+3T) + 6P(x+2T) -

- 4P(x+T) + P(x)

Последнее выражение приравняем к нулю, изменим начало отсчета аргумента на -2T и разрешим относительно центральной координаты:

P(x+2T)-4P(x+T)+6P(x)-4P(x-T)+P(x-2T) 0

P(x)=[-P(x-2T)+4P(x-T)+4P(x+T)-P(x+2T)]/6

Получена простая, устойчивая и очень удобная формула для вычисления пропущенных значений по соседним отсчетам.

_{Примечание: Коэффициент усиления помех фильтром равен 34/36.}

Отметим, что нахождение единственного пропущенного значения отличается от процесса динамического прогнозирования данных (как, например у функции predict(M,c,N) в MathCad).

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

Для получения передаточной функции, как и прежде, в качестве входного подставим дискретизированный обобщенный гармонический сигнал exp(jnT):

H()=[-e^(-j2T)+4e^(-jT)+4e^(jT)-e^(j2T)]/6

Упрощение выражения дает результат:

H()=[4cos()-cos(2)]/3

Если интерполированное значение в точности равно пропущенному, то передаточная функция принимает значение, равное 1. Легко видеть, что правильный ответ получается для нулевой частоты.

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

Рис. 4.3. Передаточная функция для интерполяции пропущенных данных полиномом 3-го порядка. Отрицательные значения означают изменение фазы на 180º или смену знака

Такие же результаты получаются при интерполяции по соседним значениям с помощью полиноминального приближения по методу наименьших квадратов.

Задача интерполяции промежуточных значений (между регулярными выборками) немного отличается от рассмотренной выше.

Пусть требуется интерполировать значение дискретного сигнала u(nT) между двумя соседними отсчетами с координатами nT и (n+1)T.

Обозначив середину текущего интервала дискретизации величиной mT, получим для полинома первой степени (случай линейной интерполяции) простейший алгоритм интерполяции:

mT=(n+1/2)T

u(mT)= [u(mT-T/2)+ u(mT+T/2)]/2