Файл: А. Б. Сергиенко минобрнауки россии санктПетербургский государственный электротехнический университет лэти им. В. И. Ульянова (Ленина) А. В. Петров а. Б. Сергиенко цифровая обработка сигналов лабораторный практикум.pdf

34 вопрос о том, как связаны друг с другом эти 2 спектральных представления, полученные на основе одних и тех же отсчетов сигнала.
Сравнение формул (1.3) и (3.2) показывает, что ДПФ представляет собой просто дискретные отсчеты спектральной функции дискретного сигнала, со- ответствующие частотам
2
n
n N
ω = π
ɶ
рад/отсчет:
2
( )
( )
n N
X n
X
ω= π
=
ω
ɶ
ɺ
ɺ ɶ
(3.3)
По этой причине значения ДПФ иногда называют
спектральными от-
счетами.
Формула (3.3) в сочетании с формулой связи абсолютных и нормиро- ванных частот (
д
(2 )
f
F
=
ω
π
ɶ
) определяет частотную шкалу ДПФ: спек- тральным отсчетам с номерами n = 0, …, N
−
1 соответствует сетка частот с шагом д
F
N , простирающаяся от нуля почти до частоты дискретизации
(рис. 3.1).
1
N
F
д
2
N
F
д
N
F
д
N-1
F
д
0 0
1 2
N-1
N
f
n
Рис. 3.1. Частотная шкала ДПФ
Из соотношения (3.3) следует еще один важный вывод: если добавить к конечному набору отсчетов некоторое количество нулей, спектральная функция дискретного сигнала (1.3), естественно, не изменится, но ДПФ (3.2) даст большее число спектральных отсчетов, соответствующих частотам, бо- лее тесно расположенным в интервале от нуля до частоты дискретизации.
Свойства ДПФ. В целом свойства ДПФ аналогичны свойствам непре- рывного преобразования Фурье, однако дискретный характер анализируемо- го сигнала привносит некоторую специфику.
Линейность. ДПФ, согласно определению (3.2), является линейной ком- бинацией отсчетов последовательности, поэтому оно подчиняется принципу суперпозиции:
1 2
( )
( )
( )
y k
ax k
bx k
=
+
↔
1 2
( )
( )
( )
Y n
aX n
bX
n
=
+
ɺ
ɺ
ɺ
Задержка.
При задержке сигнала на
∆k
отсчетов амплитудный спектр не меняется, а фазовый спектр приобретает дополнительное слагаемое, линейно зависящее от номера гармоники и от значения
∆k
:

35
( )
(
)
y k
x k
k
=
− ∆
↔
2
( )
( )
n
j
k
N
Y n
X n e
π
−
∆
=
ɺ
ɺ
Симметрия. ДПФ вещественного сигнала обладает свойством симмет- рии:
(
)
(
)
( )
X N
n
X
n
X
n
∗
−
=
− =
ɺ
ɺ
ɺ
Круговая свертка. Круговой свертке периодических последовательно- стей соответствует произведение их спектров:
1 1
2 0
( )
( )
((
) mod
)
N
n
y k
x n x
k
n
N
−
=
=
−
∑
↔
1 2
( )
( )
( )
Y n
X n X
n
=
ɺ
ɺ
ɺ
Равенство Парсеваля. Дискретный аналог равенства Парсеваля записы- вается следующим образом:
( )
1 1
2 2
0 0
1
( )
N
N
k
n
x
k
X n
N
−
−
=
=
=
∑
∑
ɺ
Матрица ДПФ. ДПФ является линейным преобразованием, трансфор- мирующим вектор временных отсчетов в вектор такой же длины, содержа- щий отсчеты спектральные. Такое преобразование может быть реализовано как умножение некоторой квадратной матрицы на входной вектор-столбец:
=
y
Ax
,
(3.4) где
A
— матрица преобразования. В случае ДПФ эта матрица имеет вид:
(
)
(
)
(
)
(
)
(
)
2 2
4 2
1 4
8 2
2 1
ДПФ
2 2
2 1
2 1
1 1
1 1
1 1
1 1
j
j
j
N
N
N
N
j
j
j
N
N
N
N
j
N
j
N
j
N
N
N
N
e
e
e
e
e
e
e
e
e
π
π
π
−
−
−
−
π
π
π
−
−
−
−
π
π
π
−
−
−
−
−
−












=












A
…
…
…
…
…
…
…
…
…
Общая формула для элемента матрицы, расположенного в n-м столбце
m-й строки, выглядит так:
(
)
(
)(
)
ДПФ
1 1
,
exp
2
, 1
, 1
m
n
m n
j
m
N
n
N
N
−
−


=
− π
≤
≤
≤ ≤




A

