МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ
ВЫСШЕГО ОБРАЗОВАНИЯ

«НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ»
Кафедра полупроводниковых приборов и микроэлектроники

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

Новосибирск
2013-2018

Задание №1
Тема: Спектральный анализ
Постановка задачи
1. Задать амплитуды Ak, частоты fk (Гц) и начальные фазы jk (–π…+π) пяти-семи
(m) гармоник полигармонического сигнала.
2. Задать частоту дискретизации fd (Гц). Частота дискретизации должна быть более чем в два раза выше частоты верхней гармоники.
3. Задать число отсчетов N и длительность сигнала T, учитывая соотношение fd =
N/T.
4. Сгенерировать полигармонический сигнал

где

= ( )=

∙ cos(

+

ti = Dt×i – время i-го отсчета, i = 0…N–1;
Dt = T/N
– шаг во времени между отсчетами, сек.

)

5. Построить график временной реализации сигнала.

Получить комплексный спектр с помощью преобразования Фурье, используя
функцию FFT (см. Приложение)
7. Выделить из комплексного спектра амплитуды и фазы гармоник. Определить
частоты гармоник в Гц.
8. Построить амплитудно-частотную и фазо-частотную характеристики.
9. Сравнить полученные характеристики с исходными значениями амплитуд,
частот и начальных фаз гармоник полигармонического сигнала.
10. Получить комплексный спектр, используя функцию СFFT. Выделить из комплексного спектра амплитуды и фазы гармоник. Совместить полученные
спектры на одном графике.
11. Используя формулу полигармонического сигнала, сгенерировать сигнал на
основании амплитуд и фаз, найденных в пп.7 и 10. Наложить результат на
график исходного сигнала. Добиться совпадения сигналов.
6.

2

Задание №2
Тема: Исследование свойств преобразования Фурье
Постановка задачи
Примечание: обозначения в изложении задания не совпадают с формой
записи обращения к массивам и их компонентам!
Используя наработки первого задания, выполнить исследование свойств ПФ
(ДПФ), для чего
а) задать исходный сигнал s(t) (или x(t)). Получить его спектр S(ω) (или
X(ω)) и построить АЧХ;
б) выполнить преобразование сигнала s(t)→f(t) (или x(i)→y(i)) согласно
формуле по свойству и найти спектр результата F(ω) (или Y(k)), построить график
спектра;
б) сопоставить полученный результат с формулой преобразования спектра.
При необходимости выполнить преобразование спектра S(ω)→F1(ω) (или
X(k)→Y1(k)) согласно формуле по свойству и наложить на одном графике спектры
F(ω) и F1(ω) (или Y(k) и Y1(k)), убедиться в идентичности.
Исследуемые свойства
1. Свойство А4. Изменение масштаба времени
( )= ( ∙ )

↔

( )=

1

ω

Чтобы удобно было сравнивать спектры, у них должно быть одинаковым
число гармоник и значение частоты дискретизации. Следовательно, число отсчетов сигналов тоже должно быть одно. Такие сигналы удобно получать из непрерывных функций подстановкой дискретного времени в качестве аргумента.
Для однозначности интерпретации результата ненулевые гармоники должны быть меньше частоты Найквиста.
2. Свойство А10. Спектр модулированного сигнала
( ) = ( )cos(Ω + φ) ↔

( )=

(ω − Ω) +
(ω + Ω)
2
2
Чтобы получить понятный результат, нужно чтобы спектр
сигнала s(t) содержал ненулевые гармоники сосредоточенные
близко к постоянной составляющей. Частота модуляции Ω должна
быть раза в 1,5-2 выше частоты верхней ненулевой гармоники сигнала s(t) и в 1,5-2 раза меньше частоты Найквиста.

3

3. Свойство Д1. Дополнение нулями.
Для этого свойства нужно сравнить спектры сигнала без и с дополнением
нулями.
,
0≤ <
=
,
= 0… −1
0,
≤ <
1
=
,
= 0… − 1
Частота дискретизации сигналов не меняется! Количество отсчетов сигналов разное, количество гармоник в спектрах тоже разное, но шаг между гармониками не одинаковый. Чтобы наложить спектры, нужно в качестве аргумента на
графике указать номер гармоники, умноженный на частоту первой гармоники.
4. Свойство Д4. Круговая свертка двух последовательностей
=

