Цифровая обработка сигналов

техника, 2012. – 448 с.
2. Ануфриев И., Смирнов А., Смирнова Е. MATLAB 7 (+CD-ROM). – СПб.: БХВ-Петербург, 2005. – 1102 с.
3. Солонина А.И., Арбузов С.М. Цифровая обработка сигналов. Моделирование в MATLAB. – СПб.: БХВ-Петербург, 2008. – 816 с.
4. Сергиенко А.Б. Цифровая обработка сигналов. 3-е издание. – СПб.: БХВ-Петербург, 2011. – 758 с.
5. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. – М.: Мир, 1978. – 848 с.
6. Stanley W.D. Digital signal processing. – Englewood, USA: Reston Publishing Co., Inc., 1975. – 323 p.

Простота перестройки алгоритмов.
3. Неограниченная временная стабильность характеристик.
4. Простота обработки низкочастотных (НЧ) сигналов.
5. Высокая стабильность и надежность работы и хорошая воспроизводимость получаемых результатов.
6. Неограниченную во времени и практически неограниченную по объему память для запоминания сигналов и параметров.
7. Снижается риск самовозбуждения системы и открывает почти неограниченные возможности микроминиатюризации.
Недостатки дискретных систем
1. Необходимо применение АЦП и ЦАП.
2. Возникают специфические ошибки обработки сигналов, связанные с их дискретизацией и квантованием.

Достоинства и недостатки дискретных систем

xн(t + kТп) ( Тп – период сигнала, k = 0, 1, 2 …),
 спектр может быть определен путем разложения в ряд Фурье
где ω0 = 2π/Тп – основная частота преобразования. Множество значений {Сk} – амплитудный спектр, множество {ϕk} – фазовый. Спектр периодического сигнала – линейчатый с дискретом ω0.

сигналом xн(t) прямым (ППФ) и обратным (ОПФ) преобразованиями Фурье:
Амплитдный спектр: Аx(ω) =|Xн(jω)|,
Фазовый спектр: ϕx(ω) = arg(Xн(jω)) = ∠ Xн(jω).

Спектр непрерывного непериодического сигнала

(1)

сигнала, полученный преобразованием Фурье,
x(nT) – дискретизированный сигнал,
X(e jωT) – спектр дискретизированного сигнала, полученный преобразованием Фурье.

6

Фурье дискретизированного сигнала

Отсчеты дискретизированного сигнала совпадают по величине с соответствующими значениями непрерывного сигнала (квантование по уровню в данном случае не рассматривается).
x(nT) = xн(t=nT) = xн(nT). (3)

размером ωд:

x(nT) = xн(t=nT) = xн(nT). (3)

(1)

(2)

точностью до множителя 1/Т из суммы бесконечного числа спектров непрерывного сигнала, следующих с периодом ωд.
Графики, приведенные на рис. 8, иллюстрируют, как формируется спектр сигнала при низкой (ωд < 2ωв) и высокой (ωд ≥ 2ωв) частотах дискретизации, где ωв – верхняя частота спектра сигнала. Как видно из графиков, при ωд < 2ωв происходит наложение спектров.
необходимо выбирать частоту дискретизации ωд из соотношения
ωд ≥ 2ωв .
В этом случае при |ω| < ωд/2 |X(e jωT)| = (1/T) |X(jω)|, - тогда спектры дискретизированного сигнала и сигнала непрерывного времени будут совпадать.

2ωв, эффект наложения спектров отсутствует, и можно восстановить сигнал непрерывного времени по его дискретным отсчетам: x(nT) → xн(t).
|ω| < ωд/2 ⇒ X(e jωT) = (1/T) Xн(jω)
⇒
По формуле прямого дискретного преобразования Фурье
⇒

отсчетов дискретного сигнала x(nT). Рассмотрим функцию ϕn(t):
T/2π = 1/ωд и по формуле Эйлера e ± j z = cos(z) ± j·sin(z)
⇒
⇒

при ⇒ ϕ0(t) = 0 при
n = 1 соответствует функции ϕ0(t), задержанной на Т, и т.д.:

Рисунок 8 – Cуммирование функций ϕn(t) c весами x(nT)

частот (ФНЧ) (см. на рис. 9 и рис. 10) с частотой среза ωc = ωд/2 и с АЧХ А(ω) = Т в полосе пропускания:

Рисунок 9 – Восстановление непрерывного сигнала из дискретизированного

Рисунок 10 – Восстановление непрерывного сигнала из дискретизированного: вверху – периодический спектр дискретизированного сигнала, посредине – амплитудно-частотная характеристика ФНЧ, внизу – результат применения ФНЧ к дискретизированному сигналу

при Т = const, – ∞ < n < ∞

Виды дискретных сигналов

Рисунок 12 – Произвольный сигнал, задержанный на m тактов, здесь m=1.

Рисунок 13 – Единичный импульс

14 – Аналоговый и дискретизированный гармонический сигнал

Непрерывный гармонический сигнал имеет вид
xн(t) = A sin(ωt) = A sin(2πf t),
где f – циклическая частота. После дискретизации с интервалом
Т = 1/fд = 2π /ωд он описывается формулой x(nT) = A sin(2πf nT).
Введем понятие относительной (нормированной) частоты :
α = f / fд = ω /ωд = fT. Тогда x(n) = A sin(2παn).

сигналами

В комплексной форме гармонический непрерывный сигнал имеет вид:
x(t) = A Re {e jωt} = A Re{cos(ωt) + j sin(ωt)} = A/2 Re{e jωt + e -jωt}.
После дискретизации сигнал описывается формулой
x(nT) = A Re{e jωnT} = A Re{cos(ωnT) + j sin(ωnT)} = A/2 Re{e jωnT + e -jωnT}.

виде линейной суммы соответствующих откликов системы:
x(n) = a1x1(n) + a2x2(n) → y(n) = a1y1(n) + a2y2(n)

Инвариантность в сдвигу

Задерживание сигнала на входе системы на m тактов вызывает задержку ее отклика на этот сигнал на m тактов:
x(n-m) → y(n-m)

Импульсная характеристика

Импульсная характеристика (ИХ) h(n) – это реакция системы на единичный импульс x0(n) :
x0(n) → y(n) = h(n)
Системы с конечной ИХ: h(n) = 0 при n < n1 и n > n2 (n1 < n2)
и бесконечной ИХ: h(n) = 0 при n < n1

можно представить в виде взвешенной суммы задержанных единичных импульсов:

По определению реакцией системы на единичный импульс x0(n) является ее импульсная характеристика h(n). Тогда реакцию y(n) системы на входной сигнал x(n) можно представить так:

Это и есть формула свертки.
Формула свертки записывается как y(n)=x(n)*h(n)=h(n)*x(n),
т.е. формула свертки симметрична.
В аналоговой записи формула свертки соответствует интегралу Дюамеля:

систем

При последовательном соединении систем ИХ h(n) эквивалентной системы определяется как свертка ИХ h1(n) и h2(n) подсистем: h(n) = h1(n) * h2(n)
При параллельном соединении систем h(n) = h1(n) + h2(n)

вызывает ограниченный по величине сигнал на выходе, т.е. при |x(n)|<∞ на входе |y(n)|<∞ на выходе.
Необходимое и достаточное условие устойчивости системы:
Докажем от противного: пусть
Сформируем сигнал на входе системы так (см. рис. 18):
Тогда сигнал на выходе при n = n0 будет
Значит, необходимость выполнения этого условия доказана.

Рисунок 18 – Формирование входного сигнала системы на основе формы ИХ

входного сигнала |x(n)| ≤ M. Тогда на выходе
что и требовалось доказать!

Физическая реализуемость системы

Система физически реализуема, если сигнал y(n0) на ее выходе на временном такте n0 зависит от текущего и ему предшествующих входных сигналов x(n ≤ n0) и не зависит от последующих входных сигналов x(n > n0). Т.е. следствие не может опережать причину.
Так как единичный импульс x0(n < 0) = 0, то и ИХ h(n<0) = 0.

