Цели работы:
1)
практическое ознакомление с реализацией процедуры измерения спектра случайного стационарного процесса (ССП) в среде MATLAB;
2)
моделирование нескольких примеров применения спектрально- го анализа случайных процессов, а именно:
– выявление периодического сигнала, который маскируется шумом;
– измерение частот аддитивной смеси нескольких гармонических колебаний;
– измерение частоты основного тона голосового сигнала.
165
1. Основные теоретические сведения
подробно описаны в гл. 2.
2. Задания и методические указания по выполнению работы
1.
Смоделируйте задачу выявления периодического сигнала, кото- рый маскируется шумом, при следующих условиях:
– ССП ( )
Y t представляет собой аддитивную смесь гармониче- ского процесса
0
( )
cos(2
)
S t
A
f t
π
ϕ
=
+
с неизвестной амплитудой
A
, частотой
0
f (значение которой находится в пределах 80 – 220 Гц), случай- ной фазой, равномерно распределенной на интервале
[0,2 ]
π
, и гауссовско- го белого шума в полосе частот 0 – 5 кГц;
– отношение сигнал-шум этой смеси равняется вх
ρ
2.
Смоделируйте в среде MATLAB аддитивную смесь с заданными параметрами.
3.
Рассчитайте объем
N экспериментальной выборки отсчетов сге- нерированной смеси, которая необходима для обеспечения отношения сигнал-шум вых
10 дБ
ρ
=
на выходе цифрового спектроанализатора, вы- числяющего периодограмму.
4.
Вычислите и постройте график разных оценок спектра ССП
( )
Y t :
– периодограмма;
– оценка Бартлетта;
– оценка Велча с 50 %-м перекрытием сегментов.
5.
Вычислите и постройте график тех же оценок спектра ССП ( )
Y t при условиях увеличения N в 100 раз, при длине сегмента
1
N , равной
«старому» значению N
6.
Осуществите спектральный анализ суммы двух гармонических колебаний (без фонового шума) при следующих условиях:
– частота первого колебания
01 0
f
f
= , частота второго колебания
02 0
5
f
f
f
=
+ Δ , где
1 2
f
T
B N
Δ =
=
;
– мощности обеих гармоник одинаковы;
–
N =1024.
7.
Осуществите натурный эксперимент по измерению частоты основного тона голосового сигнала с применением спектрального анализа.
Используя предыдущие записанные слова, введите их в рабочее простран- ство программы sptool. Осуществите спектральный анализ со следующим измерением частоты основного тона:
– периодограмма;
– оценка Бартлетта;
– оценка Велча с 50 %-м перекрытием сегментов.
166
3. Вспомогательная теория
1.
Для выполнения данной лабораторной работы введите в рабочее пространство программы MATLAB числовые значения параметров со- гласно табл. 2.
Таблица 2
Варианты значений числовых параметров
Вариант
Параметры
1 2 3 4 5 6 7 8 вх
ρ
-–10 L-–11 L-–12 L-–13 L-–14 L-–15 L-–16 –17 0
f
80 100 120 140 160 180 200 220
В
5000 5000 5000 5000 5000 5000 5000 5000
Слово
Примечание. В графу «Слово» каждый записывает свою фамилию.
2.
Сформируйте аддитивную смесь гармонического процесса и бе- лого шума ( )
( )
( )
Y t
S t
t
ξ
=
+
Для сигнала
0
( )
cos(2
)
S t
A
f t
π
ϕ
=
+
и шума
( )
t
ξ
входное отношение сигнал-шум
2 2
A
вх
D
ρ
ξ
=
, или в децибелах
2 2
10 lg
вх дБ
A
D
ρ
ξ
= ⋅
Если примем
ξ
1,
D
= тогда
10 lg 2 20 10
вх дБ
A
ρ
⋅
+
=
Если за сигнал на выходе спектроанализатора принять высоту пика математического ожидания отрезка гармонического процесса, а за шум – математическое ожидание уровня спектра шума, то
вых сп
вх
TB
ρ
ρ
=
⋅
В действительности за уровень шума принимают не уровень спектра шума, а среднеквадратичную погрешность оценки спектра шума. Можно показать, что для оценки в виде периодограммы эти два разных определе- ния отношения сигнал-шум совпадают. Выигрыш в отношении сигнал- шум за счет спектрального анализа достигает величины
,
вых сп
вх
TB
ρ
ρ
=
где
T
– длина отрезка,
B
– верхняя предельная частота шума.
167
Отсюда
(
)
0,1 10 log log log
10
log log
вых сп
вх
вых сп
вх
p
p
p
p
TB
TB
−
−
⋅
=
⇒
=
3.
Рассчитайте нужное количество отсчетов процесса по следующей формуле:
0,1(
),
log log
2 2 10
вых сп
вх
p
p
T
N
TB
t
−
=
=
= ⋅
Δ
где log
вых сп
p
и log
вх
p
– логарифмические значения входной и выход- ной спектральной мощности соответственно.
4.
Вычислите периодограмму,которую в среде MATLAB можно рассчитать с помощью функции periodogram.
5.
Найтите оценки. В данной лабораторной работе для оценки Велча следует использовать свой вариант окна по табл. 3.
6.
При вычислении и построении графиков разных оценок спектра
ССП
( )
Y t (п. 4) используйте функцию subplot.
7.
Измерьте частоты основного тона голосового сигнала с примене- нием спектрального метода. Активизируйте программу sptool, экспорти- руйте в ее рабочее пространство сигнал, который отвечает целому слову и отдельному гласному звуку из него.
Таблица 3
Типы окон
Номер варианта
Окно
Синтаксис
1
Прямоугольное
Window = boxcar (N
1
)
2
Треугольное
Window = triang (N
1
)
3
Бартлетта
Window = bartlett (N
1
)
4
Блэкмена
Window = blackman (N
1
, ‘symmetric’)
5
Чебышева (40db)
Window = chebyshev (N
1
, 40)
6
Хэмминга
Window = hamming (N
1
, ‘symmetric’)
7
Хэннинга
Window = hanning (N
1
, ‘symmetric’)
8
Кайзера (beta = 6)
Window = kaiser (N
1
, 6)
Вычислите (в среде программы sptool) те же разновидности оценок спектра: периодограмму, оценку Бартлетта, оценку Велча. Параметр Nfft для оценки в виде периодограммы вычислите по формуле
(
)
2 1
d
N f f t
t
t F
=
−
, где
1
t и
2
t – моменты начала и конца фрагмента (
1
t = 0);
d
F
=11025 Гц.
168
ЛАБОРАТОРНАЯ РАБОТА № 5. РАСЧЕТ ЦИФРОВЫХ ФИЛЬТРОВ
В СРЕДЕ MATLAB
Цель работы:
провести исследования по применению фильтров для обработки речевых сигналов в системе MATLAB.
1. Основные теоретические сведения
1.1. Расчет коэффициентов цифрового фильтра
Для расчета коэффициентов
(
)
при
0, ,
k
a
k
N
…
=
цифрового фильтра
(ЦФ) используют следующие расчетные формулы: c
c c
c c
s in
(
)
a a
a
π
π
k
k
k
S
k
k
Ω
Ω
Ω
=
=
=
Ω
−
Ω
, где c
c
2
c
d
d
f
f
f
ω
π
=
=
Ω
,
( )
sin
x
x x
S a
=
, c
ω
– частота среза (иногда ее обозначают
â
ω
и называют «верхняя граничная частота»),
d
f
– частота дискретизации и сf– частота среза.
Рассмотрим пример ФНЧ 6-го порядка (
3
N
= ) для сf=25 Гц,
fd =100 Гц, коэффициенты фильтра:
0
a
0,5
=
,
1 1
a a
1 π 0,3183
−
=
=
=
,
2 2
a a
0
−
=
= ,
3 3
a a
1 3π
0,1061
−
=
= −
= −
Интересно и полезно проверить, действительно ли эти коэффициен- ты с точностью до коэффициента
1
t
f d
Δ =
совпадают со значениями
ИПХ
(
)
h k t
Δ :
(
)
a
k
th k t
= Δ
Δ .
Обратное преобразование Фурье от функции
( )
rect
2
c
f
H f
f
⎛
⎞
=
⎜
⎟
⎜
⎟
⎝
⎠
дает выражение для ИПХ аналогового фильтра:
( )
(
)
c c
2 2
a
h t
S
t
f
f
π
=
, откуда следует
(
)
(
)
(
)
c c
c
1 2
2
a
a
k
c
d
h k t
S
k t
S
k
a
f
f
f
t
π
π
Ω
Δ
=
Δ
=
=
Ω
Δ
Сравнение позволяет обнаружить различие «по вертикали» ИПХ аналогового и цифрового фильтров.
Рассмотрим различие «по горизонтали».
169
Расстояние между нулями функции
h(t) c
1 2
fτ
Δ =
. В рассмотренном примере было принято
4 1 4
cdtff=
=
Δ , откуда следует
2
tτ
Δ = Δ . Дей- ствительно, коэффициент
2
a
0
= . Ясно также, что все четные коэффициен- ты должны быть равны нулю.
Для иллюстрирования примера, построим график функции
h(t) (рис. .4).
% График ИПХ непрерывного ФНЧ dt=1; % шаг дискретизации; dtau=2*dt; %
расстояние между нулями; fc=1/(2*dtau); % частота среза t=-3*dtau:0.1:3*dtau; % время h=2*fc*sinc(2*fc*t); % ИПХ plot(t,h) % график grid on % сетка на графике
Рис. 8.4. График ИПХ непрерывного ФНЧ
Построим частотную характеристику рассчитанного фильтра
(рис. 8.5):
( )
0 1
2
c o s
kNkaak tH dω
ω
=
=
+
Δ
∑
a0=0.5; ak=[0.3183 0 -0.1061]; % коэффициенты фильтра c 1 по 3 dt=1; N=3; df=0.02;
170
f=-0.5:df:1.5; % диапазон частот sum=0; for k=1:N, sum=sum+ak(k)*cos(2*pi*f*k*dt); end;
H=a0+2*sum; % ЧХ; plot(f,H); grid on
Рис. 8.5. Частотная характеристика ФНЧ 3-го порядка (N = 3)
Как видно из графика на рис. 8.5, частотная характеристика фильтра
3-го порядка существенно отличается от прямоугольника с относительной частотой среза c
1 4
d
f
f
=
. Полагая, что причина тому – малый порядок фильтра, рассмотрим случай
9
N
= и досчитаем недостающие коэффици- енты ИПХ с нечетными номерами (коэффициенты с четными номерами равны нулю) (рис. 8.6):
5
a
0,0637
=
;
7
a
0,0455
= −
;
9
a
0,0354
=
% Частотная характеристика ФНЧ a0=0.5; ak=[0.3183 0 -0.1061 0 0.0637 0 -0.0455 0 0.0354]; % коэффициенты фильтра c 1 по 9 dt=1; N=9; df=0.02; f=-0.5:df:1.5; for k=1:N, sum=sum+ak(k)*cos(2*pi*f*k*dt); end
H=a0+2*sum; plot(f,H); grid on
171
Как следует из рис. 8.6, частотная характеристика действительно ста- новится более прямоугольной, но степень осцилляций слева и справа час- тоты среза не уменьшается.
Применим треугольное окно к ИПХ фильтра 18-го порядка (
9
N
= ) (рис. 8.7).
% Частотная характеристика ФНЧ a0=0.5; ak=[0.3183 0 -0.1061 0 0.0637 0 -0.0455 0 0.0354]; fd=100; dt=1/fd;
N=9; r=1:N; k(r)=ak(r).*(1-r/N); df=0.02*fd; f=-0.5:df:1.5; sum=0; for k=1:N, sum=sum+ak(k)*cos(2*pi*f*k*dt); end
H=a0+2*sum; plot(f,H); grid on
Рис. 8.6. Частотная характеристика ФНЧ 3-го порядка (N=9)
Рис. 8.7. Частотная характеристика ФНЧ 18-го порядка (N=9)
172
Как и следовало ожидать, осцилляции уменьшились, однако за это пришлось «заплатить» уменьшением крутизны склонов частотной харак- теристики в районе частоты среза.
1.2. Расчет цифровых фильтров в среде MATLAB
В среде MATLAB цифровые фильтры можно рассчитывать, по меньшей мере, тремя способами: из командного окна; с помощью пакета
sptool; с помощью пакета fdatool.
1.
Расчет цифровых фильтров из командного окна
Функция fir1 реализует вычисления по методу обратного преобразо- вания Фурье с использованием окон: a=firl(n,Wn,’ftype’,window,’normalization’), где n – порядок фильтра (количество коэффициентов равно n+1), его лучше задавать четным (так как для некоторых типов фильтров нечетное n хотя и можно задавать, но результат будет такой, как если бы задавался порядок на единицу больше);
Wn – относительная частота среза (по отношению к частоте Найкви- ста, которая принимается равной единице); представляет собой вектор из двух чисел, если фильтр полосовой или режекторный и вектор из m пар чисел, если фильтр многополосный, из m полос;
‘
ftype’ – тип фильтра (если отсутствует – ФНЧ;
‘
high’ – ФВЧ;
‘
stop’ – режекторный;
‘
DC-1’ – многополосный пропускающий;
‘
DC-0’ – многопо- лосный режекторный); window – вектор-столбец из n+1 элементов (по умолчанию применя- ется окно Хэмминга hamming(n+1));
‘normalization’ – нормировка ИПХ (по умолчанию значение ‘scale’ – единичное значение коэффициента передачи в центре полосы пропуска- ния; ‘noscale’ – нормировка не производится).
Пример
% расчет коэффициентов КИХ-фильтра с нормализацией window=rectwin(7) % синтез прямоугольного окна из 7 отсчетов a=fir1(6,0.5,window);
Результат
a =[–0,1148 0,0000 0,3443 0,5409 0,3443 0,0000 –0,1148]
Сравнивая эти результаты с рассчитанными вручную коэффициен- тами, нетрудно увидеть разницу. Например, расчеты вручную дают
0
а
0,5
=
, тогда как в MATLAB мы получили
0
а
0,5409
=
. Естественно предположить, что причина тому – проводимая по умолчанию нормировка
ИПХ. Проверяем это предположение, задавая в программе значение
’noscale’ для параметра нормализации: window=rectwin(7); a=fir1(6,0.5,window,’noscale’).
173
Результат: a = [– 0,1061 0,0000 0,3183 0,5000 0,3183 0,0000 – 0,1061].
2.
Расчет цифровых фильтров с помощью пакета sptool
Для активизации пакета нужно в командном окне набрать команду
sptool или открыть его следующим образом: start
→
toolboxes
→
More
…
→
Signal
Processing
→
Signal
Processing tool
(sptool). Затем в поя- вившемся окне в колонке кнопок Filters необходимо нажать кнопку New.
В возникшем окне Filter Designer нужно:
- задать частоту дискретизации 100 Гц;
- выбрать в позиции Design Method выбрать FIR-значение Kaiser
Window FIR (выбираем из трех вариантов: Equiripple FIR, Least Square FIR и Kaiser Window FIR);
- отключить флажок Minimum Order; задать Order = 6; Type =
= lowpass; Passband Fc = 25;
- отключить флажок Autodesign; закрыть окно Filter Designer;
- в окне sptool в колонке Filters нажать кнопку View.
В появившемся окне Filter Viewer наблюдаем графики АЧХ, ФЧХ,
ИПХ (ИПХ наблюдаем после активизации соответствующего флажка).
Мы выбрали в позиции Design Method значение Kaiser Window FIR.
Кроме данного алгоритма есть еще два: Equiripple FIR, Least Square FIR.
Из всех трех алгоритмов только алгоритм Кайзера реализует метод обрат- ного преобразования Фурье с весовым окном Кайзера. При значении пара- метра «0» окно Кайзера превращается в обычное прямоугольное окно.
Как показывает эксперимент, рассчитанные таким образом коэффи- циенты ФНЧ оказываются ненормированными, т.е. в точности равными вычисленным вручную.
Чтобы узнать значения коэффициентов, нужно:
- активизировать график ИПХ, щелкнув по нему мышкой;
- активизировать вертикальные маркеры (кнопкой, расположенной под меню);
- поместить один из маркеров (всего имеется два маркера: первый изображается сплошной вертикальной линией, второй пунктирной) напро- тив нужного отсчета ИПХ;
- рассчитать значение отсчета ИПХ в специальном окне.
Другой способ в зависимости от версии MATLAB – нажать на спе- циальную кнопку Filter coefficients.
174 3.
Расчет фильтра
{ }
k
a с помощью пакета fdatool
Для активизации пакета нужно в командном окне набрать команду
fdatool, затем в появившемся окне в разделе Design Filter задать:
- Designe Method: FIR=Window;
- Window Specifications: Window = Rectangular;
- Filter order: Specify order = 6;
- частоту дискретизации (Fs = 100 Гц);
- Filter Type = lowpass;
- Passband Fc=25.
С помощью кнопок под меню включаем режим просмотра коэффи- циентов фильтра.
В результате проведенных расчетов убеждаемся, что здесь по умол- чанию проводится нормирование ИПХ (см. выше fir1).
2. Задания и методические указания по выполнению работы
Необходимо выполнить следующее:
1.
С помощью микрофонной гарнитуры введите в компьютер рече- вой сигнал (фамилию студента) с параметрами: длительность – несколько секунд; частота дискретизации – 11 025 Гц. Для ввода сигнала в компью- тер удобно использовать программу «Звукозапись» из раздела «Стандарт- ные – Развлечения».
2.
С помощью команды y = wavread('filename') импортируйте рече- вой сигнал в среду MATLAB, постройте график сигнала с помощью ко- манды plot.
3.
Запустите программу sptool и импортируйте сигнал y в её среду.
4.
Визуализируйте и прослушайте введенный сигнал с помощью ин- струментария программы sptool
.
5.
Синтезируйте НЧ фильтр с окном Кайзера минимального поряд- ка; граничная частота полосы пропускания 1 500 Гц; неравномерность в полосе пропускания 1 дБ; граничная частота полосы задержания 2 000 Гц; минимальное затухание в полосе задержания 40 дБ.
6.
Примените синтезированный НЧ фильтр к вашему речевому сигналу.
7.
Прослушайте сигнал, полученный в результате НЧ фильтрации, постройте его график и ответьте на следующие вопросы.
Как и почему изменилось звучание речевого сигнала после фильтрации?
Что произойдет, если граничные частоты полос пропускания и за- держания уменьшить вдвое? Экспериментально проверьте свои предполо- жения.
8.
Подберите оптимальные с вашей точки зрения параметры НЧ фильтра. Попробуйте обосновать свой выбор.
175
ЛАБОРАТОРНАЯ РАБОТА № 6. КОДИРОВАНИЕ РЕЧЕВЫХ ДАННЫХ
НА ОСНОВЕ ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ
Цель работы: получить навыки применения метода линейного предсказания и расчета коэффициентов фильтра предсказания речевого сигнала в системе MATLAB.
1. Основные теоретические сведения Основной принцип метода линейного предсказания состоит в том, что текущий отсчет речевого сигнала можно аппроксимировать линейной комбинацией предшествующих отсчетов. Коэффициенты предсказания – это весовые коэффициенты, используемые в линейной комбинации. Они определяются однозначно из условия минимизации среднего квадрата раз- ности между отсчетами речевого сигнала и их предсказанными значениями
(на конечном интервале).
Основные положения метода линейного предсказания хорошо согла- суются с моделью речеобразования, где речевой сигнал представляется в виде сигнала на выходе линейной
системы с переменными во времени па- раметрами, возбуждаемой квазипериодическими импульсами (в пределах вокализованного сегмента) или случайным шумом (на невокализованном сегменте). Метод линейного предсказания позволяет точно и надежно оце- нить параметры этой линейной системы с переменными коэффициентами.
Модель имеет следующие параметры: классификатор вокализован- ных и невокализованных звуков, период основного тона для вокализован- ных звуков, коэффициент усиления
g , коэффициенты
{ }
ka цифрового фильтра. Все эти параметры медленно меняются во времени. Остановимся подробно на задаче определения коэффициентов цифрового фильтра
{ }
ka .
Предположим, что отсчеты речевого сигнала
( )
s n связаны с сигна- лами возбуждения
( )
u n разностным уравнением
( )
(
)
( )
1
kps na s n kg u nk=
− + ⋅
∑
=
. (8.7)
В этом случае передаточная функция линейной системы с входом
( )
u n и выходом
( )
s n имеет вид
( ) ( ) ( )
1 1
pkkkH zS z U zga z=
⎛
⎞
−
=
=
−
⋅
∑
⎜
⎟
⎝
⎠
, где
z – формальная переменная;
( )
S z и
( )
U z –
Z-преобразования речево- го сигнала
( )
s n и сигнала возбуждения
( )
u n .
176
Линейный предсказатель с коэффициентами
{ }
k
a определяется как система, на выходе которой в момент времени n имеем
( )
(
)
1
k
p
k
s n
s n k
α
=
=
−
∑
. (8.8)
Системная функция предсказателя p-го порядка представляет собой полином вида
( )
1
p
k
k
k
P z
z
α
=
−
=
⋅
∑
Погрешность предсказания определяется как
( ) ( )
(
)
1
k
p
k
e n
s n
s n k
α
=
=
−
−
∑
Другими словами, погрешность предсказания представляет собой сигнал на выходе системы с передаточной функцией
( )
1 1
α
p
k
k
k
A z
z
=
−
= −
⋅
∑
Таким образом, если сигнал точно удовлетворяет модели (8.7) и
k
k
a
α
=
, то
( )
( )
e n
g u n
= ⋅
. Отсюда следует, что фильтр погрешности пред- сказания
( )
A z – обратный фильтр для системы с передаточной функцией
( )
H z , соответствующей уравнению (8.7), т.е.
( )
( )
H z
g A z
=
Основная задача анализа на основе линейного предсказания заклю- чается в определении параметров
{ }
k
a по речевому сигналу. При этом предполагается, что полученные параметры являются параметрами сис- темной функции
( )
H z в модели речеобразования. Вследствие изменения свойств речевого сигнала во времени коэффициенты предсказания должны оцениваться на коротких сегментах речи – кадрах.
В качестве критерия, по которому производится оптимизация синте- за фильтра
( )
A z , удобно взять минимум суммы квадратов погрешностей линейного предсказания на сегменте (кадре) речевого сигнала. Основные причины для выбора такого критерия следующие: получающиеся уравне- ния – линейные, они относительно просто решаются и дают хорошие ре- зультаты.
Пусть
1 0
,
n n
⎡
⎤
⎣
⎦
– некоторый интервал. Сумма квадратов погрешно- стей линейного предсказания определяется следующим образом:
( )
2 1
0
n
n n
E
n
e
=
=
∑
177
Параметры ak можно получить, минимизируя E. Подставим (8.8) в выражение для Е и приравняем к нулю производные
,
1,2, ,
k
E a
k
p
∂ ∂
=
… .
Получаем
( )
(
)
(
)
(
)
( )
( ) (
)
(
) (
)
1 0
0 0
0 1
1 1
1 1
1 1
2 1
2 2
p
n n
n
n
p
n n
i
n n
n
p
p
i j
i
j
n n
n
E
s n
s n
s n
p
a
a
n
s n
s n
i
s
ai
s n
i
s n
j
a a
=
=
=
=
= =
=
=
−
⋅
−
− … −
⋅
−
=
∑
=
−
⋅
− +
∑
∑
∑
+
− ⋅
−
∑ ∑
∑
(8.9)
Дифференцируем (8.9) по
,
1, 2, ,
k
a
k
p
…
=
:
( ) (
)
(
) (
)
1 0
1 1
0 0
k
n
E a
s n s n k
n n
n
p
a
s n k s n i
i
i
n n
∂ ∂
=
⋅
− −
∑
=
−
− ⋅
− =
∑
∑
=
=
(8.10)
Заменив k на j , получим систему p линейных уравнений относи- тельно p неизвестных
1 2
,
,...,
p
a a
a
0
,
1, 2, , ,
1
j
i j
p
a c
c
j
p
i
i
=
=
…
∑
=
(8.11) где
(
) (
)
1 0
i j
j i
n
c
c
s n i s n
j
n n
=
=
− ⋅
−
∑
=
. (8.12)
Данная система называется системой уравнений Юла – Волкера. Ре- шив её, нетрудно оценить и минимальную достижимую погрешность предсказания. Для этого подставим (8.11) в (8.9)
00 0
2 1
1 1
i
i j
i
i
j
p
p
p
E c
a c
a
a c
i
i
j
=
−
+
∑
∑
∑
=
=
=
и, используя (8.11), упростим это выражение. В результате получим
00 0
1
i
p
E c
a c i
i
=
− ∑
=
Для определения коэффициентов
k
a из уравнений Юла – Волкера необходимо знать величины
,
0,1,
, ,
1,2,
,
i j
c
i
p j
p
=
…
=
…
. Имеется два подхода к вычислению этих величин. Один называется ковариационным методом, второй – автокорреляционным.
178 1. Автокорреляционный метод
Для этого метода примем пределы анализа равными (
,
)
−∞ + ∞ , ин- тервал анализа (0 , )
N , причем сигнал обнуляется вне интервала анализа, т.е.
( )
0 при
0,
s n
n
n N
=
<
≥ . Такие пределы позволяют упростить выра- жение для
,
1,2, , ,
0,1, ,
c
i
p j
p
i j
=
=
…
…
:
(
) (
)
( )
(
)
1 1
0 0
i j
N
i
j
N
p
c
s n i s n
j
s n s n i
j
n
n
− − −
+ −
=
− ⋅
− =
⋅
+ −
∑
∑
=
=
В этом случае cij являются функциями величины i j
− и с точно- стью до множителя совпадают с оценками автокорреляционной функции
( )
ˆR
τ
сигнала
( )
s n , вычисленными при
i
j
τ
= − ,
(
)
( )
(
)
1
ˆ
1 0
i j
N
i
j
R i
j
c
N
N
s n s n i
j
n
− − −
−
=
=
⋅
+ −
∑
=
Разделив уравнения в системе (8.11) на N, получаем систему уравне- ний Юла – Волкера для автокорреляционного метода
(
)
( )
ˆ
ˆ
,
1, 2,
,
1
i
p
a R i
j
R j
j
p
i
⋅
−
=
=
…
∑
=
(8.13)
В матричном виде система может быть записана как
a R b
× =
, где
(
)
1 2
ˆ
ˆ
ˆ
( ,
,...,
),
(1),
(2),
,
( )
p
a
a a
a
b
R
R
R p
=
=
…
,
ˆ
ˆ
ˆ
(0)
(1)
(
1)
ˆ
ˆ
ˆ
(1)
(0)
(
2)
ˆ
ˆ
ˆ
(
1)
(
2)
(0)
R
R
R p
R
R
R p
R
R p
R p
R
⎡
⎤
−
⎢
⎥
−
⎢
⎥
= ⎢
⎥
⎢
⎥
⎢
⎥
−
−
⎣
⎦
…
…
…
…
…
…
…
Матрица R в автокорреляционном методе обладает двумя важными свойствами. Она симметрическая (ее элементы, симметричные относи- тельно главной диагонали, равны) и теплицева (каждая следующая строка получается из предыдущей сдвигом вправо). Структура теплицевой матри- цы позволяет решить систему (8.13) особенно просто. Для ее решения по алгоритму Левинсона – Дарбина, описание которого приводится ниже, тре- буется порядка 2
p операций. Решение произвольной системы из p урав- нений с p неизвестными потребовало бы порядка 3
p операций.
179 2. Ковариационный метод
В этом методе выбирается
0 1
0,
1
N
n
n
=
= − , а сигнал
( )
s n не ограни- чивается. При этом для величин
,
1, 2, , ,
0, 1,
,
c
i
p j
p
i j
…
…
=
=
(
) (
)
1 0
N
c
s n i s n
j
i j
n
−
=
− ⋅
−
∑
=
Изменив индекс суммирования, это выражение можно представить в виде
( ) (
)
1
N i
c
s n s n i
j
i j
n
i
− −
=
⋅
+ −
∑
= −
, при
1, 2, , ,
0, 1,
, .
i
p j
p
…
…
=
=
(8.14)
Выражение (8.14) похоже на выражение для c i j для автокорреля- ционного метода, но имеет другие пределы суммирования. В (8.14) ис- пользуются значения сигнала вне интервала 0 1
n N
≤ ≤ − . Другими слова- ми, для вычисления c i j в ковариационном методе необходимо знать зна- чения сигнала
( ) (
)
(
)
,
1 , ,
1
s
p s
p
s N
−
− +
−
…
, т.е. сигнал должен быть извес- тен на несколько большем интервале, чем в автокорреляционном методе.
Однако, как правило, p
N и данное требование не очень существенно.
Этот метод приводит не к автокорреляционной, а к взаимной корреляци- онной функции между двумя очень сходными, но не одинаковыми сегмен- тами речевого сигнала конечной длительности
( )
(
) (
)
1
ˆ ,
1 0
N
R i j
c
N
N
s n i s n
j
i j
n
−
=
=
− ⋅
−
∑
=
Нетрудно увидеть, что
( )
( )
ˆ
ˆ
,
,
R i j
R j i
=
, однако
( )
ˆ ,
R i j не является функцией от i
j
− , как это было в автокорреляционном методе. Разделив все уравнения в системе (8.10) на N, получаем систему уравнений Юла –
Волкера для ковариационного метода
( )
( )
ˆ
ˆ
,
0, ,
1, 2,
,
1
p
a R i j
R
j
j
p
i
i
⋅
=
=
…
∑
=
. (8.15)
В матричном виде система уравнений (8.15) имеет вид
a P c
× =
, где
(
)
(
)
ˆ
ˆ
ˆ
,
,
,
,
(0,1),
(0,2),
,
(0, )
1 2
a
c
R
R
R
p
a
a
a p
=
…
=
…
,
ˆ
ˆ
ˆ
(1,1)
(1,2)
(1, )
ˆ
ˆ
ˆ
(2,1)
(2,2)
(2, )
ˆ
ˆ
ˆ
( ,1)
( ,2)
( , )
R
R
R
p
R
R
R
p
P
R p
R p
R p p
⎡
⎤
⎢
⎥
⎢
⎥
= ⎢
⎥
⎢
⎥
⎢
⎥
⎣
⎦
…
…
…
…
…
…
…
180
В отличие от матрицы
R в автокорреляционном методе матрица
P будет симметрической, но не теплицевой. Решение такой системы в общем виде требует 3
p операций.
Алгоритм решения уравнений линейного предсказания для автокор-реляционного метода. Система уравнений Юла – Волкера имеет вид (8.13), матрица коэффициентов
R является теплицевой и симметрической. Это по- зволяет найти решение за 2
p операций с помощью алгоритма Левинсо- на – Дарбина.
Данный алгоритм был предложен Левинсоном в 1948 г. и усовер- шенствован Дарбиным в 1960 г. Особенность алгоритма – его итеративный характер. В нем последовательно решается система уравнений вида (8.13) порядка 1, 2, ,
lp…
=
, причем решение системы порядка
l выражается че- рез решение системы порядка
1
l− .
Решение системы порядка
l будем обозначать через
(
)
( )
( )
( )
( )
1 2
,
,
,
lllllaaaa=
…
. На каждом шаге алгоритма вычисляется также ошибка предсказания
lE для решения системы
l -го порядка и вспомога- тельный коэффициент
lk .
Ниже приводится формальное описание алгоритма.
Начальные условия:
( )
(0)
0 0,
0 ,
0
lRaE=
=
= .
Итеративная процедура: при
1,
,
lp=
…
вычисляются
( ) ( )
1 1
1 1
lllilR l iR lkaEi−
−
⎛
⎞
=
⋅
− −
∑
⎜
⎟
−
⎜
⎟
=
⎝
⎠
,
( )
a
lllk= − ,
( )
( 1)
1
, 1 1
lllljjl jj laak a−
−
−
=
+ ⋅
≤ ≤ − ,
(
)
2 1
1
lllkEE−
=
−
На последнем шаге алгоритма при
l =
p получается искомое решение
(
)
( )
1 2
,
,
,
,
pppaE Ea aaa=
…
=
=
Пример. Проследим работу алгоритма Левинсона – Дарбина на пер- вых шагах.
Шаг 1:
(
)
1
(1)
2
(2)
1 1
(1) / (0),
(1) / (0),
(0)
(1) / (0).
kRRRRRaERR= −
=
=
−
181
Шаг 2:
(
) (
)
(
)
(
)
(
)
(
)
(
)
(
)
2 2
(2)
2
(2)
2 2
1 2
2 2
2 2
(1)
(0) (2)
(0)
(1) ,
(1) (0)
(1) (2)
(0)
(1) ,
(0)
(2)
(0)
(2) (0) 2
(1)
(0)
(1) .
R
R
k
R
R
R
R R
R R
a
R
R
R
R
R
R
E
R
R
R
R
=
−
−
=
−
−
=
−
+
−
−
При синтезе случайного процесса с помощью полученного фильтра необходимо предварительно рассмотреть вопрос об устойчивости. На практике отсутствие устойчивости проявляется в том, что моделирование с построенным фильтром не даст нужного результата, т.е. малым сигналам возбуждения на входе фильтра могут соответствовать большие выходные сигналы. Ниже рассмотрим некоторые вопросы, связанные с устойчиво- стью авторегрессионного фильтра и порождаемого им процесса.
Последовательность
{ }
( )
(1), (2), , ( ),
x n
x
x
x n
…
=
называется устойчи- вой, если
1
( )
( )
n
n
X z
x n z
∞
=
= ∑
абсолютно сходится внутри единичного кру- га, т.е. при
1
z
< . Фильтр ( )
H z называется устойчивым, если все полюсы
( )
H z лежат внутри единичного круга. Из этих определений следует, что если на вход фильтра поступает устойчивая последовательность и фильтр устойчив, то на выходе будет также устойчивая последовательность.
Действительно, мы имеем
( )
( )
( )
X z
Y z H z
=
< ∞ при всех
1
z
< то- гда и только тогда, когда ( )
Y z
< ∞ и ( )
H z
< ∞ при
1
z
< . Следовательно, для проверки устойчивости фильтра с передаточной функцией
( ) 1 ( )
H z
A z
=
нужно вычислить все корни полинома
1 1
( ) 1
p
p
A z
a z
a z
−
−
= −
− −
…
и убедиться в том, что они удовлетворяют ус- ловию
1, 1
i
z
i
p
<
≤ ≤ (для того чтобы привести ( )
A z к полиномиальному виду, можно умножить ( )
A z на
p
z ).
Процедура вычисления всех комплексных корней многочлена доста- точно трудоемка и на практике применяется редко. Однако если фильтр синтезируется по алгоритму Левинсона – Дарбина, то условием устойчи- вости будет выполнение на каждом шаге неравенства
1
l
k
< .
Итерационный алгоритм Левинсона – Дарбина реализуется в
MATLAB функцией levinson. Функция rlevinson решает обратную задачу – позволяет найти вектор отсчетов корреляционной функции сигнала по за- данным коэффициентам линейного предсказания.
Функция lpc реализует расчет коэффициентов линейного предсказа- ния автокорреляционным методом и является аналогом функции aryule.
Эти две функции различаются лишь MATLAB-кодом, используемым для вычисления оценки корреляционной матрицы. В итоге результаты совпа- дают с точностью до вычислительных погрешностей.
182
2. Задания и методические указания по выполнению работы
1. С помощью микрофонной гарнитуры введите в компьютер рече- вой сигнал (фамилию студента).
2. Самостоятельно изучите описание функций levinson, rlevinson и
lpc в среде MATLAB с помощью раздела меню Help.
3. Реализуйте с помощью функции levinson алгоритм
Левинсона – Дарбина для кадров речевого сигнала размером
180 – 240 отсчетов.
4. С помощью функции lpc рассчитайте коэффициенты линейного предсказания для речевого сигнала.
5. С помощью функции rlevinson найдите вектор отсчетов корреля- ционной функции речевого сигнала по полученным коэффициентам ли- нейного предсказания.
6. Повторите эксперименты для множества сегментов речевого сигнала.
7. Проанализируйте полученные результаты и представьте их в отчете.
183
Заключение
Несмотря на значительный прогресс в области теории и практики цифровой обработки речи и других сигналов усилия, направленные на ре- шение проблем ЦОС, не снижаются, а, напротив, возрастают. Успехи дос- тигнуты благодаря революции в области проектирования микроэлектрон- ных устройств с высочайшим уровнем интеграции. Этот прогресс наглядно демонстрирует, например, мобильная связь, причем значительная доля ус- пеха связана с развитием методов обработки и передачи сигналов.
В то же время остается нерешенной главная стратегическая задача в области обработки речевых сигналов – понимание речевого сообщения.
Надо учесть, что чисто механическое представление речевых сигналов яв- ляется односторонним, черно-белым. На его основе может быть решен лишь ограниченный круг простейших задач. Обучение человека речевому общению происходит в течение многих лет. Но ведь человек получает спо- собность к речевым коммуникациям от природы!
Речь – это способ передачи мысли. В ней существенную роль играют эмоции, интонация и контекст, а восприятие речи субъектом опирается на его предыдущий опыт. Эти богатейшие составляющие при обработке рече- вого сигнала в технической системе либо утрачиваются, либо пока не ис- пользуются. В диалоге речь сопровождается также общением по визуаль- ному каналу.
Хотя в области обработки речи уже немало сделано, предстоит сде- лать намного больше. В связи с перспективами коммерческого применения основные фирмы-разработчики «придерживают» информацию по новей- шим исследованиям. Надеемся, что достижения в области распознавания речи в ближайшее время будут множиться и станут доступны в учебной литературе.
184
Библиографический список*
1. Gaurang Kishor Parikh, B.E. The effect of noise on the spectrum of speech: thesis.– Texas Un-ty, 2002.
2. Haykin, S. Adaptive Filter Theory. – 4-th edition. – Prentice Hall, 2002.
3. Picone, J. Fundamentals of speech recognition / Department of Electri- cal and Computer Engineering, Mississippi State University. – Режим доступа: http://www.isip.msstate.edu/resources/cources/ece_8463.
4. Tokuda, K. Speech coding based on adaptive mel-cepstral analysis : proc.
ICASSP'94 / K. Tokuda, H. Matsumura, T. Kobayashi and S. Imai. – 1994. – Ре- жим доступа: http://ktlab.ics.nitech.ac.jp/tokuda/selected_pub/pdf/confe- rence/tokuda_icassp1994.pdf.
5. Gales, M. The Theory of Segmental Hidden Markov Models. – Cam- bridge University, 1993.
6. Mouly, M. The GSM System for Mobile Communications / M. Mouly,
M.B. Pautet. – 1992. – 702 p.
7. MATLAB 6.5 SP1/7.0 (R14). Пакет программ.
8. Mixed-Signal and DSP design techniques. Analog Devices (Walt Ke- ster) 2000.
9. Proakis, Y.G. Digital Communication. Mc Graw Hill. – 3-rd ed. – New
York, 1995.
10. Imai, S. Mel log spectrum approximation (MLSA) filter for speech synthesis : transactions of the IECE of Japan, February 1983 / S. Imai,
K. Sumita, and C. Furuichi.
11. Robinson, T.Speech Analysis. Lent Term 1998. – Режим доступа: http://svr-www.eng.cam.ac.uk/ajr/SpeechAnalysis/SpeechAnalysis.html.
12. Wiener, N. Extrapolation, interpolation and smoothing of stationary time series. – John Wiley, New York, 1949.
13. Сергиенко, А. Б. Алгоритмы адаптивной фильтрации : особенно- сти реализации в MATLAB // EXPonenta Pro. Математика в приложениях. –
2003. – № 1.
14. Новосельский, А. Интернет-публикация. – Режим доступа: http://dox.sbnet.ru:8082/serge/speech.ru/cplusp/6n96y2a.htm#vv.
15. Шульгин, В. И. Основы теории связи. Ч. 1. Теория и практика ко- дирования : учеб. пособие. – Харьков, 2005. – 194 с.
16. Введение в цифровую фильтрацию / под ред. Р. Богнера, А. Кон- стантинидиса. – М. : Мир, 1976. – 216 с.
17. Веденисов, Д. Синтез речи. – 2004. – Режим доступа: http://www.temator.ru/section/10/1.html.
* Приводится в авторской редакции.
185 18. Венцов, А. В. Современные модели восприятия речи : критиче- ский обзор / А. В. Венцов, В. Б. Касевич. – СПб : Изд-во Санкт-Петербург. ун-та, 1994. – 316 с.
19. Галунов, В. И. Аналитический обзор по проблеме кодирования речевых сигналов / В. И. Галунов, А. Б. Викторов. – Режим доступа: www.auditech.ru.
20. Гробман, М. З. Выделение скрытых периодичностей и формант- ный анализ речи. Распознавание образов : теория и приложения /
М. З. Гробман, В. И. Тумаркин. – М. : Наука, 1977.
21. Гудонавичюс, Р. В. Распознавание речевых сигналов по их струк- турным свойствам / Р. В. Гудонавичюс, П. П. Кемешис, А. Б. Читавичюс. –
Л. : Энергия, 1977. – 300 с.
22. Давыдов, А. В. Сигналы и линейные системы : Интернет- публикация. – Режим доступа: http://prodav.narod.ru/signals/index.html.
23. Дьяконов, В. П. MATLAB 6.5 SP1/7 +Simulink 5/6. Обработка сигналов и проектирование фильтров. – М. : СОЛОН-пресс, 2005. – 576 с. –
ISBN 5-98003-206-1.
24. Дэйвид, Г. Порядковые статистики. – М. : Наука, 1979. – 336 с.
25. Золотухин, И. П. Цифровые звуковые магнитофоны /
И. П. Золотухин, А. А. Изюмов, М. М. Райзман. – Томск : Радио и связь,
Том. отд., 1990. – 160 с. – (Массовая радиобиблиотека. Вып. 1153).
26. Каппелини, В. Цифровые фильтры и их применение / В. Каппе- лини, Дж. Константинидис, П. Эмилиани. – М. : Энергоатомиздат, 1983. –
360 с.
27. Кривошеев, В. И. Digital Signal Processing : курс лекций. – ННГУ,
2004. – Режим доступа: www.wl.unn.ru.
28. Рабинер, Л. Р. Скрытые марковские модели и их применение в избранных приложениях при распознавании речи // ТИИЭР. – 1989. –
Т. 77. – № 2.
29. Маркел, Дж. Линейное предсказание речи / Дж. Маркел,
А. Грей. – М. : Связь, 1980. –308 с.
30. Методы и стандарты кодирования и сжатия речи в цифровой те- лефонии : Интернет-публикации. – Режим доступа: http://dox.sbnet. ru:8082/serge/speech.ru/coder/.
31. Михайлов, В. Г. Измерение параметров речи / В. Г. Михайлов,
Л. В. Златоустова ; под ред. М. А. Сапожкова. – М : Радио и связь, 1979. –
416 с.
32. Назаров, М. В. Цифровая реализация устройств первичной обра- ботки речевых сигналов с линейным предсказанием / М. В. Назаров,
Р. В. Шафер // Тезисы докладов 11-го Всесоюзного семинара АРСО-11. –
Ереван, 1980.
186 33. Научно-исследовательская группа «Phrase research group» : Рече- вые технологии третьего тысячелетия : Интернет-публикация. – Режим доступа: http://phrase.ru/rus/npes.htm.
34. Огородников, А. Н. Материалы VIII Всерос. науч.-практ. конф. «На- учное творчество молодежи». – Томск : Изд-во Том. ун-та, 2004. – С. 52 – 53.
35. Пауков, Д. П. Импульсно-кодовая модуляция, использующая за- кон μ и Α. – ПМИ ДонНТУ, 2002. – Режим доступа: http://www.uran.do- netsk.ua/masters/2003/fvti/paukov/library/zakon2.htm.
36. Попов, В. И. Основы сотовой связи стандарта GSM / В. И. По- пов. – М.: Эко-Трендз, 2005. – 296 с. – ISBN 5-88405-068-2.
37. Потапова, Р. К. О типологических особенностях слога. Распозна- вание образов: теория и приложения / Р. К. Потапова. – М. : Наука, 1977. –
296 с.
38. Продеус, А. Н. Методы обработки акустических сигналов :
курс лекций. – Режим доступа: http://aprodeus.narod.ru.
39. Рабинер, Л. Р. Цифровая обработка речевых сигналов : пер. с англ. / Л. Р. Рабинер, Р. В. Шафер ; под ред. М. В. Назарова и Ю. Н. Про- хорова. – М. : Радио и связь, 1981. – 496 с.
40. Рабинер, Л. Р. Теория и применение цифровой обработки сигна- лов / Л. Р. Рабинер, В. Гоулд. – М. : Мир, 1978. – 848 с.
41. Иконин, С. Ю.Система автоматического распознавания речи
SPIRIT ASR Engine / С. Ю. Иконин, Д. В. Сарана // Цифровая обработка сигналов. – 2003. – № 3. – Режим доступа: www.spirit.ru и www.spiritdsp.com.
42. Секунов, Н. Ю. Обработка звука на РС / Н. Ю. Секунов. – СПб. :
БХВ – Петербург, 2001. – 1248 с. – ISBN 5-94157-037-6.
43. Сергиенко, А. Б. Цифровая обработка сигналов: учеб. для вузов /
А. Б. Сергиенко. – СПб. : Питер, 2003. – 604 с. – ISBN 5-318-00666-3.
44. Сорокин, В. Н. Элементы кодовой структуры речи. Распознава- ние образов: теория и приложения. – М. : Наука, 1977. – с. 42 – 60.
45. Смит, С. В. Научно-техническое руководство по цифровой обра- ботке сигналов / С. В. Смит ; пер. с англ. В. Н. Покровского, В. И. Силан- тьева. – СПб. : АВТЭКС, 2001. – 630 с.
46. Хайкин, С. Спектральный анализ радиолокационных мешающих отражений методом максимальной энтропии / С. Хайкин, Б. У. Карри,
С. Б. Кеслер // ТИИЭР. – №9. – 1982. – с. 51 – 62.
47. Хемминг, Р. В. Цифровые фильтры / Р. В. Хемминг. – М.: Недра,
1987. – 224 с.
48. Громаков, Ю. А. Сотовые системы подвижной радиосвязи. Тех- нологии электронных коммуникаций / Ю. А. Громаков. – М. : Эко-Трендз,
1994. – 302 с.
187
Оглавление Введение…………………………………………………………………...
3