1

↔

, при ( −

=

=

)<0

Чтобы увидеть понятный результат, сигналы должны быть сформированными из гармоник с одинаковыми частотами, но разными амплитудами, в том
числе и нулевыми.
5. Свойство Д5. Свойство частотного сдвига
=
,
= 0… − 1
=
, при ( − ) < 0
=

∙

,

4

= 0…

−1

Дополнительная информация
К умножению сигнала на комплексную экспоненту ejΩt.
Сдвиг спектра умножением на комплексную экспоненту сдвигает все части спектра (и положительные и отрицательные частоты) в одну сторону. При этом ненулевые
гармоники из отрицательной части спектра могут переходить в область положительных
частот.
Если требуется выполнить сдвиг спектра действительного сигнала только в рамках положительных частот, сигнал из действительного (с симметричным спектром) преобразуется в комплексный с несимметричным спектром с нулевой частью отрицательных частот и удвоенными амплитудами положительных частот. Мнимая добавка к сигналу есть результат преобразования Гильберта от действительного сигнала. После выполнения сдвига (во временной области) несимметричного спектра, сигнал из комплексного преобразуется в действительный путем отбрасывания мнимой части.
Преобразование Гильберта [4, тема 17]:
̇( ) = ( ) + ∙

( )

Ортогональное дополнение сигнала (http://www.dsplib.ru/content/hilbert/hilbert.html):
( )
( − )

( )=
Обратное преобразование:

=−

Выделение действительного сигнала из комплектного:
̇( )
( )=

Функция Hilbert() выполняет преобразование Гильберта в MathCAD

5

Задание №3
Тема: Фильтрация с помощью КИХ-фильтра
Постановка задачи
1. Задать АЧХ полосового фильтра в виде массива (M = 8÷128 или более) амплитуд равноотстоящих гармоник в диапазоне от 0 до частоты дискретизации
(или частоты Найквиста). Построить график.

2. Задать ФЧХ фильтра. Допускается испробовать варианты с равными (нулевыми) фазами, фазами, распределенными по закону (например, линейному), случайные фазы, заданные с помощью генератора случайных чисел.
3. Сформировать комплексный спектр, выполнить обратное преобразование
Фурье и найти импульсную функцию фильтра. Построить график.

4. Сгенерировать полигармонический сигнал продолжительностью в 10÷20 раз
больше продолжительности импульсной функции. Частота дискретизации
сигнала и импульсной функции должна быть одинаковой (одинаковый шаг во
времени между отсчетами сигналов). Ненулевые гармоники следует расположить как внутри полосы пропускания, так и внутри полосы подавления, сосредоточенные преимущественно в районе границ перехода между полосами
пропускания и подавления. Амплитуды ненулевых гармоник лучше сделать
одинаковыми, чтобы проще было оценить результат. Построить совмещенный
график спектров фильтра и сигнала.

6

Построить график временной реализации сигнала.

5. Выполнить фильтрацию с помощью операции циклической свертки сигнала и
импульсной функции.
6. Построить спектр результата фильтрации, наложить на график АЧХ фильтра.

(первый график — с нулевой ФЧХ фильтра, второй — со случайными фазами)

7

Задание №4
Тема: Преобразование фильтров из аналогового в цифровой
Постановка задачи
1) Задать передаточную функцию аналогового фильтра в виде
( )
+⋯+
+
( )=
=
( )
+
+ ⋯+
+
Например:
1,6 + + 1
( )=
+2 +2 +1

2) Получить обобщенную частотную характеристику, построить АЧХ фильтра.
3) Выполнить преобразование фильтра из аналогового в цифровой БИХ-фильтр