сигнал.
Пусть на входе линейной системы гармонический сигнал частоты ω амплитудой Ax дискретизированный c интервалом Т:
xω(nT) = Ax(ω) e jωnT,
тогда сигнал на выходе может отличаться от входного амплитудой Ay и фазовым сдвигом ϕ:
yω(nT) = Ay(ω) e j[ωnT+ϕ(ω)].
Для конкретной заданной частоты ЧХ H(e jωT) связывает значения выходного и входного сигналов:
yω(nT) = xω(nT) H(e jωT).
⇒ H(e jωT) = yω(nT) / xω(nT) = [Ay(ω) / Ax(ω)] e jϕ(ω).

/ Ax(ω)] e jϕ(ω)
включает амплитудночастотную характеристику (АЧХ), которая определяется как модуль ЧХ и характеризует коэффициент усиления системы на различных частотах:
K(ω) = |H(e jωT)| = Ay(ω) / Ax(ω), 0 ≤ K(ω) < ∞;
и фазочастотную характеристику (ФЧХ), характеризующую фазовый сдвиг выходного сигнала по отношению к входному:
ϕ(ω) = arg [H(e jωT)] = ∠ H(e jωT), –π < ϕ(ω) < π .

= e jωnT (Ax(ω) = 1).
Сигнал y(nT) на выходе системы определяется с помощью свертки входного сигнала с ее ИХ h(mT):
Но по определению частотной характеристики
.
Сравнивая последние две формулы получаем
т.е. ЧХ системы вычисляется как преобразование Фурье от ее ИХ.

через их Фурье-спектры X(e jωT) и Y(e jωT) :
Сигнал на выходе системы можно также определить, если каждую составляющую спектра входного сигнала умножить на соответствующее значение ЧХ:
Из сравнения двух последних формул следует, что
⇒
ЧХ вычисляется как отношение спектров сигналов на входе и выходе

функция с периодом, равным частоте дискретизации ωд.
Так как спектры дискретизированных входного и выходного сигналов – периодические функции с периодом ωд, то аналогичным свойством обладает и ЧХ. Значит, ЧХ достаточно определить на интервале частот [0, ωд] (для нормированных частот – на интервале [0, 1]).
3. Если ИХ h(mT) системы – вещественная функция, то АЧХ К(ω) – __четная функция, а ФЧХ ϕ (ω) – нечетная. Докажем это.
Положим ω = ω1, тогда
Пусть теперь ω = ω2 = –ω1, тогда
Функции и – комплексно сопряженные. При смене знака аргумента модуль не меняет знак а фаза - меняет. Следовательно, АЧХ четная функция, а ФЧХ - нечетная, ч.т.д. Их достаточно вычислять на интервале частот [0, ωд/2]

переменная Лапласа, С – контур, включающий все особенности X(s).

Аналогом преобразования Лапласа при работе с системами дискретного времени является Z-преобразование. Z-преобразование последовательности x(n) описывается выражением
Например, для x(n) =[ 2 –3,5 1,4 –0,6….] X(z) = [2 –3,5z –1 1,4z –2 –0,6z -3….]. Существует и обратное Z-преобразование, позволяющее по X(z) определить x(n).

то
Z[а1x1(n) + а2x2(n)] = а1 X1(z) + а2 X2(z) (это правило распространяется на любое число слагаемых).
2. Z-преобразование задержанной на m тактов последовательности
Если x(n) = x1(n–m), то X(z) = z –m X1(z).

Передаточная функция

Передаточная функция H(z) системы есть Z-преобразование ее ИХ:
Если y(n), x(n) и h(n) связаны между собой формулой свертки
то
Значит, передаточную функцию определяется как отношение Z-преобразований сигналов на выходе и входе: H(z) = Y(z) / X(z) .

из передаточной функции H(z) путем подстановки z = e jω:

дифференциальными уравнениями связывающим входной и выходной сигналы x(n) и y(n) :
или ,
где N – порядок уравнения, ak и bk – постоянные коэффициенты. Вспомним формулу вычисления производной :
Допустим, переходя к случаю дискретизированного времени, что интервал дискретизации Δ t равен 1, тогда в качестве производной можно использовать разность двух соседних отсчетов сигнала, и от дифференциальных уравнений переходят к разностным (в общем случае количества слагаемых в первой и второй суммах могут различаться):

уравнению

Порядок N ЦФ определяется максимальным числом элементов задержки в прямом или обратном регистрах. Если все коэффициенты в регистре обратной связи равны нулю (ak = 0 для k = 1…N), то такой ЦФ называется нерекурсивным (регистр обратной связи отсутствует). Если хотя бы один коэффициент аk не равен нулю, то ЦФ является рекурсивным.

разностным уравнением
y(n) = x(n) + a y(n–1), (N = 1, b0 =1, b1 = 0, a1 = –a).

Рисунок 20 – Пример цифрового фильтрр 1-го порядка

1. Свойство линейности. Так как схема ЦФ состоит из линейных элементов, то для нее справедлив принцип суперпозиции: если на входе фильтра x(n) = c1·x1(n) + ... + cp·xp(n), то на его выходе y(n) = c1·y1(n) + ... + cp·yp(n).
2. Система инвариантна к сдвигу (b0 = const, a1 = const), то есть это система с постоянными во времени параметрами.
3. ИХ h(n) фильтра равна:
следовательно, это БИХ-фильтр.

устойчивости
Если 0 < a < 1, то по определению суммы геометрической прогрессии сумма S модулей отсчетов ИХ (см. рис. 21) определяется как

Рисунок 21 – Импульсный отклик устойчивой системы

В этом случае система устойчива.

геометрической прогрессии сумма S модулей отсчетов ИХ определяется как
Поэтому в этом случае система неустойчива. Ее отклик на единичный импульс показан на рис 22.

Рисунок 22 – Импульсный отклик неустойчивой системы

6. Система физически реализуема, т.к. h(n < 0) = 0.
7. ЧХ системы описывается выражением

–ФЧХ исследуемого фильтра ϕ(ω) = arg[H(e jω)] для угловых ω, циклических f, и нормированных α частот

ЧХ – непрерывная функция частоты ω.
ЧХ – периодическая функция с периодом, равным ωд.
К(ω) – четная функция частоты (К(ω1) = К(–ω1)).
ϕ(ω) – нечетная (ϕ(ω1) = – ϕ(–ω1)).
8. ИХ может быть вещественной только при вещественных ak и bk.
9. Передаточная функция:

Пример цифрового фильтра

⇒ .
Выполнив Z-преобразование от обеих частей разностного уравнения, получим
⇒

Связь передаточной функции с разностным уравнением

на множители:

Передаточная функция системы:

Здесь K – коэффициент усиления, вещественная величина,
zк – нули ПФ (zero), pk – полюса ПФ (pole), – вещественные или комплексно сопряженные величины. Таким образом, система задается множествами К, {zk}, {pk}.

справедливо: |pk| < 1. Для устойчивой системы полюса ПФ должны лежать внутри окружности единичного радиуса (r = 1). См. примеры, приведенные на рис. 25.

Рисунок 25 – Слева – устойчивая система, справа - неустойчивая

x(n) + a y(n-1) .
Передаточная функция этой системы имеет вид:

z1 = 0. Если подобрать такое значение параметра a, что p1′= 0,5, система будет устойчива, если p1″= 1,5 – система неустойчива (см. на рис. 26).

Рисунок 26 – Положение полюсов ПФ для устойчивой (p1′) и неустойчивой (p1″) системы

N + 1, причем bk = h(k), k = 0 … N, где N и bk – соответственно порядок фильтра и коэффициенты разностного уравнения.
Рекурсивные фильтры (РФ), как правило, обладают бесконечной ИХ, но могут иметь и конечную. Например, фильтру на рис. 27 соответствует уравнение
y(n) = x(n) – x(n–2) + y(n–1) (b0 = 1, b1 = 0, b2 = –1, a1 = –1),
но этот РФ имеет ИХ длины 2, как это проиллюстрировано рисунком 28.

