Нет ничего, сколь бы великим и изумительным оно не показалось с первого взгляда, на что мало-помалу не начинаешь смотреть с меньшим изумлением.
Лукреций. Римский поэт и философ, I в.д.н.э.
Цифровая обработка сигналов оперирует с дискретными преобразованиями сигналов и обрабатывающих данные сигналы систем. Математика дискретных преобразований зародилась в недрах аналоговой математики еще в 18 веке в рамках теории рядов и их применения для интерполяции и аппроксимации функций, однако ускоренное развитие она получила в 20 веке после появления первых вычислительных машин. В принципе, в своих основных положениях математический аппарат дискретных преобразований подобен преобразованиям аналоговых сигналов и систем. Однако дискретность данных требует учета этого фактора, и его игнорирование может приводить к существенным ошибкам. Кроме того, ряд методов дискретной математики не имеет аналогов в аналитической математике.
Дискретное преобразование Фурье может быть получено непосредственно из интегрального преобразования дискретизаций аргументов (tk = kDt, fn = nDf):
S(f) = s(t) exp(-j2pft) dt, S(fn) = Dts(tk) exp(-j2pfnkDt), (6.1.1)
s(t) =S(f) exp(j2pft) df, s(tk) = DfS(fn) exp(j2pnDftk). (6.1.2)
Напомним, что дискретизация функции по времени приводит к периодизации ее спектра, а дискретизация спектра по частоте - к периодизации функции. Не следует также забывать, что значения (6.1.1) числового ряда S(f
n) являются дискретизаций непрерывной функции S'(f) спектра дискретной функции s(tk), равно как и значения (6.1.2) числового ряда s(tk) являются дискретизацией непрерывной функции s'(t), и при восстановлении этих непрерывных функций S'(f) и s'(t) по их дискретным отсчетам соответствие S'(f) = S(f) и s'(t) = s(t) гарантировано только при выполнении теоремы Котельникова-Шеннона.Для дискретных преобразований s(k
Dt) Ы S(nDf), и функция, и ее спектр дискретны и периодичны, а числовые массивы их представления соответствуют заданию на главных периодах Т = NDt (от 0 до Т или от -Т/2 до Т/2), и 2fN = NDf (от -fN до fN), где N – количество отсчетов, при этом:D
f = 1/T = 1/(NDt), Dt = 1/2fN = 1/(NDf), DtDf = 1/N, N = 2TfN. (6.1.3)Соотношения (6.1.3) являются условиями информационной равноценности динамической и частотной форм представления дискретных сигналов. Другими словами: число отсчетов функции и ее спектра должны быть одинаковыми. Но каждый отсчет комплексного спектра представляется двумя вещественными числами и, соответственно, число отсчетов комплексного спектра в 2 раза больше отсчетов функции? Это так. Однако представление спектра в комплексной форме -
не более чем удобное математическое представление спектральной функции, реальные отсчеты которой образуются сложением двух сопряженных комплексных отсчетов, а полная информация о спектре функции в комплексной форме заключена только в одной его половине - отсчетах действительной и мнимой части комплексных чисел в частотном интервале от 0 до fN, т.к. информация второй половины диапазона от 0 до -fN является сопряженной с первой половиной и никакой дополнительной информации не несет.При дискретном представлении сигналов аргумент t
k обычно проставляется номерами отсчетов k (по умолчанию Dt = 1, k = 0,1,…N-1), а преобразования Фурье выполняются по аргументу n (номер шага по частоте) на главных периодах. При значениях N, кратных 2:S(fn) є Sn = sk exp(-j2pkn/N), n = -N/2,…,0,…,N/2. (6.1.4)
s(tk) є sk = (1/N)Sn exp(j2pkn/N), k = 0,1,…,N-1. (6.1.5)
Главный период спектра в (6.1.4) для циклических частот от -0.5 до 0.5, для угловых частот от -
p до p. При нечетном значении N границы главного периода по частоте (значения ± fN) находятся на половину шага по частоте за отсчетами ± (N/2) и, соответственно, верхний предел суммирования в (6.1.5) устанавливается равным N/2.В вычислительных операциях на ЭВМ для исключения отрицательных частотных аргументов (отрицательных значений номеров n) и использования идентичных алгоритмов прямого и обратного преобразования Фурье главный период спектра обычно принимается в интервале от 0 до 2f
N (0 Ј n Ј N), а суммирование в (6.1.5) производится соответственно от 0 до N-1. При этом следует учитывать, что комплексно сопряженным отсчетам Sn* интервала (-N,0) двустороннего спектра в интервале 0-2fN соответствуют отсчеты SN+1-n (т.е. сопряженными отсчетами в интервале 0-2fN являются отсчеты Sn и SN+1-n).Пример:
На интервале Т= [0,99], N=100, задан дискретный сигнал s(k) =d(k-i) - прямоугольный импульс с единичными значениями на точках k от 3 до 8. Форма сигнала и модуль его спектра в главном частотном диапазоне, вычисленного по формуле S(n) =s(k)Ч exp(-j2pkn/100) с нумерацией по n от -50 до +50 с шагом по частоте, соответственно, Dw=2p/100, приведены на рис. 6.1.1.Рис. 6.1.1. Дискретный сигнал и модуль его спектра.
На рис. 6.1.2 приведена огибающая значений другой формы представления главного диапазона спектра. Независимо от формы представления спектр периодичен, в чем нетрудно убедиться, если вычислить значения спектра для большего интервала аргумента n с сохранением того же шага по частоте, как это показано на рис. 6.1.3 для огибающей значений спектра.
Рис. 6.1.2. Модуль спектра. Рис. 6.1.3. Модуль спектра.
На рис. 6.1.4. показано обратное преобразование Фурье для дискретного спектра, выполненное по формуле s'(k) =(1/100)
S(n)Ч exp(j2pkn/100), которое показывает периодизацию исходной функции s(k), но главный период k={0,99} этой функции полностью совпадает с исходным сигналом s(k).Рис. 6.1.4. Обратное преобразование Фурье.
Преобразования (6.1.4-6.1.5) называют дискретными преобразованиями Фурье (ДПФ). Для ДПФ, в принципе, справедливы все свойства интегральных преобразований Фурье, однако при этом следует учитывать периодичность дискретных функций и спектров. Произведению спектров двух дискретных функций (при выполнении каких-либо операций при обработке сигналов в частотном представлении, как, например, фильтрации сигналов непосредственно в частотной форме) будет соответствовать свертка периодизированных функций во временном представлении (и наоборот). Такая свертка называется циклической (см. раздел 6.4) и ее результаты на концевых участках информационных интервалов могут существенно отличаться от свертки финитных дискретных функций (линейной свертки).
Из выражений ДПФ можно видеть, что для вычисления каждой гармоники нужно N операций комплексного умножения и сложения и соответственно N
2 операций на полное выполнение ДПФ. При больших объемах массивов данных это может приводить к существенным временным затратам. Ускорение вычислений достигается при использовании быстрого преобразования Фурье.Быстрое преобразование Фурье (БПФ, fast Fourier transform - FFT). Оно базируется на том, что при вычислениях среди множителей (синусов и косинусов) есть много периодически повторяющихся значений (в силу периодичности функций). Алгоритм БПФ группирует слагаемые с одинаковыми множителями в пирамидальный алгоритм, значительно сокращая число умножений за счет исключения повторных вычислений. В результате быстродействие БПФ в зависимости от N может в сотни раз превосходить быстродействие стандартного алгоритма. При этом следует подчеркнуть, что алгоритм БПФ даже точнее стандартного, т.к. сокращая число операций, он приводит к меньшим ошибкам округления.
Допустим, что массив чисел s
k содержит N = 2r отсчетов (r - целое). Разделим исходный массив на два первых промежуточных массива с четными и нечетными отсчетами:sk' = s2k, sk" = s2k+1, 0 Ј k Ј N/2-1.
Выполним ДПФ каждого массива с учетом того, что шаг функций равен 2 (при
Dt=1), а период промежуточных спектров будет соответственно равен N/2:sk' Ю Sn', sk" Ю Sn", 0 Ј n Ј N/2-1.
Для получения одной половины искомого спектра S
n сложим полученные спектры с учетом теоремы запаздывания, т.к. отсчеты функции sk" сдвинуты относительно sk' на один шаг дискретизации:Sn = Sn'+Sn"Ч exp(-j2pn/N). (6.1.6)
Вторая половина спектра, комплексно сопряженная с первой, с учетом периода повторения N/2 промежуточных спектров определяется выражением:
Sn+N/2 = Sn'+Sn"Ч exp(-j2p(n+N/2)/N) = Sn'- Sn"Ч exp(-j2pn/N). (6.1.7)
Нетрудно видеть, что для вычисления полного спектра в данном случае потребуется N
2/4 операций для вычисления промежуточных спектров плюс еще N операций комплексного умножения и сложения, что создает ощутимый эффект по сравнению с ДПФ.Но деление массивов на две части может быть применено и к первым промежуточным массивам, и ко вторым, и т.д. до тех пор, пока в массивах не останется по одному отсчету, фурье - преобразование которых равно самому отсчету. Тем самым, алгоритм преобразования превращается в пирамидальный алгоритм перестановок со сложением/вычитанием и с единичным умножением на значение exp(-j2
pn/N) соответствующего уровня пирамиды. Первый алгоритм БПФ на данном принципе (из множества модификаций, существующих в настоящее время) был разработан Кули-Тьюки в 1965 г. и позволил повысить скорость вычислений в N/r раз по сравнению с ДПФ. Чем больше N, тем больше эффект БПФ. Так, при N = 1024 имеем r = 10 и соответственно N/r » 100. Что касается условия по количеству точек N = 2r, то оно рассматривается в варианте Nk Ј 2r, где r - минимальное целое. Массивы с Nk < 2r дополняется до 2r нулями, что не изменяет форму спектра. Изменяется только шаг Dw по представлению спектра (Dw = 2p/2r < 2p/N), который несколько избыточен по адекватному представлению сигнала в частотной области. В настоящее время существуют и алгоритмы БПФ с другими основаниями и их комбинациями, при которых не требуется дополнения сигналов нулями до 2r.Заметим, что в соответствии с (6.1.7) отсчеты, сопряженные с правой половиной главного частотного диапазона (0,
p ), относятся не к диапазону (-p ,0), а к диапазону (p ,2p ), что, учитывая периодичность спектра дискретных данных, значения не имеет. Т.е. выходной частотный диапазон БПФ равен (0, 2p ). Общее количество отсчетов комплексного спектра в этом условно главном диапазоне равно количеству точек исходного сигнала (с учетом нулевых точек при дополнении сигнала до N=2r). Алгоритм быстрого обратного преобразования (ОБПФ) тождественен алгоритму прямого БПФ.Алгоритмы прямого и обратного БПФ широко используются в современном программном обеспечении для анализа и обработки цифровых данных. Пример выполнения БПФ приведен на рис. 6.1.5.
Рис. 6.1.5. Пример БПФ.
Основная область использования ДПФ – спектральный анализ физических данных. При этом интерес обычно представляют только амплитуды отдельных гармоник, а не их фазы, и спектр отображается в виде графика зависимости амплитуды (модуля спектра) от частоты. Часто шкала амплитуд градуируется в децибелах. Децибелы - логарифмы отношения амплитудных значений. Например, разница на 20 дБ означает различие амплитуд в 10 раз, разница на 40 дБ - 100 раз. Различию амплитуд в 2 раза отвечает разница примерно в 6 дБ. Шкала частот также часто градуируется в логарифмическом масштабе.
Перед вычислением спектра из сигнала, как правило, вырезается отрезок сигнала. Число последовательных отсчетов отрезка для использования БПФ должно быть степенью двойки, если в программном обеспечении вычислительной системы не оговорена ее способность выполнять БПФ по произвольным числовым рядам. В противном случае числовой ряд дополняется нулями до необходимого размера, что не изменяет формы спектра и сказывается только на увеличении частотного разрешения по спектру.
При вычислении спектра возможен следующий нежелательный эффект. При разложении участка сигнала в ряд Фурье мы тем самым принимает этот участок за один период Т, который периодически повторяется за пределами участка с фундаментальной частотой 1/Т. При ДПФ, а равно и при БПФ, вычисляется спектр именно такого периодического сигнала. При этом на границах периодов такая функция наверняка будет иметь разрывы или скачки, тем самым существенно искажая спектр. Для устранения этого эффекта применяются так называемые весовые окна, похожие на гауссиан, размер которых равен размеру участка. Анализируемый участок умножается на весовое окно, что плавно сводят сигнал на нет вблизи краев анализируемого участка и в значительной степени устраняют рассмотренные искажения спектра. Методика применения весовых окон подробно рассматривается в курсе цифровой обработки сигналов.
Дискретное преобразование Лапласа (ДПЛ), как и ДПФ, может быть получено из интегрального преобразования дискретизаций аргументов (tk = kDt, wn = nDw):
Y(p) =y(t) exp(-pt) dt, Y(pn) = Dty(tk) exp(-pntk), (6.2.1)
где p =
s+jw - комплексная частота, s і 0.y(t) = (1/2pj) Y(p) exp(pt) dp. y(tk) = DtY(pn) exp(pntk). (6.2.2)
Функцию Y(p) называют изображением Лапласа функции y(t) - оригинала изображения. Нетрудно заметить, что при
s = 0 преобразование Лапласа превращается в одностороннее преобразование Фурье, а для каузальных сигналов - в полную аналогию ПФ. Наиболее существенной особенностью преобразования Лапласа является возможность его применения для спектрального анализа функций, не имеющих фурье-образов из-за расходимости интегралов Фурье. Для понимания последнего запишем интеграл Лапласа в развернутой форме:Y(p) =y(t) exp(-st-jwt) dt =y(t) exp(-st) exp(-jwt) dt =y'(t) exp(-jwt) dt.
Правый интеграл для каузальных сигналов представляет собой преобразование Фурье, при этом сам сигнал y'(t) за счет экспоненциального множителя exp(-
st) соответствующим выбором значения s>0 превращается в затухающий, конечный по энергии, что и требуется для существования его фурье-образа. Все свойства и теоремы преобразований Фурье имеют соответствующие аналоги и для преобразований Лапласа.Пример сопоставления преобразований Фурье и Лапласа приведен на рис. 6.2.1.
Рис. 6.2.1. Сопоставление преобразований Фурье и Лапласа.
6.3. Z - преобразование сигналов [2,13,21].
Определение преобразования.
Распространенным способом анализа дискретных цифровых последовательностей является z-преобразование (z-transform). Оно играет для дискретных сигналов и систем такую же роль, как для аналоговых – преобразование Лапласа.Произвольной непрерывной функции s(t), равномерно дискретизированной и отображенной отсчетами s
k = s(kDt), равно как и непосредственно дискретной функции, можно поставить в соответствие степенной полином по z, последовательными коэффициентами которого являются значения sk:sk = s(kDt) Ы TZ[s(kDt)] =sk zk = S(z). (6.3.1)
где z =
s+jv = rЧ exp(-jj) - произвольная комплексная переменная. Полином S(z) называют z-образом или z-изображением функции s(kDt). Преобразование имеет смысл для области тех значений z, в которой ряд S(z) сходится, т.е. сумма ряда представляет собой аналитическую функцию переменной z, не имеющую полюсов и особых точек.Пример:
sk = {1, 2, 0, -1, -2, -1, 0, 0}.S(z) = 1z0+2z1+0z2-1z3-2z4-1z5+0z6+0z7 = 1+2z-z3-2z4-z5.
Впервые z-преобразование введено в употребление П.Лапласом в 1779 и повторно "открыто" В.Гуревичем в 1947 году с изменением символики на z
-1. В настоящее время в технической литературе имеют место оба вида символики. На практическое использование преобразования это не влияет, т.к. смена знака только зеркально изменяет нумерацию членов полинома (относительно z0), числовое пространство которых в общем случае от -Ґ до +Ґ . В дальнейшем в качестве основной будем использовать символику положительных степеней z, давая пояснения по особенностям отрицательной символики, если таковая имеется.По заданному или полученному в результате анализа какой-либо системы z-полиному однозначно восстанавливается соответствующая этому полиному функция путем идентификации коэффициентов степеней при z
k с k-отсчетами функции.Пример:
S(z) = 1+3z2+8z3-4z6-2z7 = 1z0+0z1+3z2+8z3+0z4+0z5-0z6-2z7.sk = {1, 0, 3, 8, 0, 0, -4, -2}.
Смысл величины z в z-полиноме заключается в том, что она является оператором единичной задержки по координатам функции. Умножение z-образа сигнала s(k) на величину z
n означает задержку сигнала на n интервалов: znS(z) Ы s(k-n).Z-образы с положительными степенями z соответствуют каузальным (физически реализуемым) процессам и системам, которые работают в реальном масштабе времени с текущими и "прошлыми" значениями сигналов. При обработке информации на ЭВМ каузальность сигналов не относится к числу ограничений и возможно использование отрицательных степеней z, соответствующих отсчетам сигналов "вперед". Последнее применяется, например, при синтезе симметричных операторов фильтров, что позволяет производить обработку информации без внесения в сигнал фазовых искажений. При использовании символики z
-1 "прошлым" значениям соответствуют значения с отрицательными степенями z, "будущим – с положительными.Основное достоинство z-преобразований заключается в простоте математических операций со степенными полиномами, что имеет немаловажное значение при расчетах цифровых фильтров и спектральном анализе.
Примеры z-преобразования
часто встречающихся на практике дискретных сигналов.Импульсы Кронекера.
В общем случае, в произвольной точке числовой оси:d
(k-n) =1 при k=n, d(k-n) = 0 при k ≠ n.Xd(z) =d(k-n) zk = zn.
Для импульса Кронекера в нулевой точке соответственно X
d(z) = z0 =1.Функция Хевисайда
(единичный скачок).x(k) = 0 при k < 0, x(k) = 1 при k і 0.
X(z) =zk = zk.
Ряд сходится при |z| < 1, при этом его сумма равна:
X(z) = 1/(1-z), |z| < 1.
При использовании символики z
-1:X(z) = 1/(1-z-1) = z/(z-1), |z| > 1.
Экспоненциальная функция:
x(k) = 0 при k < 0, x(k) = ak при k і 0.
X(z) =x(k) zk = ak zk = (az)k.
Как и в предыдущем случае, ряд сходится при |az| > 1, т.е. при |z| < |a|, при этом:
X(z) = 1/(1-az), |z| > |a|.
Связь с преобразованиями Фурье и Лапласа.
Запишем дискретный сигнал sk в виде суммы весовых импульсов Кронекера:sk = s(kDt) =s(nDt) d(kDt-nDt).
Определим спектр сигнала по теореме запаздывания:
S(w ) =s(kDt) exp(-jwkDt).
Выполним замену переменных, z = exp(-j
wDt), и получим:S(w ) =s(kDt)Ч zk = S(z).
Отсюда следует, что дискретное преобразование Фурье является частным случаем z-преобразования при z = exp
(-jwDt). Аналогичной подстановкой z = exp(-p) может осуществляться переход к дискретному преобразованию Лапласа. В общем виде:S(w) = S(z), z = exp(-jwDt); S(p) = S(z), z = exp(-pDt). (6.3.2)
Обратное преобразование:
S(z) = S(w), w = ln z/jDt; S(z) = S(p), p = ln z/Dt. (6.3.3)
При отрицательной символике z связь между представлениями осуществляется соответственно подстановками z
-1 = exp(jwDt) и z-1 = exp(p).Свойства z-преобразования.
Без углубления в теорию, можно констатировать, что все свойства ДПФ действительны и для z-преобразования. Отметим некоторые из них.Линейность
: Если S(k) = a·x(k)+b·y(k), то S(z) = aX(z)+bY(z). Соответственно, z-преобразование допустимо только для анализа линейных систем и сигналов, удовлетворяющих принципу суперпозиции.Задержка
на n тактов: y(k) = x(k-n).Y(z) =y(k)Ч zk =x(k-n)Ч zk =znx(k-n)Ч zk-n = zn x(m)Ч zm = zn X(z).
Соответственно, умножение z-образа сигнала на множитель z
n вызывает сдвиг сигнала на n тактов дискретизации.Для z-преобразования действительны все известные теоремы о спектрах. В частности, свертка двух сигналов отображается в z-области произведением их z-образов, и наоборот:
s(k) * h(k) Ы S(z)H(z), s(k)·h(k) Ы S(z) * H(z).
При z = exp(-j
wDt) z-преобразование представляет собой особую форму представления дискретных сигналов, при которой на полином S(z) можно ссылаться как на временную функцию (по значениям коэффициентов kDt), так и на функцию частотного спектра сигнала (по значениям аргумента w).Отображение z-преобразования
выполняют на комплексной z-плоскости с Re z и Im z по осям координат (рис. 6.3.1). Спектральной оси частот w на z-плоскости соответствует окружность радиуса:|z| = |exp(-jwDt)| = = 1.
Подстановка значения какой-либо частоты
w в z = exp(-jwDt) отображается точкой на окружности. Частоте w = 0 соответствует точка Re z = 1 и Im z = 0 на правой стороне оси абсцисс. При повышении частоты точка смещается по окружности против часовой стрелки, и занимает крайнее левое положение на частоте Найквиста wN = p/Dt (Re z = -1, Im z = 0). Отрицательные частоты спектра отображаются аналогично по часовой стрелке на нижней полуокружности. Точки wN совпадают, а при дальнейшем повышении или понижении частоты значения начинают повторяться в полном соответствии с периодичностью спектра дискретной функции. Проход по полной окружности соответствует одному периоду спектра, а любая гармоника спектра сигнала задается на плоскости двумя точками, симметричными относительно оси абсцисс.Z-преобразование позволяет производить разложение сигналов и функций, например передаточных функций фильтров, на короткие составляющие операции свертки, для чего достаточно приравнять z-полином к нулю, найти его корни a
i, и переписать полином в виде произведения двучленов:S(z) = a0(z-a1)(z-a2)...,
где а
0- последний отсчет сигнала (коэффициент при старшей степени z).Но произведению в z-области соответствует свертка в координатной области, и при обратном преобразовании двучлены (z-a
i) превращаются в двухточечные диполи {-ai,1}, а сигнал длиной N представляется сверткой (N-1) диполей:sk= a0{-a1,1}*{-a2,1}*{-a3,1}* ...
Пример
. sk = {1.4464, -2.32, 3.37, -3, 1}. S(z) = z4-3z3+3.37z2-2.32z+1.4464. a0 = 1.Корни полинома S(z): a
1 = 0.8+0.8j, a2 = 0.8-0.8j, a3 = 0.7+0.8j, a4 = 0.7-0.8j,S(z) = (z-0.8-0.8j)(z-0.8+0.8j)(z-0.7-0.8j)(z-0.7+0.8j).
Корни полинома представлены на z-плоскости на рис. 6.3.1. Корни полинома комплексные и четыре двучлена в координатной области также будут комплексными. Но они являются сопряженными, и для получения вещественных функций следует перемножить сопряженные двучлены и получить биквадратные блоки: S(z) = (z
2-1.4z+1.13)(z2-1.6z+1.28).При переходе в координатную область: s
k = {1.13, -1.4, 1} * {1.28, -1.6, 1}.Таким образом, исходный сигнал разложен на свертку двух трехчленных сигналов (функций).
Аналитическая форма z-образов
существует для z-преобразований, если возможно свертывание степенного ряда в аналитическое выражение. Выше, в примерах z-преобразования, уже приводилось приведение к аналитической форме z-образов функции Хевисайда и экспоненциальной функции.Обратное z-преобразование
в общем случае производится интегрированием по произвольному замкнутому контуру C, расположенному в области сходимости и окружающему все особые точки (нули и полюсы) z-образа:sk = (1/2pj)
Способом, удобным для практического применения, является разложение рациональных S(z) на простые дроби. С учетом линейности преобразования:
S(z) =an/(1-bnz) Ы an(bn)k = sk.
Пример.
S(z) = 1/(1-5z+6z2) = 3/(1-3z)-3/(1-2z) Ы 3Ч 3k -3Ч 2k = s(k).При разложении функции S(z) по степеням z обратное z-преобразование не вызывает затруднений.
6.4. Дискретная свертка (конволюция) [5,17,21].
Свертка – основной процесс в цифровой обработке сигналов. Поэтому важно уметь эффективно ее вычислять.
Уравнение дискретной свертки
двух функций (сигналов) может быть получено непосредственно из интегрального уравнения свертки при замене интегрирования суммированием мгновенных значений функций с шагом Dt:y(kDt) = Dth(nDt) s(kDt-nDt). (6.4.1)
При выполнении дискретной свертки мы имеем дело с цифровыми массивами, при этом шаг дискретизации для массивов по физическому аргументу свертки должен быть равным и принимается за 1, а в качестве аргумента используется нумерация отсчетов в массивах:
y(k) =h(n) s(k-n) є hn sk-n є yk. (6.4.1')
y(k) = h(n) * s(k-n) є s(k) * h(n) є sk * hn.
Техника свертки
приведена на рис. 6.4.1. Для вычисления свертки массив одной из функций (sk- входного сигнала) располагается по ходу возрастания номеров. Массив второй функции (hn - более короткой, оператор свертки), строится параллельно первому массиву в обратном порядке (по ходу уменьшения номеров, в режиме обратного времени). Для вычисления yk значение h0 располагается против sk, все значения sk-n перемножаются с расположенными против них значениями hn и суммируются. Результаты суммирования являются выходным значением функции yk, после чего оператор hn сдвигается на один номер k вперед (или функция sk сдвигается ему навстречу) и вычисление повторяется для номера k+1 и т.д.В начальный момент свертки при вычислении значений y
k оператор hn, построенный в режиме обратного времени, "зависает" для значений k-n при n>k против отсутствующих отсчетов входной функции. "Зависание" исключают либо заданием начальных условий - дополнительных отсчетов, чаще всего нулевых или равных первому отсчету входной функции, либо началом свертки с отсчета входной функции k = n с соответствующим сокращением интервала выходной функции. Для операторов со значениями -n (вперед по времени) такой же момент может наступать и в конце входного массива.Пример.
Уравнение свертки: yk =bn xk-n = bo xk + b1 xk-1 + b2 xk-2. Значения оператора bn:bo = 5, b1 = 3, b2
= 2. Входной сигнал: xk = {0,1,0,0,0}, начальные условия: x-n = 0.Расчет выходного сигнала:
yo = 5xo + 3x-1+ 2x-2 = 5 · 0 + 3 · 0 + 2 · 0 = 0, y1 = 5x1 + 3xo + 2x-1 = 5 · 1 + 3 · 0 + 2 · 0 = 5,
y2 = 5x2 + 3x1 + 2xo = 5 · 0 + 3 · 1 + 2 · 0 = 3, y3 = 5x3 + 3x2 + 2x1 = 5 · 0 + 3 · 0 + 2 · 1 = 2,
y4 = 5x4 + 3x3 + 2x2 = 5 · 0 + 3 · 0 + 2 · 0 = 0, y5 = 5x5 + 3x4 + 2x3 = 5 · 0 + 3 · 0 + 2 · 0 = 0
Выходной сигнал: yk = {0, 5, 3, 2, 0}
Заметим: свертка функции оператора с единичным входным сигналом представляет собой повторение функции оператора свертки на выходе.
На рис. 6.4.2 приведен пример выполнения дискретной свертки каузальным (односторонним) и четным (симметричным, двусторонним) оператором одного и того же сигнала.
Рис. 6.4.2. Примеры выполнения дискретной свертки.
Прямое вычисление свертки требует K·N умножений, где K – длина исходного сигнала, а N – длина ядра свертки. Как длина сигнала, так и длина ядра свертки может достигать нескольких тысяч точек, и число умножений становится огромным.
Для дискретной свертки действительны все свойства и теоремы интегральной свертки. В частности, свертка функций в координатной области отображается произведением их спектров в частотной области, а умножение в координатной области эквивалентно свертке в частотной области. Это значит, что для выполнения свертки двух сигналов можно перевести их в частотную область, перемножить их спектры, и перевести результат обратно во временную область, т.е. действовать по следующей схеме:
s(k) Ы S(w), h(n) Ы H(w), Y(w) = S(w)Ч H(w), Y(w) Ы y(k).
С появлением алгоритмов БПФ, позволяющих быстро вычислять преобразования Фурье, вычисление свертки через частотную область стало широко использоваться. При значительных размерах сигналов и длины ядра свертки такой подход позволяет в сотни раз
сократить время вычисления свертки.
Выполнение произведения спектров может производиться только при одинаковой их длине, и оператор h(n) перед ДПФ необходимо дополнять нулями до размера функции s(k).
Второй фактор, который следует принимать во внимание, это цикличность свертки при ее выполнении в спектральной области, обусловленная периодизацией дискретных функций. Перемножаемые спектры являются спектрами периодических функций, и результат на концевых интервалах может не совпадать с дискретной линейной сверткой, где условия продления интервалов (начальные условия) задаются, а не повторяют главный период.
На рис. 6.4.3 приведены результаты свертки сигнала s
k, заданного на интервале k=(0-50), с функцией hn = aЧ exp(-aЧ n), a = 0.1. Свертка, выполненная через ДПФ, в левой части интервала резко отличается от линейной свертки. Характер искажения становится понятным, если дополнить главный интервал с левой стороны его периодическим продолжением (на рисунке показана часть левого бокового периода, свертка с которым заходит в главный период). Для операторов hn со значениями n, вперед по положению, аналогичные искажения появятся и в правой части главного периода. Для устранения таких искажений сигнальная функция должна продлеваться нулями на размер оператора h(n), что исключит наложение боковых периодов главной трассы функции. Следовательно, перед переводом функций s(k) и h(n) в спектральную область их размер должен продлеваться нулями до длины K+N при односторонних операторах h(n), или до длины K+2N при двусторонних операторах h(n).При выполнении свертки через БПФ ощутимое повышение скорости вычислений появляется только при большой длине функций и операторов (например, M>1000, N>100). Следует также обращать внимание на разрядность результатов, т.к. перемножение чисел дает увеличение разрядности в 2 раза. При ограниченной разрядности числового представления с соответствующим округлением это может приводить к погрешностям суммирования.
В системах оперативной обработки данных часто возникает потребность вычислить свертку очень сигнала, поступающего на вход системы последовательными порциями (например, при получении данных от датчиков скважинных приборов). В таких случаях применяется так называемая секционная свертка. Суть ее состоит в том, что каждая из этих частей сворачивается с ядром отдельно, а затем полученные части объединяются. Для объединения достаточно размещать их друг за другом с наложением (перекрытием) в N-1 точку (N – длина ядра свертки), и производить суммирование в местах перекрытия.
Copyright ©2006 Davydov А.V.