Вариант 1. Записать дифференциальное уравнение, перейти от производных
к конечным разностям, получить разностное уравнение, записать передаточную
функцию W(z)
Вариант 2. Используя билинейное преобразование
2+ ∙
2 −1
=
= ∙
2− ∙
+1
записать передаточную функцию W(z)
4) Получить АЧХ цифрового фильтра
5) Совместить АЧХ аналогового и цифрового фильтров

8

Задание №5
Тема: Фильтрация с помощью БИХ-фильтра
Цель: сравнить работу КИХ фильтра из задания 3 и одного из БИХ фильтров:
Баттерворта, Чебышева первого рода, Чебышева второго рода, Кауэра (эллиптический фильтр)
Постановка задачи
1.
Задать частоту среза фильтра низких частот w0 равную верхней частоте среза фильтра из задания 3.

2.

Задать параметры коридоров АЧХ фильтра Gp, Gs и w1.

Вычислить значения параметров
= −20 ∙ log
,
=

3.

10

⁄

=

= −20 ∙ log( ) ,

− 1,

=

,

=

Рассчитать порядок фильтра

Баттерворта
Чебышева первого и второго
рода
Кауэра (эллиптический фильтр)

10

⁄

log ⁄
log( ⁄ )
arch ⁄
=
arch( ⁄ )
=

′( ) ∙ ( )
=
( ) ∙ ′( )
9

−1

[дБ]

arch() – арккосинус гиперболический

K(k) – Полный эллиптический интеграл
K’(k) – Комплиментарный
эллиптический интеграл

Полный эллиптический интеграл может быть вычислен с применением преобразования
Ландена:

( )=

2

(1 +

=⎛

) , где

⎞ ,

=

1+ 1−
⎝
⎠
Количество сомножителей M определяется выполнением условия (xi–1 – xi) < dx,
где dx — определенная заданная величина, например 10–9.
Комплиментарный эллиптический интеграл
′( ) =
1−
Порядок фильтра N округляется в большую сторону до ближайшего целого числа.
4.

5.

Определить величины L и r (r=0 или 1) в выражении
=2 +

Для эллиптического фильтра пересчитать параметр k:
=

1−

, где

=

∙

sn

∙

,

;

=

1−

Определить функцию расчета квадрата модуля АЧХ фильтра-прототипа.
1
| ( )| =
1+ ∙ ( )
Аппроксимирующая функция F определяется в зависимости от способа аппроксимации
6.

Баттерворта
Чебышева первого рода
Чебышева второго рода
Кауэра (эллиптический
фильтр)
7.

( )=

( )=

( )

1

( )=

( )=

1

( )

( ) = cos

∙ arccos( )

многочлен Чебышева N-го порядка.

RN(ω) – эллиптическая дробнорациональная функция

Построить график квадрата модуля АЧХ фильтра и допустимые границы.

10

8.

Записать передаточную функцию
8.1. для фильтра Баттерворта
( )=

∙( + ) ∙∏
1
=
;

1
( + 2 ∙ ∙ sin(
2 −1
= ∙
2

8.2.для фильтра Чебышева первого рода
∙ (− ) ∙ ∏
( )=
( − ) ∙∏ ( −2∙

+
∙ +

+

)

1
cos
−2∙ ∙ +

+

)

= − sin( ) ∙ sh( ) ;

= cos( ) ∙ ch( ) ;

=

2

8.3.для фильтра Чебышева второго рода
( )=

=

∏

( −

) ∙∏

)∙ +

= − sh( ) ;
1
1
(2 − 1);
= arsh

+

(

)

− sin( ) ∙ sh( )
1
;
=−
;
cos
∙ ch
+ sin
∙ sh
sh( )
cos( ) ∙ ch( )
1
1
(2 − 1);
;
=
= arsh
∙ ch + sin
∙ sh
2

=
cos

8.4.для фильтра Кауэра (эллиптический фильтр)
( )=

(0) =

= ∙ cd (

(0) ∙ ( −

∏
+
( − ) ∙∏ (

− ∙ ν ) ∙ ( ),

ν =−

1

∙ sn

∙

1

;

∙∏
) ∙∏ (
+

)

+
−2∙

∙ +

= Re( );

;

= ∙ cn( ∙ ν ∙ ( ), ) ;

∙ ( ),

;

=

+

)