Рисунок 27 – Рекурсивный фильтр 2-го порядка

Рисунок 28 – Входной (сверху) и выходной (снизу) сигналы этого ЦФ

ЦФ можно описать:
Разностным уравнением:
y(n) = f [ bk, ak, x(n–m), y(n–m) ].
Сверткой:
Частотной характеристикой: x(n) → ППФ → X(e jω),
Y(e jω) = X(e jω)H(e jω),
Y(e jω) → ОПФ → y(n).
Здесь используются прямое (ППФ) и обратное (ОПФ) преобразования Фурье.
Передаточной функцией:
x(n) → Z-преобразование → X(z),
Y(z) = X(z)H(z),
Y(z) → обр. Z-преобр. → y(n).

единичный импульс:
x(n) = x0(n)→ разностн. уравнение → y(n) = h(n).
6. Передаточной функции с разностным уравнением:
H(z) =
= Z{f [bk, ak, x(n–m), y(n–m)]}.
7. Импульсной и частотной характеристик:
h(n) → ППФ → H(e jω).
8. Импульсной характеристики и передаточной функции:
h(n)→ Z-преобразование → H(z).
9. Передаточной функции и частотной характеристики:
H(z = e jω) = H(e jω).

Соотношение параметров, характеризующих ЦФ

Рисунок 29 – Соотношение параметров ЦФ

полосы пропускания: а – ФНЧ, б – ФВЧ, в – ПФ, г – РФ

По параметрам полосы пропускания ЦФ делятся на следующие типы:
1. Фильтры низких частот (ФНЧ) – пропускают частоты ω < ω0, где ω0 – частота среза (см. на рис. 30.а).
2. Фильтры высоких частот (ФВЧ) – пропускают частоты ω > ω0. (рис. 30.б).
3. Полосовые фильтры (ПФ) – пропускают частоты ω01 < ω < ω02 (ω01, ω02 – частоты среза) (см. на рис. 30.в).
4. Режекторные фильтры (РФ) – пропускают частоты ω < ω01 и ω > ω02 (рис. 30.г).
5. Другие ЦФ, например, имеющие несколько полос пропускания или подавления.

Реальные АЧХ отличаются от приведенных идеальных:
1. Переходные полосы не имеют нулевой ширины.
2. АЧХ в полосах пропускания и подавления может не быть строго постоянной.

пропускания и задерживания изменяется монотонно. (см. рис. 31 и 32)

На рис. 31 и рис. 32 использованы следующие обозначения:
ωp – граница полосы пропускания,
ωp = {ωp(1), ωp(2)} – границы полосы пропускания,
ωs – граница полосы подавления,
ωs = {ωs(1), ωs(2)} – границы полос подавления,
Rp – уровень пульсаций в полосе пропускания в дБ,
Rs – минимальное затухание в полосе задерживания в дБ.

полосе пропускания, крутизна АЧХ в переходной полосе выше, чем у фильтра Баттерворта (см. рис. 33).
Чебышева 2 рода – имеются пульсации в полосе задерживания, крутизна АЧХ в переходной полосе выше, чем у фильтра Баттерворта (см. рис. 34).
Эллиптический (Кауэра) – имеются пульсации в полосах пропускания и задерживания, крутизна АЧХ в переходной полосе наиболее высокая (см. рис. 35).

ЦФ

Прямая форма непосредственно соответствует разностному уравнению
При реализации этой схемы требуется от N до 2N элементов задержки (N – порядок ЦФ).

функции в виде

Это соответствует следующей укрупненной вычислительной схеме (см. рис. 37)

Рисунок 37 – Представление ЦФ в виде пары последовательно соединенных фильтров

выходного y(n) и «промежуточного» w(n) сигналов:

схема, представленная на рис. 38

Рисунок 39 – Каноническая форма представления ЦФ

Получаем каноническую форму реализации ЦФ, используя вместо двух регистров сдвига один (см. на рис. 39).
Таким образом можно сократить число элементов задержки до двух раз.

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

Рисунок 40 – «Промежуточная» вычислительная схема ЦФ

– «Промежуточная» вычислительная схема ЦФ

Преобразуем схему, представленную на рис. 40:
поменяем местами схемы сложения и задержки сигналов,
объединим задержки в соответствующих каналах,
будем использовать только схемы задержки на 1 такт.

имеющих непрерывный по частоте Фурье-спектр. Рассмотрим теперь случаи, когда и частота имеет дискретную форму.
Сигнал x(nT), дискретизированный с частотой ωд = 2π/T, при –∞ < n < ∞ связан со своим спектром X(e jωT) прямым и обратным преобразованиями Фурье:
Спектр X(e jωT) – это непрерывная периодическая (с периодом ωд) функция частоты.
Рассмотрим функции:
бесконечную во времени периодическую: xp(nT) c периодом Тп = NT,
xр(nT) = xр[(n+kN)T] при k = 0, 1, 2 … ;
конечную: x(nT) = 0 при n<0 и n >N–1, имеющую длительность NT = Tп.

Фурье (ПДПФ и ОДПФ):
k – отсчеты частоты
n – отсчеты времени
Они получены из непрерывных преобразований путем дискретизации оси частот с шагом Ω = 2π/Tп = 2π/(NT) = ωд/N.
X(kΩ) представляет собой дискретную периодическую (с периодом NΩ = ωд) функцию, состоящую из отсчетов непрерывного спектра X(e jωT), взятых с интервалом Ω :
X(kΩ) = X(e j kΩ nT).
Если Т = const, то ωд = const и Ω = const. Поскольку Ω = 2π ⁄ NT, запишем:

для прямого и обратного ДПФ:

с соответствующими значениями преобразования Фурье и следуют с шагом Ω = ωд/N по частоте.
ДПФ – периодическая функция с периодом в N отсчетов:
X(k) = X(k + mN), m = 0, 1, 2 …
3. Свойство линейности:
при любом числе слагаемых.
4. ДПФ задержанной на m тактов последовательности:
__Для конечных последовательностей осуществляется циклический сдвиг.
5. Свойство симметрии для действительных последовательностей:
X(k) = X(N–k).

__преобразование Фурье по результатам вычисления z-преобразования.
7. Квадратичная зависимость числа умножений от длины __последовательности N является основным недостатком ДПФ. Требуется __N2 комплексных умножений (4N2 умножений действительных чисел) и __N(N–1) комплексных сложений (2N(N–1) действительных).
8. Размер оперативной памяти: для записи N отсчетов входной __последовательности и N отсчетов ДПФ необходимо 2N ячеек для хранения __комплексных чисел (4N – для хранения действительных чисел).

Свойства дискретного преобразования Фурье

применив ДПФ один раз.
Пусть x(n) и y(n) – действительные последовательности (n = 0, … , N−1). Последовательности разной длины, их можно выровнять, дополнив нулями.
Определим спектры X(k) и Y(k) последовательностей:
1. Сформируем последовательность v(n) = x(n) + j y(n) и вычислим её ДПФ:

Определим значения вещественных и мнимых частей V(k) и V(N–k) учитывая, что для ДПФ действительных последовательностей x(n) и y(n) справедливы следующие соотношения:
Re [X(k)] = Re [X(N–k)],
Re [Y(k)] = Re [Y(N–k)],
Im [X(k)] = – Im [X(N–k)],
Im [Y(k)] = – Im [Y(N–k)].

Re [X(k)] – Im [Y(k)], (1)
Im [V(k)] = Im [X(k)] + Re [Y(k)], (2)
Re [V(N–k)] = Re [X(N–k)] – Im [Y(N–k)] = Re [X(k)] + Im [Y(k)], (3)
Im [V(N–k)] = Im [X(N–k)] + Re [Y(N–k)] = – Im [X(k)] + Re [Y(k)]. (4)
Попарно складывая или вычитая равенства (1) – (4), получим значения вещественных и мнимых частей X(k) и Y(k)
(1) + (3) → Re [X(k)] = ½ { Re [V(k)] + Re [V(N–k)]},
(2) – (4) → Im [X(k)] = ½ {Im [V(k)] – Im [V(N–k)]},
(2) + (4) → Re [Y(k)] = ½ {Im [V(k)] + Im [V(N–k)]},
(3) – (1) → Im [Y(k)] = ½ {Re [V(N–k)] – Re [V(k)]}.
Таким образом вычислено ДПФ двух действительных последовательностей путем выполнения единственной операции ДПФ.