36
Вычисление ДПФ путем умножения матрицы на вектор полностью со- ответствует формуле (3.2). Этот метод требует большого количества вычис- лительных операций, поэтому на практике вместо него применяются быст- рые алгоритмы.
Алгоритм быстрого преобразования Фурье. Для вычисления одного коэффициента ДПФ по (3.2) необходимо выполнить N комплексных умноже- ний и сложений. Таким образом, расчет всего ДПФ, содержащего
N коэффициентов, потребует N
2
пар операций «умножение—сложение» (то же самое можно увидеть и из (3.4), представляющей ДПФ как умножение матрицы на вектор). Число операций возрастает пропорционально квадрату размерности ДПФ.
Однако, если N не является простым числом и может быть разложено на множители, процесс вычислений можно ускорить, разделив анализируемый набор отсчетов на части, вычислив их ДПФ и объединив результаты. Такие способы вычисления ДПФ называются быстрым преобразованием Фурье
(БПФ; английский термин — Fast Fourier Transform, FFT) и повсеместно ис- пользуются на практике. Например, при размерности ДПФ, равной степени двойки, широко используется алгоритм Кули—Тьюки (Cooley—Tukey), чис- ло операций в котором пропорционально N log
2
(N).
Растекание спектра. Как уже говорилось, при ДПФ предполагается, что последовательность отсчетов анализируемого сигнала является периодически продолженной вперед и назад во времени. При этом, если значения началь- ных и конечных отсчетов сигнала сильно различаются, при периодическом повторении на стыках сегментов возникают скачки, из-за которых спектр сигнала расширяется.
Это явление, называемое растеканием спектра (spectrum leakage), можно наглядно проиллюстрировать на простейшем примере вычисления спектра дискретного гармонического сигнала:
( )
cos(
),
0, 1,
,
1
x k
A
k
k
N
=
ω + ϕ
=
−
ɶ
…
(3.5)
Если анализируемая последовательность содержит целое число перио- дов гармонического сигнала (т. е. если отношение
(2 )
Nω
π
ɶ
является целым числом), то периодически продолженный сигнал представляет собой гармо- нические колебания (без скачков), а подстановка (3.5) в формулу ДПФ (3.2) показывает, что вычисленное ДПФ содержит лишь 2 спектральных отсчета, отличных от нуля:

37
,
,
2 2
( )
,
1
,
2 2
0,
в остальных случаях.
j
j
AN
e
n
N
AN
X n
e
n
N
ϕ
− ϕ
ω

=

π

ω
=



=
−



π




ɶ
ɺ
ɶ
Таким образом, аналогично спектру непрерывного гармонического сиг- нала, ДПФ отличается от нуля всего для двух значений
n
. Однако если отно- шение
(2 )
N
ω
π
ɶ
не является целым числом, спектр оказывается значительно более богатым. Этому можно дать простое объяснение: в данном случае пе- риодически продолженная последовательность уже не будет являться набо- ром отсчетов непрерывной синусоиды. Поэтому в полном соответствии со свойствами преобразования Фурье в спектре появляются дополнительные со- ставляющие.
−10 0
10 20 30
−1 0
1 0
5 10 15 0
5 10
−10 0
10 20 30
−1 0
1 0
5 10 15 0
5 10
Рис. 3.2. Происхождение растекания спектра: слева — периодически продолженные сигналы, справа — амплитудные спектры одиночных сигналов (штриховые линии) и модули ДПФ (кружки)
Необходимо подчеркнуть, что причиной растекания спектра является именно периодическое продолжение анализируемого сигнала. Спектр оди- ночного фрагмента дискретной синусоиды является периодической непре- рывной функцией частоты. Эта функция имеет лепестковую структуру неза-

38 висимо от того, целое или нецелое число периодов укладывается в анализи- руемом сегменте. Однако дискретный ряд частот (3.3), на которых вычисля- ется ДПФ, может быть по-разному расположен относительно лепестков спек- тральной функции. В случае целого числа периодов все анализируемые час- тоты (кроме двух) попадают как раз в нули между лепестками. При нецелом числе периодов такого не происходит.
Иллюстрация растекания спектра представлена на рис. 3.2, где построе- ны графики периодически продолженных сигналов, а также модули их ДПФ и непрерывных амплитудных спектров одиночных фрагментов синусоиды.
Верхняя пара графиков соответствует отсутствию растекания спектра
(на интервале анализа укладывается ровно 4 периода колебания). Нижние графики демонстрируют растекание спектра (в данном случае
(2 )
8 3
N
ω
π =
ɶ
не является целым числом).