= Im( );
=

1

1
∙ cd( ∙ ( ), )

(2 − 1)

Эллиптические функции cd(u·K(k),k), sn(u·K(k),k) могут быть рассчитаны через
рекуррентные соотношения (М получается из заданной точности х)
=⎛
⎝

1+ 1−

⎞ ;
⎠

11

= ;

= 1,2, … ,

=

При начальном значении

1

1+
+

∙

,

=

,

− 1, … , 2, 1;

= cos( ∙ /2) ;

cd( ∙ ( ), ) =

= sin( ∙ /2) ;

sn( ∙ ( ), ) =

получается значение w0, которое и есть искомое значение функции cd(u·K(k),k).
При начальном значении

получается значение w0, которое есть искомое значение функции sn(u·K(k),k).
Обратные эллиптические функции cd–1(ω,k), sn–1(ω,k) могут быть рассчитаны через рекуррентные соотношения
2∙
=
;
= ; = 1,2,3, … ,
(1 + ) ∙ 1 + 1 +
∙
cd ( , ) =

9.

2

sn ( , ) =

2

arcos(
arsin(

)

)

Построить график обобщенной частотной характеристики W(jω).

10. Определить коэффициенты рекурсивного цифрового фильтра в соответствии с алгоритмом по заданию 4.
11. Реализовать программно алгоритм фильтрации рекурсивным цифровым
фильтром.
12. Выполнить фильтрацию сигнала из задания 3. Получить спектр отфильтрованного сигнала, построить график АЧХ. Сравнить с АЧХ нерекурсивного фильтра.

12

Литература
1. Сергиенко А.Б. Цифровая обработка сигналов учебное пособие для вузов. —
СПб. Питер, 2007. (Расстановочный шифр 621.39 С323, Bib Id 65599)
2. Стивен Смит Цифровая обработка сигналов практическое руководство для
инженеров и научных работников; пер. с англ. Ю. А. Линовича, С. В. Витязева, И. С. Гусинского. — М. Додэка-XXI, 2011. (Расстановочный шифр 621.37
C509, Bib Id 162582)
3. Яковлев А.Н. Основы теории сигналов в примерах, упражнениях и заданиях
[учебное пособие для радиотехнических направлений и специальностей]. —
Новосибирск Изд-во НГТУ, 2012. (Расстановочный шифр 621.37 Я474, Bib Id
174664)
4. Давыдов А.В. Сигналы и линейные системы: Тематические лекции. / Екатеринбург: УГГУ, ИГиГ, кафедра геоинформатики. – 2008. /
http://prodav.exponenta.ru/signals/index.html
5. Давыдов А.В. Цифровая обработка сигналов: Тематические лекции. / Екатеринбург: УГГУ, ИГиГ, кафедра геоинформатики. – 2007. /
http://prodav.exponenta.ru/dsp/index.html
6. Глинченко, А. С. Цифровая обработка сигналов. Версия 1.0 [Электронный
ресурс] : курс лекций / А. С. Глинченко. – Электрон. дан. (3 Мб). – Красноярск : ИПК СФУ, 2008. – url: . http://files.lib.sfukras.ru/ebibl/umkd/50/u_lectures.pdf
7. Васюков В.Н. Цифровая обработка сигналов и сигнальные процессы в системах подвижной радиосвязи [учебник]. — Новосибирск Изд-во НГТУ, 2006.
(Расстановочный шифр 621.37 В201, Bib Id 62388)
8. Баскаков С.И. Радиотехнические цепи и сигналы: Учеб для вузов по спец «Радиотехника». — М.: Высш. школа, 2000.
9. Гоноровский И.С. Радиотехнические цепи и сигналы. — М.:Сов.радио, 1967
10.Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения: Пер. с
англ. – М.: Мир, 1990. – 584 с.

13

Приложения
Справочник по функциям MathCAD
Преобразование Фурье
Прямое преобразование
Функция
Формула
c := FFT(x)
c := fft(x)

сk =
сk =

1
N

åxe