ДПФ X(k) последовательности x(n) для k = 0, … , N–1.
Формула ОДПФ имеет вид
Обозначим как x*(n) и X*(k) функции, комплексно-сопряженные с x(n) и X(k). Тогда
Умножение левой и правой части равенства на N дает
где X(n) – ДПФ функции X*(k).
Следовательно, ОДПФ вычисляется с помощью прямого ДПФ следующим образом:
X(k) →*→ X*(k) →ПДПФ→ X(n) →/N→ x*(n) →*→ x(n).

n = 0, …, N–1, N = 2L, L >0 – целое число. Определим её ДПФ X(k), k = 0, … , N–1.
Разобьем x(n) на две подпоследовательности (см. на рис. 42), в одну из которых попадут четные члены (x1(n)), а в другую – нечетные (x2(n)):
x1(n) = x(2n), n = 0, … , N/2–1.
x2(n) = x(2n+1), n = 0, … , N/2–1.

Рисунок 42 – Деление входной последовательности (сверху) на подпоследовательности четных (посредине) и нечетных (снизу) отсчетов

подставим результат преобразования в формулу для X(k):

0, …, N–1, но X1(k) и X2(k) определены для k = 0, …, N/2–1. Будем вычислять X (k) так:

Выполним во второй формуле следующее преобразование:

⇒

⇒

Вычисляем ДПФ последовательности так:
x(n) → {x1(n), x2(n)} → {X1(k), X2(k)} → X(k).

⇒

Число комплексных умножений вместо N2 равно 2 (N/2)2+N/2=N/2 (N+1)≈N2/2.

из одних и тех же элементов:

и так далее.

Для получения результатов каждой пары вычислений может быть использована единая базовая функция, показанная на рис. 43 и рис. 44:

Рисунок 43 – Вычислительный граф для единой базовой вычислительной функции ДПФ

Рисунок 44 – Вычислительный граф «бабочка» для БПФ с прореживанием по времени

осуществив аналогичное разбиение подпоследовательностей x1(n) и x2(n) прореживанием по времени.
Пусть, например, последовательность x(n) имеет длину N = 23 = 8. Тогда окончательный алгоритм формируется в 3 этапа (L = 3).

Рисунок 45 – Вычислительный граф, полученный после 1-го этапа преобразований

1-й этап. Полученный вычислительный граф можно представить в виде, показанном на рис. 45. Требуется 2×16+4=36<64=N2 комплексных умножений.

x2(n) длины N/2 = 4 делятся соответственно на подпоследовательности a(n), b(n), c(n), d(n) длины N/4=2:
x1(n) → a(n), b(n); x2(n) → c(n), d(n).
⇒
где A(k), B(k) – ДПФ от a(n), b(n). Аналогично преобразуем X2(k). Тогда вычислительный граф ДПФ на 2-м этапе принимает вид, приведенный на рис. 46.

Рисунок 46 – Вычислительный граф ДПФ, полученный на 2-м этапе преобразований

ДПФ можно выполнитьс помощью базовой операции БПФ. На рис. 47 представлен граф двухточечного ДПФ, полученный с учетом того, что

Рисунок 47 – Вычислительный граф двухточечного ДПФ

Рисунок 48 – Вычислительный граф базовой двухточечного операции БПФ

Тогда граф базовой операции БПФ для 2-точечной последовательности можно представить в виде, показанном на рис. 48.

8-точечного БПФ с прореживанием по времени

n = 0, …, N–1, N = 2L, L >0 – целое число. Определим её ДПФ X(k), k = 0, … , N–1.
Разобьем x(n) на две подпоследовательности (см. на рис. 50), в одну из которых попадает первая половина членов x(n), а в другую – вторая:

Рисунок 50 – Деление входной последовательности (вверху) на подпоследовательности первой половины (посредине) и второй половины (внизу) отсчетов

имеет вид:

Далее выполняются преобразования, аналогичные разделу, посвященному БПФ с прореживанием во времени.
Основное отличие:
- определяются по разным формулам четные и нечетные отсчеты БПФ;
- и в качестве базовой вычислительной операции используется граф “бабочка” показанный на рис. 51.

Рисунок 51 – Вычислительный граф «бабочка» для БПФ с прореживанием по частоте

8-точечного БПФ с прореживанием по частоте

Как и при применении прореживания по времени, количество комплексных умножений М = 12, или при исключении умножений на W80 = 1 число комплексных умножений М1 = 5.

Проводя соответствующие преобразования для последовательности длины N = 8, получим результирующий граф вычисления БПФ, показанный на рис. 52.

Быстрое преобразование Фурье с прореживанием по частоте

можно еще уменьшить, принимая во внимание, что

Для получения коэффициентов , можно применять:
Табличный метод – требует N/2 комплексных ячеек памяти.
Программный метод – требует использования подпрограмм sin(x) и cos(x), __что замедлит работу.
Рекуррентный метод: основан на использовании соотношения
при Его применение возможно, т.к. на каждом этапе показатель степени m меняется с одним и тем же шагом. Например, для N = 8 коэффициенты на
трех этапах имеют вид и ,
соответственно. Необходимо запомнить только шаги и .

Оба алгоритма (с прореживанием по времени и по частоте) выполняют дискретное преобразование Фурье числовой последовательности :
x(n), n = 0, …, N–1; → X(k), k = 0, …, N–1; N = 2L;
где L – целое положительное число. Если это условие не выполняется, то можно дополнить последовательность нулями до длины, равной ближайшему N = 2L.

граф 8-точечного БПФ с прореживанием по частоте

Рисунок 49 – Вычислительный граф 8-точечного БПФ с прореживанием по времени

Меняется порядок следования входных данных

Меняется порядок следования входных данных

выполнена с помощью двоичной инверсии. Если номер элемента массива до пересортировки записывается L-разрядным двоичным числом aL-1 ... a1 a0 , то его двоичная инверсия имеет вид a0 a1 ... aL-1 . Она и определяет номер элемента массива данных после
пересортировки, как это показано на рис 53.

Рисунок 53 – Пример изменения (пересортировки) номеров элементов последовательности методом двоичной инверсии

Например, элемент x(6) массива после пересортировки окажется на месте элемента x(3):
x(610) = x(1102) → x(0112) = x(310).

последовательности N=2L, где L – целое число. Если это не выполнялось, последовательность дополнялась нулями до ближайшей длины N=2L . Основой подхода была возможность представления одномерной последовательности в виде многомерного массива, например:

0, 1, 2, 3, 4, 5, 6, 7

Длина: N = 8 = 2×2 ×2

2

2

2

Если длина последовательности простое число, представить совокупность ее данных в виде многомерного массива не получится. Непростые числа можно разложить на множители, например, 60 = 4×3×5, тогда одномерную последовательность можно представить как многомерную и выполнить БПФ.

4

3

5

= 12×5 = 4×3×5, разложена на строки длиной 12, каждая строка также раскладывается на подстроки длиной 4.

также двумерную последовательность, строки которой нумеруются переменной s, а столбцы нумеруются переменной r. Тогда номер k текущего элемента массива вычислим как

Тогда от одномерного ДПФ (1) можно перейти к многомерному следующим образом:

(1)

(2)

формуле (2) и учтем при этом, что длина последовательности N=М×L, а
WMLlr = WNlr = 1, т.к. по определению:

(2)

.

.

⇒

,

где q(s, m) – L-точечное ДПФ m-того столбца исходного массива с ядром WM преобразования.
Теперь найдем новый массив h(s, m) = q(s, m) × Wms, где Wms – поворачивающий множитель. Тогда формула (3) примет вид:

(3)

.

.

Это M-точечное преобразование каждой из строк матрицы h(s, m) с ядром преобразования WL.