39
Справка
Для вычисления ДПФ в MATLAB имеется функция fft (она вычисляет
ДПФ с использованием быстрых алгоритмов — FFT, Fast Fourier Transform).
В простейшем случае она вызывается без дополнительных параметров:
y = fft(x)
. Здесь x — вектор отсчетов сигнала; y — вектор результатов вычисления ДПФ.
Замечание
Чтобы улучшить внешний вид частотных графиков, целесообразно изменить принятые по умолчанию пропорции графического окна, растянув его в ширину.
3. Оценка ширины спектра сигнала. Используя обратное ДПФ, опре- делите минимальное число низкочастотных гармонических составляющих сигнала, содержащих не менее 90 % его энергии. Примерная последователь- ность действий:
1. Рассчитайте энергию исходного сигнала:
2 0
( )
k
E
x k
=
∑
(3.6)
2. Создайте копию вектора результатов ДПФ под новым именем.
3. Обнулите в созданной копии результатов ДПФ те элементы, которые соответствуют гармоникам с номерами, превышающими некоторое порого- вое значение
N
max
. Исходное значение
N
max следует принять равным 0.
Замечание
При определении диапазона номеров обнуляемых элементов ДПФ следует принять во внимание следующее:
• нумерация элементов массивов в MATLAB начинается с единицы;
• первое значение в массиве результатов ДПФ соответствует постоянной со-
ставляющей сигнала (
(0)
X
ɺ
);
• для получения корректного спектра вещественного сигнала в векторе ре- зультатов ДПФ должны сохраняться пары значений, соответствующие парам одинаковых по модулю положительных и отрицательных частот;
• вторая половина вектора результатов ДПФ с учетом периодичности спек- тров дискретных сигналов соответствует отрицательным частотам. В каче- стве примера в табл. 3.1 приведены разные способы нумерации элементов 8- точечного ДПФ. Двойной рамкой и полужирным шрифтом в таблице выделен диапазон элементов вектора, который необходимо обнулить в случае N
max
= 2.

40
Таблица 3.1
Номер элемента вектора MATLAB
1 2 3 4 5
6
7 8
Номер n из теоретической формулы ДПФ
0 1 2 3 4
5
6 7
Нумерация с использованием отрицательных частот 0 1 2 3
±4 −3 −2 −1 4. Вычислите обратное ДПФ и постройте график получившегося сигна- ла функцией stem.
Справка
Вычисление обратного ДПФ производится с помощью функции ifft:
x = ifft(y);
5. Вычислите энергию получившегося сигнала по (3.6). Если она мень- ше, чем 90 % энергии исходного сигнала E
0
, увеличьте N
max на единицу и повторите шаги 3–5.
Поиск значения N
max можно реализовать вручную (задавая новое значе- ние, запуская программу и анализируя результаты) либо написать MATLAB- программу, которая будет это делать автоматически.
В результате выполнения данного пункта работы должно быть получено значение N
max
, при котором сигнал после обратного БПФ содержит не менее
90 % энергии исходного сигнала. Кроме того, должен быть построен график указанного сигнала вместе с исходным в общих координатных осях.
4. Дополнение сигнала нулями. Скопируйте исходный дискретный сигнал в переменную с новым именем и добавьте к концу этой копии сигнала нулевые отсчеты в количестве, равном длине сигнала (длина вектора должна, таким образом, увеличиться в 2 раза). Вычислите ДПФ для дополненного ну- лями сигнала, постройте (с помощью функции stem) графики модуля и фазы спектральных отсчетов в одном графическом окне друг под другом (с помо- щью функции subplot).
Замечание
Чтобы улучшить внешний вид частотных графиков, целесообразно изменить принятые по умолчанию пропорции графического окна, растянув его в ширину.
5. Измерение скорости расчетов при вычислении ДПФ непосредст-
венно по теоретической формуле. Функция fft гибко выбирает алгоритм вычисления в зависимости от длины преобразуемого сигнала, поэтому изме- рить с ее помощью скорость вычислений по прямой формуле не удастся. Для реализации такого измерения придется вычислять ДПФ как произвольное

41 линейное преобразование — путем перемножения вектора сигнала и квад- ратной матрицы преобразования.
Справка
Матрица преобразования для ДПФ рассчитывается функцией dftmtx (N — размер вычисляемой матрицы ДПФ):
D = dftmtx(N);
Для измерения скорости вычислений необходимо реализовать следую- щий код (предполагается, что вектор сигнала x — строка). Для измерения времени выполнения фрагмента кода используется пара функций tic и toc
(первая запускает таймер, вторая выводит результат), вычисления повторяют- ся в цикле много раз, чтобы время выполнения кода было достаточно велико.
Внимание!
Число повторений цикла (в приведенном примере — 1000) необходимо по- добрать экспериментально, чтобы время расчета при N = 1024 составляло
примерно 1 секунду. Это время ощутимо варьируется от запуска к запуску, так что стремиться к большой точности здесь не следует.
N = 1024; % размер ДПФ
% дополнение сигнала нулями до длины N
x1 = [x zeros(1, N-length(x))];
D = dftmtx(N); % матрица ДПФ
y = zeros(1, N); % массив для результатов ДПФ
tic % старт таймера
for k = 1:1000 % цикл для измерения времени
y = x1 * D; % вычисление ДПФ по прямой формуле
end
toc % отображение измеренного времени
Запишите подобранное число повторений цикла K и занесите измерен- ное время выполнения кода при N = 1024 в табл. 3.2.
Таблица 3.2
N
Параметр
64 128 256 512 1024 2048 4096 8192
Время, с
Время одно- кратного вычис- ления ДПФ, мкс
Повторите эксперимент для остальных значений N, указанных в табли- це, не меняя подобранное число повторений цикла K.