Обратное преобразование
Функция
Формула

- j 2pki / N

xi =

x := IFFT(c)

i

åc e
k

i

1
N

j 2pki / N

k

åxe

j 2pki / N

x := ifft(c)

i

xi =

i

1
N

åc e

- j 2pki / N

k

k

x – действительный массив временных отсчетов размерностью N;
N – число временных отсчетов, N = 2m ;
i – номер временного отсчета от 0 до N–1;
c – массив комплексных коэффициентов спектра размерностью M+1;
M=N/2 – число гармоник половины спектра (до частоты Найквиста);
k – номер гармоники от 0 до M;
ak = Ak cosjk – синфазная компонента;
bk = Ak sinjk – квадратурная компонента;
Ak = ck = ak2 + bk2 – амплитуда k-той гармоники;
jk = atan( bk / ak )
ak = Re(ck)
bk = ±Im(ck)

– фаза k-той гармоники
– связь синфазной компоненты с комплексным коэффициентом;
– связь квадратурной компоненты с комплексным коэффициентом
“+” в случае Ak cos(kω1t+jk), “–” в случае Ak cos(kω1t–jk);
T
– длительность сигнала, сек;
Dt = T/N – шаг во времени между отсчетами, сек;
ti = Dt×i – время i-го отсчета;
f1 = 1/T – частота первой гармоники;
fk = f1×k – частота k-той гармоники, Гц;
fd = N/T – частота дискретизации, Гц;
ω = 2π×f – круговая частота, рад/с;
ωk×ti = ω1×k × Dt×i = 2π/T×k × T/N×i = 2π/N × k × i безразмерная величина.
Прямое преобразование
Функция
Формула
c := CFFT(x)
c := cfft(x)

сk =
сk =

1
N

åxe

1
N

Обратное преобразование
Функция
Формула

- j 2pki / N

x := ICFFT(c)

i

xi =

åc e
k

i

åxe

j 2pki / N

k

j 2pki / N

x := icfft(c)

i

i

xi =

1
N

åc e
k

k

x – комплексный массив временных отсчетов размерностью N;
N – число временных отсчетов, N – произвольное целое число ;
i – номер временного отсчета от 0 до N–1;
c – массив комплексных коэффициентов спектра размерностью M;
M=N – число гармоник спектра от нулевой до частоты дискретизации (не включая);
k – номер гармоники от 0 до M–1;
14

- j 2pki / N

Функции работы со звуковыми файлами в MathCAD
Функция GETWAVINFO(filename) возвращает вектор из четырех чисел:
- number of channels — число каналов (1-моно или 2-стерео)
- the sample rate — частота дискретизации, Гц (11025, 22050, 44100)
- the bit resolution — число бит на отсчет (8, 16)
- the average bytes per second for a WAV file — число байт в секунду
Функция READWAV(filename) возвращает вектор отсчетов для монофонического
сигнала или матрицу из двух столбцов для стереосигнала.
Функция WRITEWAV(filename, rate, res) := x выводит массив x в файл с именем
из переменной filename, rate — частота дискретизации (Гц), res — число бит на отсчет
(8, 16).
Примечание
Для разрешения 8 бит на отсчет значения отсчетов должны быть в диапазоне от 0
до 255.
Для разрешения 16 бит на отсчет значения отсчетов должны быть в диапазоне от
–32768 до +32767.
Пример 1
fd:=4096
dt:=1/fd
T:=1 N:=T/dt
xi := 127 + 127 · sin( 2π · 400 · dt · i )

i:=0…N

WRITEWAV(“440.wav”, fd, 8) := x
Пример 2
fd:=4096
dt:=1/fd
T:=1 N:=T/dt
xi,0 := 32767 · sin( 2π · 400 · dt · i )
xi,1 := 32767 · cos( 2π · 400 · dt · i )
WRITEWAV(“Stereo.wav”, fd, 16) := x
Пример 3
2
4096
GETWAVINFO(“Stereo.wav”) =
16
16384
x := READWAV(“Stereo.wav”)
N := rows(x) = 4097

15

i:=0…N