N2 = М2×L2 ) было заменено на (3) - последовательным применением ДПФ размерности M и размерности L и дополнительным N-кратным умножением на поворачивающий множитель (вычислительная сложность пропорциональна M2+L2+N << N2) .

Единый подход к алгоритмам БПФ

(1)

→

,

(3)

.

0, 1, 2, 3, 4, …, 59 →

суммирования.

(4)

Варианты (3) и (4) можно уподобить выполнению алгоритма Кули-Тьюки соответственно с прореживанием по времени или по частоте.

последовательности: N = L×M ×P

L

M

P

Рассмотрим случай представления одномерной последовательности трехмерной:

В этом случае ДПФ будет выполняться в виде следующей последовательности операций:

Вычислительная сложность:

– Пример выполнения БПФ с непостоянным основанием

с одинаковыми периодами, равными N. Свертка yp(n) этих последовательностей имеет вид:
yp(n) – периодическая последовательность с периодом, равным N.
Запишем выражения для ДПФ последовательностей-операндов свертки:
⇒
Это теорема о свертке: Фурье-спектр свертки равен произведению спектров ее операндов. Таким образом, алгоритм вычисления быстрой свертки имеет вид:
xp(n), hp(n)→БПФ→Xp(k), Hp(k)→×→Yp(k)→ОБПФ→yp(n).

доказана с применением преобразования Лапласа как теорема Бореля.
Для функций непрерывного времени свертка описывается интегралом Дюамеля:
Преобразование Лапласа от y(t) имеет вид:
Здесь был изменен порядок интегрирования, были сделаны замены переменных (t1 = t–s; t = t1+s; dt = dt1) и изменены пределы интегрирования (см. рисунок).
⇒

Рисунок – Изменение пределов интегрирования

…, Nx–1), h(n) – ИХ ЦФ (n = 0, …, Nh–1). Тогда сигнал на выходе ЦФ определяется с помощью прямой свертки:
Сигнал y(n), вычисленный с помощью свертки, имеет N = Nx + Nh – 1 ненулевых отсчетов.
Вычисление быстрой сверткой осуществляется в следующем порядке:
Последовательности x(n) и h(n) удлиняются до длины N путем добавления в каждую соответствующего числа нулей и рассматриваются как периодические xp(n) и hp(n) с периодом N; при выполнении БПФ длина N обеих последовательностей должна быть дополнительно увеличена до Ne – ближайшего значения, требуемого применяемым алгоритмом БПФ, например, при применении алгоритма Кули-Тьюки до Ne=2L, где L – целое число.
Для функций xp(n) и hp(n) вычисляются их БПФ Xp(k) и Hp(k).
Определяется произведение спектров:Yp(n) = Xp(k) Hp(k).
С помощью обратного ДПФ по Yp(n) находится функция yp(n). Подразумевается, что это один период длины Ne выходнго сигнала yp(n).

x длиной Nx и h длиной Nh требуется дополнить оба операнда до длины N=Nx+Nh-1 (для простоты опущено дополнение до длины 2L при использовании алгоритма Кули-Тьюки).
Вычислительная сложность выполнения БПФ каждого операнда и обратного БПФ будет N/2×log2 (N). Учитывая выполнение N операций умножения при перемножении спектров, получаем суммарную вычислительную сложность быстрой свертки:
CБ = 3[N/2 × log2 (N)]+N = 3[(Nx+Nh–1) /2×log2 (Nx+Nh–1)]+(Nx+Nh–1) .
При вычислении прямой свертки выполняется СП=(Nx+Nh-1)×Nh умножений (для простоты оценки считается что Nh< Однако при существенном различии длины операндов использование быстрой свертки может оказываться невыгодным. Например, при Nh=3, Nx >>Nh CП≈3Nx. Тогда CП

ИХ, вторая половина периода заполняется N нулями
Деление сигнала на “периоды” длиной 2N. Соседние периоды налагаются на N отсчетов. Первый период начинается с N нулевых отсчетов
Вычисление свертки каждого периода сигнала с расширенной ИХ, отбрасывание первых N отсчетов каждой свертки
Конкатенация результатов сверток в единый сигнал.

ограниченным числом разрядов.

Рисунок 66 – Ошибки квантования методом усечения

Рисунок 65 – Этапы обработки сигнала, на которых сказываются эффекты квантования

На различных этапах обработки входного сигнала xн(t) квантуются следующие сигналы и параметры (см. рис . 65:
входной сигнал xн(t);
коэффициенты ЦФ bk, ak;
результаты арифметических операций в ЦФ;
коэффициенты и результаты ДПФ и БПФ;
и др.
Рис. 66 иллюстрирует ошибку, возникающую при квантовании исходного сигнала, подлежащего обработке.

с шагом Q в b-разрядном АЦП (см. рис. 66):
Количество уровней квантования N = 2b;
xн(n) –точное значение непрерывного входного сигнала в момент времени t = nT;
x(n) – значение отсчетов квантованного сигнала.
e(n) = x(n) – xн(n) – ошибка квантования (см. рис. 67).

e(n) – дискретный стационарный случайный процесс,
e(n) и x(n) – не коррелированные случайные величины,
e(n) и e(m ≠ n) – не коррелированные случайные величины,
e(n) равномерно распределена в интервале Q = U / 2b, где U – динамический диапазон сигнала, в котором выполняется квантование.

Рисунок 67 – Формирование квантованного сигнала

– Формирование квантованного сигнала

Кm = mQ – уровни квантования,
m – целое число,
Q – шаг квантования

Метод квантования:
Km ≤ xн(n) < Km+Q ⇒ x(n) = Km

e(n) = x(n)–xн(n) – случайная величина, равномерно распределенная в интервале (–Q, 0];
m e = –Q/2 – математическое ожидание ошибки;
D e = Q 2/12 – дисперсия ошибки;
|e max| = Q – максимальное по модулю значение ошибки.

mQ – уровни квантования,
m – целое число,
Q – шаг квантования

Метод квантования:
Km –Q/2 ≤ xн(n) < Km+Q/2 ⇒ x(n) = Km

e(n) = x(n)–xн(n) – случайная величина, равномерно распределенная в интервале (–Q/2, Q/2];
m e = 0 – математическое ожидание ошибки;
D e = Q 2/12 – дисперсия ошибки;
|e max| = Q/2 – максимальное по модулю значение ошибки.

квантователем по методу округления c шагом квантования Q. Дисперсия ошибки квантования De = Q2/12. Отношение амплитуды сигнала к среднеквадратической ошибке квантования:
Необходимое число уровней квантования:
Тогда отношение сигнал/(шум квантования) имеет вид:
Учитывая, что N = 2b, предыдущую формулу можно представить в виде:

x(n) = xн(n)+e(n), где xн(n) – сигнал до квантования. Статистические характеристики шумов квантования e(n) при квантовании с шагом Q по методу округления имеют следующие значения:
m e = 0,
D e = Q 2/12,
|emax| = Q/2,
Для ЦФ с ИХ h(n) сигнал на выходе определяется с помощью свертки:
где yн(n) – выходной сигнал без квантования, евых(n) – шумы квантования на выходе.

Можно при известной ИХ h(n) и допустимой дисперсии Dвых шумов квантования на выходе определить необходимую разрядность b квантователя, учитывая, что Q = 2–b.

на выходе цифрового рекурсивного фильтра 1-го порядка, показанного на рис. 70. и описываемого следующим разностным уравнением:
y(n) = x(n) + a·y(n–1), a < 1.
Пусть De – дисперсия шумов квантования сигнала на входе ЦФ.

Рисунок 70 – Пример рекурсивного фильтра 1-го порядка

ИХ этого фильтра имеет вид:
Дисперсия шумов квантования на выходе ЦФ имеет вид:
При а → 1, De → ∞ .

ошибки квантования произведения

Рассмотрим шумы квантования, возникающие при умножении в различных схемах ЦФ. Пусть квантование производится по методу округления.

При умножения целых чисел одинаковой разрядности разрядность результата может удвоиться.
Результат yн(n) умножения приводится к допустимой разрядности мантиссы отбрасыванием ряда значащих разрядов.
Отбрасывание младших разрядов приводит к ошибкам, сходным с ошибкам квантования (см. рис. 71) и описывается выражением
y(n) = x(n)·a + e(n) = yн(n) + e(n).

Общий вид нерекурсивного фильтра

Рисунок 73 – Нерекурсивный фильтр с ошибками квантования при умножении

На рис. 73 ek(n) – ошибка квантования произведения в k-ом канале с математическим ожиданием me = 0 и дисперсией Dе = Q2/12.

Нерекурсивный фильтр с ошибками квантования при умножении

Рисунок 74 – Эквивалентная схема нерекурсивного ЦФ с ошибками квантования произведений

Эквивалентная схема нерекурсивного ЦФ с ошибками квантования произведений

– сумма ошибок в каналах, равная ошибке квантования на выходе ЦФ

– дисперсия ошибки на выходе ЦФ вычисляется суммированием дисперсий Dе в N каналах

Дисперсия Dвых шума квантования произведений на выходе фильтра зависит от порядка N фильтра и шага квантования Q и не зависит от коэффициентов bk.

Рекурсивный фильтр с ошибками квантования при умножении

каналах, D∑=(2N+1)Q2/12 – ее дисперсия;
, eвых(n)≠e∑(n), т.к. выходная ошибка зависит и от сигналов на предыдущих тактах времени;

зависит от N, Q, ak, hос.

Рисунок 76 – Эквивалентная схема рекурсивного фильтра с ошибками квантования произведений

Обозначим как hос ИХ системы обратной связи ЦФ. Тогда:

ДПФ, учитывающий ошибки квантования произведений, показан на рис. 77. Этот вычислительный граф может быть приведен к виду, представленному на рис. 78, причем суммарная по всем каналам ошибка e∑=eвых вычисляется по формуле:

Рисунок 77 – Вычислительный граф алгоритма ДПФ, учитывающий ошибки квантования произведений

Рисунок 78 – Приведенный вычислительный граф алгоритма ДПФ, учитывающий ошибки квантования произведений

ошибки ek в каналах также являются комплексными величинами:
ek = e′k + je″k .
Дисперсия каждой из составляющих e′k и e″k равна De = Q2/12, поэтому дисперсия комплексного числа еk равна D = 2De = Q2/6, а дисперсия ошибки квантования при вычислении каждого значения X(k) равна
DX = N·D = N·Q2/6.

Рисунок 77 – Вычислительный граф алгоритма ДПФ, учитывающий ошибки квантования произведений

базовой операции БПФ с прореживанием по времени

Рисунок 80 – Вычислительный граф базовой операции БПФ, учитывающий ошибки квантования произведений

Введем обозначения: W = W′+jW′′, b = b′+jb′′. Тогда
bW+e = b′W′+e1 + j(b′W′′+e2) + j(b′′W′+e3) – b′′W′′–e4 .
Ошибка на выходе равна
e = (e1 – e4) + j (e2 + e3).
Дисперсия каждой составляющей (действительной и мнимой) ошибки каждого из выходов базовой операции равна Dek = Q2/12. Дисперсия полной ошибки при выполнении каждой базовой операции равна De = 4Dek = Q2/3.

8-точечного БПФ с прореживанием по времени

Количество B базовых операций, выполняемых последовательно при расчете каждого отсчета БПФ (см. рис. 81, где N=2L=8, L=3), определяется как сумма членов геометрической прогрессии :

Для рассматриваемого примера В = 7. Таким образом, дисперсия ошибки, возникающей вследствие квантования результатов умножения при вычислении каждого отсчета БПФ, может быть рассчитана по формуле
D = De B = (N–1) Q2/3.

выходной сигнал при k0 = 0; снизу – выходной сигнал при k0 ≠ 0

Если ЦФ устойчив, при прекращении сигнала на входе выходной сигнал должен стремиться к 0, т. е. при x(n ≥ n0) = 0 y(n→∞) → 0. При квантовании чисел возможен эффект «предельных циклов»: y(n→∞) → k0 ≠ 0; k0 – мертвая зона.

y(n) = x(n)+a·y(n–1), а = 0,9 , n0 = 0 , y(–1) = 7. x(n ≥ 0) = 0
=> y(n ≥ 0) = a·y(n–1).
Для шага квантования Q = 1 расчет выходного сигнала приведен в таблице.
Таблица – Расчет сигнала на выходе фильтра
Определим величину k0 мертвой зоны. При квантовании округлением:
y(n) = k если a·k ≥ k – 0,5 ;
=> k (1 – а) ≤ 0,5 ;
=> k ≤ 0,5 / (1 – а) . (*)
Величина мертвой зоны находится как максимальное целое число k, удовлетворяющее условию (*). Например, для рассматриваемого фильтра:
а = 0,9 => k ≤ 0,5 / 0,1 = 5 => k0 = 5; => y(n) = 7, 6, 5, 5, 5 ;
а = 0,1 => k ≤ 0,5 / 0,9 = 0,55 => k0 = 0; => y(n) = 7, 1, 0, 0, 0 .

а < 0 значения k0 знакопеременны, как показано на рис. 83.

Рисунок 83 – Предельные циклы при a > 0 (слева) и a < 0 (справа).

и результат его дискретизации по времени и квантования по уровню, слева – плотность распределения вероятности значений аналогового сигнала

При равномерном квантовании усечением шум квантования не превышает шага квантования Q.
Неравномерное квантование учитывает статистические свойства сигнала, его применяют для минимизации среднеквадратической величины шума квантования. Шаги квантования меньше для более вероятных уровней сигнала (см. на рис. 84).

квантование потребует 12 двоичных разрядов (4096 уровней),
неравномерное – только 8 двоичных разрядов (256 уровней квантования).

Приравнивание к нулю частных производных этого выражения по ak и bk позволяет определить оптимальные параметры квантования:

Диапазон возможных значений сигнала x, описываемый плотностью распределения вероятностей p(x), делится на N неодинаковых зон квантования а0…а1, а1…а2, … аN-1...aN, каждой из которых ставится в соответствие квантованное значение bk ( bk∈[ak-1,ak] ).
Средний квадрат ошибки рассчитывается по формуле

обнаруживать на фоне шумов гармонического сигнала с неизвестными __амплитудой и частотой;
измерять эти его параметры с помощью ДПФ (БПФ).
Наблюдаем аналоговый сигнал xн(t) на интервале времени длиной Тs (0

Полезный сигнал us(t) гармонический с неизвестными амплитудой Us и частотой fs:
us(t) = Us cos(2π fs t) = Us cos(ωst).
Известно :
ω1 ≤ ωs ≤ ω2 .
Uш (t) - аддитивный гауссов шум с математическим ожиданием 0, и дисперсией __Dш = σш2 = Pш ≠ 0, где σш - СКО, Pш – мощность шума.
Дискретизируем сигнал xн(t) с интервалом Т = 2π /ωд (с частотой ωд = 2π fд) :
xн(t) →x(nT), n=0 … N-1,
N = [ Ts/T ], где [.] – результат округления.

- при отсутствии полезного сигнала это шум;
- или это смесь полезного сигнала и шума.

xн(t), полученный преобразованием Фурье;
X(e j ω t) – спектр дискретизированного сигнала x(nT), полученный с помощью _преобразования Фурье;
X(kΩ) – спектр дискретизированного сигнала x(nT), полученный с помощью _ДПФ, где Ω – основная частота ДПФ. Ω = 2π /(NT) = 2π /Ts = ωд /N.
Номера отсчетов дискретного спектра, соответствующих границам _анализируемого частотного диапазона, k1 и k2 (k1 = [ω1/Ω], k2 = [ω2/Ω]).

Спектральный анализ

…, k2.

Рисунок 85 – Спектр сигнала при σш = 0 : вверху – аналогового; посредине – дискретизированного; внизу –полученный с помощью ДПФ

в канале связи рассматриваются две статистические гипотезы :
|X(kΩ)| → ? → H0:(Us = 0) – гипотеза H0: сигнала нет;
|X(kΩ)| → ? → H1:(Us ≠ 0) – гипотеза H1: сигнал есть.
Если хотя бы один отсчет спектра в анализируемом диапазоне частот k превышает по амплитуде |X(kΩ)| установленный порог Uп (см. также на рис. 86), принимается гипотеза H1, в противном случае принимается гипотеза H0 :
|X(kΩ)| > Uп → H1,
|X(kΩ)| < Uп → H0,
k ∈[k1; k2], (σш ≠0).

Рисунок 86 – Принятие решения о наличии полезного сигнала; здесь km – номер максимального спектрального отсчета:

– вероятность принятия гипотезы H1 при Us=0:
F = P{|X(kmΩ)| > Uп} при Us = 0;
вероятность D правильного обнаружения – принятия гипотезы H1 при Us ≠ 0:
D = P{|X(kmΩ)| > Uп} при Us ≠ 0.
В общем случае F + D ≠ 1, так как суммируются условные вероятности, взятые при разных условиях.
Выбор порога Uп определяется принятым критерием обнаружения. В подобных задачах часто используется критерий Неймана-Пирсона:
сначала фиксируется допустимая вероятность ложной тревоги (F=F0=const);
затем максимизируется вероятность правильного обнаружения (D → max).

порога:
F = func(Uп) при Us = 0 (как показано на рис. 87);
По полученной кривой определяется порог Uп0, обеспечивающий допустимое __значение F = F0 вероятности ложной тревоги.
Находится зависимость D = func(Us) вероятности правильного обнаружения __от амплитуды сигнала при F = F0 (Uп = Uп0).
При заданной требуемой вероятности правильного обнаружения D = D0 (см. __на рис. 88) определяют необходимое отношение сигнал/шум q = Us0 /σш, __обеспечивающее заданные параметры F0 и D0. Обычно задается D0 ≈ 0,9. 

Рисунок 87 – Определение порогового уровня, соответствующего допустимой вероятности ложной тревоги

Рисунок 88 – Определение уровня сигнала, обеспечивающего нужную вероятность его обнаружения

ωs сигнала производится, только если __на предыдущем этапе была принята гипотеза H1.
Измеряемый сигнал – случайный процесс, следовательно, определяются не __сами измеряемые величины, а их оценки Us* и ωs* .
Вся информация, необходимая для проведения измерений, содержится в __положении km на частотной оси и величине |X(kmΩ)| максимального __спектрального отсчета в частотном диапазоне k1Ω ÷ k2Ω (см. на рис. 89).
Оценки частоты ωs* и амплитуды сигнала Us* находятся по ниже __приведенным формулам:
_________ω s* = km Ω ,
______Us*=2 |X(kmΩ)| / N .

Рисунок 89 – Измеряемые параметры сигнала

величина, равномерно __распределенная в интервале [–Ω/2, Ω/2]. Дисперсия ошибки D(δωs)=Ω 2/12.
Для уменьшения ошибки необходимо уменьшать основную частоту Ω ДПФ. __Т.к. Ω = 2π/Ts, это достигается удлинением интервала Ts наблюдения сигнала.
Ошибка измерения амплитуды δX = |X(jωs)| – |X(jωs*)| зависит от случайной _ величины – ошибки измерения частоты, ⇒ расчет дисперсии ошибки осложнен.
Ошибки измерения частоты и амплитуды не зависят от мощности шума в канале. __С ростом амплитуды сигнала суммарная дисперсия ошибки может стремиться __не к нулю, а к некоторому значению, определяемому методической ошибкой.

Если реальная частота сигнала ωs не попадает на дискретные точки ωk = k Ω оси частот, возникают методические ошибки измерения частоты и амплитуды, как показано на рис. 89 для σш = 0.

Спектральный анализ: методические ошибки измерения

Рисунок 89 – Методические ошибки изменения

характеристик обнаружения и измерений можно использовать метод Монте-Карло следующим образом в контексте данной задачи:
Проводится многократное формировании реализаций дискретизированного __случайного входного сигнала x(nT) (M испытаний) c постоянными значениями __параметров.
По каждой сформированной реализации принимается решение о наличии или __отсутствии полезного сигнала. Отношение числа К обнаружений к общему __числу M испытаний определяет частоту p* = K/M обнаружений при данных __значениях параметров.
Эта частота стремится при бесконечном возрастании М к вероятности p __обнаружения (M → ∞ => p* → p).
При наличии шумов величина и положение максимального отсчета ДПФ на частотной оси являются случайными величинами, и возможно определение только оценок Us* и ωs* амплитуды и частоты. Вследствие этого возникают:
ошибка измерения амплитуды – δUs = Us* – Us ,
ошибка измерения частоты – δω s = ωs* – ωs .

испытания, в котором обнаружен сигнала при данном наборе его параметров, определяются оценки амплитуды Usi* и частоты ωsi*.
Они позволяют вычислить следующие статистические характеристики:
математические ожидания оценок амплитуды и частоты
и
дисперсии ошибок
и
среднеквадратические ошибки
и
относительные среднеквадратические ошибки
γUs = σUs / Us и γωs = σωs / ω s .
Функциональные зависимости указанных характеристик от параметров сигнала и шума строятся путем изменения значения одного из параметров и повторения испытания. Наглядные зависимости получаются при числе испытаний М ≥ 1000.

данных (будем рассматривать изображения) и применяются для:
Извлечения полезных признаков анализируемых данных
Сжатия данных
Кодирования данных
Сокращения объема вычислений
…

Унитарные преобразования являются частным случаем линейных преобразований, когда:
линейный оператор обратим
ядро преобразования удовлетворяет условиям ортогональности

Матрица A унитарного преобразования обладает следующими свойствами:
A-1A = E,
где индекс A-1 – матрица, обратная матрице A; E – единичная матрица.

N1×N2 преобразуется в матрицу F (n1,n2) того же размера:

где AC (или BC) и AR (или BR) – соответственно одномерные операторы преобразования столбцов и строк.

в два этапа:
Сначала выполняется одномерное преобразование по всем столбцам матрицы изображения, и образуется матрица с элементами
Затем выполняется второе унитарное преобразование по всем строкам полученной матрицы:

Унитарные преобразования можно записывать с помощью векторных операций. Если F и f – соответственно матричное и векторное представление массива пикселов, а F и f – матричное и векторное представление результата преобразования A, тогда двумерное преобразование в векторной форме запишется так:
f=Af.

Представление двумерной матрицы пикселов в виде одномерного вектора путем растрового сканирования содержимого матрицы

номера {n1, n2} строк и столбцов матриц на пространственные координаты – номера {j, k} строк и столбцов пикселов.
Заменяем номера {m1, m2} __на номера частотных __отсчетов {u, v}

Подставив соответствующие комплексные экспоненты в матрицы преобразований получаем прямое и обратное дискретные преобразования Фурье
[здесь i = √(-1) ].

Ниже мы рассмотрим и другие полезные двумерные унитарные преобразования

пространственной частоты 2 по одной координате и частоты 3 – по другой

косинусов и синусов.

Получаем набор гармонических базисных функций разложения, как показано на рисунке.

Рисунок – Косинусные и синусные базисные функции дискретного преобразования Фурье последовательности длиной 16 отсчетов

⇒Значение нулевой гармоники – это N-кратно увеличенная средняя яркость изображения

постоянные, получим

т.е. вычисленный спектр периодический (см. на рисунке).

Рисунок – Справа – периодическое исходное изображение, слева – его периодический двумерный спектр, рассчитанный дискретным преобразованием Фурье

части спектра можно рассчитать по его ранее вычисленным значениям (см на рисунке)

Преобразование Фурье изображения

Рисунок – Часть спектральных отсчетов, для вычисления которых не требуется применение преобразования Фурье

центре, а ___ДПФ располагает их по углам спектральной плоскости (см. на рисунке). Чтобы ___разместить нулевую гармонику в центре «изображения» спектра, можно перед ___вычислением преобразования Фурье каждый пиксел преобразуемого ___изображения умножить на (–1)j+k, где jϵ[0; N–1], kϵ[0; M–1] – координаты ___пиксела, а N и M – количество пикселов по абсциссе и ординате.
ДПФ вычисляет действительную и мнимую чассти комплексного спектра, что ___увеличивает время вычислений и объем памяти, необходимый для хранения ___их результатов.
ДПФ медленно сходится из=за возникновения скачков яркости на границах ___«периодов» преобразуемых изображений.

Преобразование Фурье изображения

Рисунок – Спектральная плоскость: слева – естественное размещение начала координат, справа начало координат расположено в левом верхнем углу спектральной плоскости, рассчитанной с помощью ДПФ

___компоненты и мнимые синусные компоненты.
Для преобразуемых изображений, симметричных относительно начала координат в результате ДПФ должны остаться только спектральные компоненты, соответствующие четным (косинусным) базисным функциям.
Для создания симметричных преобразуемых изображений применяют два ___варианта отражения исходного изображения относительно «начала ___координат» как показано на рисунке.

Рисунок – два варианта построения симметричных изображений:
а – прикладывание друг к другу четырех зеркально отраженных изображений;
б – прикладывание зеркально отраженных изображений с наложением (в обоих случаях одноименные углы наложенных изображений помечены крестиками)

(j, k) номера строк j и столбцов k меняются в диапазоне [0; N–1], а в преобразованном (симметричном) изображении j, k∈[–N; N–1].

Массив Fs(j, k) симметричен относительно точки j=–1/2, k=–1/2. Поместив в точку симметрии начало координат, вычислим от него ДПФ:

В соответствии с рисунком пикселы симметричного изображение Fs(j, k) для симметричного четного ДКП формируется из исходного изображения F(j, k) согласно следующим формулам:

где пространственные частоты u и v меняются в диапазоне [–N; N–1].

из действительных чисел, следовательно, разложение комплексной экспоненты ДПФ по формулам Эйлера будет содержать только косинусные компоненты:

Четное симметричное ДКП вычисляют с дополнительной нормировкой:

где C(0) = 2-1/2, а C(w) = 1, при w = 1, 2, …, N–1.

C(w):

На рисунке слева приведен пример базисных функций, используемых в четном симметричном ДКП.

Рисунок – Базисные функции симметричного четного ДКП для преобразуемой последовательности , содержащей N = 16 отсчетов.

длиной 16 отсчетов

рисунком пикселы симметричного изображение Fs(j, k) для симметричного нечетного ДКП формируется из исходного изображения F(j, k) согласно следующим формулам:

Вычисление двумерного ДПФ для такого массива дает:

где u, v = –N+1, …, -1, 0, 1, …, N-1.

то для реальных изображений

Поэтому согласно определению нечетного симметричного ДКП для его вычисления используют следующее нормированное изображение:

производится по следующим формулам:

можно выполнить в виде последовательности одномерных.

k=0, …, N–1 – номера столбцов пикселов преобразуемого изображения F(j, k); u=0, …, N–1 – номера частотных отсчетов вдоль столбцов пикселов, а v=0, …, N–1 – номера частотных отсчетов вдоль строк пикселов преобразуемого изображения.
Обратное преобразование вычисляется по такой же формуле, только F(j, k) и (u, v) поменяются местами.

длиной 16 отсчетов

вектора
mk – математическое ожидание k-того компонента вектора
p – функция плотности распределения вероятности

Ковариационная матрица множества векторов

Собственный вектор матрицы

Собственным вектором невырожденной матрицы A размером N×N является такой вектор x, для которого выполняется следующее соотношение:
Ax = λx,
где λ – собственное число, соответствующее этому вектору. Невырожденная матрица A размером N×N имеет N линейно-независимых собственных векторов.

Рассчитаем для них ковариационную матрицу C и вычислим N ее собственных векторов yi и собственных чисел λi, i=1…N.
Упорядочим собственные векторы xi в порядке убывания величины соответствующих им собственных чисел λi . Выберем M собственных векторов, соответствующих M наибольшим собственным числам. Транспонируем выбранные собственные векторы и составим из них матрицу Y размером M×N:
где y`i = (yi1, yi2, …, yiN) – транспонированный i-тый собственный вектор.
Рассчитаем преобразованные векторы zi :
zi=Yxi,
это представление векторов xi в новом ортогональном базисе выбранных собственных векторов y`i, поскольку скалярное умножение вектора строки y`i на вектор столбец xi - это проекция вектора xi, на базисный вектор y`i .

представления векторов так, чтобы оси эллипсоида рассеяния векторов в исходном пространстве признаков совпали с направлением ортов {y1, y2, …, yM} в новом базисе.

Если разброс значений в направлении y2 пренебрежимо мал, можно сократить размерность пространства, оставив только орт y1, так можно сжимать данные.

отсчетов

в цифровой реализации
ортогональны
но не образуют полную систему

ортогональны
образуют полную систему

матрицу Адамара H, строки которой являются функциями Уолша (см. на рисунке):
F=Hf.

Рисунок – Матрица Адамара: а – представление маатрицы в виде изображения; b – графики функций Уолша, соответствующих строкам матрицы; Двоичное представление матрицы Адамара

Матрицы Адамара генерируются следующим образом:

,

,

HN-1=HN.

преобразование Адамара DWT по всем столбцам матрицы изображения B. Затем выполняется второе одномерное преобразование Адамара по всем строкам полученной матрицы:
F=DWT×B×DWT.
Аналогично выполняется обратное двумерное преобразование:
B=DWT×F×DWT.

Рисунок – Двумерное преобразование Адамара выполняет разложение изображения по таким двумерным функциям

Рисунок – Одна из аналогичных базисных функций преобразования Фурье

преобразование Адамара DWT по всем столбцам матрицы изображения B. Затем выполняется второе одномерное преобразование Адамара по всем строкам полученной матрицы:
F=DWT×B×DWT.
Эти сложные многомерные операции можно описать следующими формулами:
где j, k – соответственно номера строк и столбцов элементов двумерного преобразуемого массива данных (например, пикселов); u, v – соответственно номера базисных функций Адамара по строкам и столбцам преобразованного массива; p(j, k, u, v) вычисляется по следующему правилу:
где переменные ui, vi, ji и ki равны цифрам в двоичном представлении чисел u, v, j и k соответственно. Напимер, если u=13, то u3=1, u2=1, u1=0 и u0=1.

Для сравнения. Для преобразования Фурье:

отсчетов

Вычисление одномерного преобразования Хаара F выполняется умножением преобразуемого вектора b на матрицу преобразования H:
F = Hb; b=H-1F=HTF.

Рекурсивное правило генерации матриц Хаара:
генерация начинается с матриц H2 второго порядка.
Матрицы последующих порядков вычисляются по матрицам предыдущих порядков на основе произведения Кронекера:
,
где E2k – единичная матрица размерами 2k×2k, N=2k+1.

одномерное преобразование Хаара DHT по всем столбцам матрицы изображения B. Затем выполняется второе одномерное преобразование Хаара по всем строкам полученной матрицы:
F=DHT×B×DHT.
Вычисление двумерного преобразования Хаара двумерного массива можно представить как его разложение по двумерным базисным функциям, представленным на рисунке.

Рисунок – Двумерное преобразование Хаара выполняет разложение изображения по таким двумерным базисным функциям

из 16 отсчетов

Вычисление одномерного наклонного преобразования F выполняется умножением преобразуемого вектора b на матрицу преобразования S:
F = Sb.

Рекурсивное правило генерации матриц преобразования: генерация начинается с матриц S2 второго порядка.

Матрицы последующих порядков вычисляются по матрицам предыдущих порядков согласно формуле:

где Ik – единичная матрица k-того порядка,

длиной 16 отсчетов

Рассмотренные унитарные преобразования «концентрируют энергию» в начальных компонентах «спектра», что позволяет отбрасывать часть высокочастотных компонент. Осуществляя таким образом сжатие данных.