А.Н. Огородников
ВЫБОР ИНТЕРВАЛОВ АНАЛИЗА СИГНАЛА ПРИ РАСПОЗНАВАНИИ РЕЧИ
Рассматривается проблема выбора интервалов анализа сигнала при распознавании речи. Исследуются преимущества и недостатки прямоугольного временного окна при кратковременном спектральном анализе. Предлагаются методы выбора интервалов анализа для речевых сигналов, которые позволяют устранить недостатки прямоугольного окна.
При распознавании речи обрабатываемый сигнал сначала разбивается на небольшие отрезки (интервалы анализа), которые затем классифицируются. Полученная последовательность отрезков речевого сигнала может затем подвергаться дальнейшей обработке, например, с целью получения последовательности слов. Интервалы анализа обычно имеют одинаковую длину, а для их классификации переходят к спектральному представлению сигнала. Перед нахождением спектра сигнал подвергается преобразованию с помощью непрямоугольного в общем случае временного окна (как правило, Хэмминга или Ханна). Такое преобразование искажает сигнал и может затруднить его классификацию. Использование прямоугольного окна, не искажающего сигнал, также имеет свои недостатки. Однако возможно выбирать интервалы анализа таким образом, что недостатки прямоугольного окна устраняются и при этом сохраняются его полезные свойства.
ТРЕБОВАНИЯ К ДЛИНЕ ИНТЕРВАЛОВ АНАЛИЗА
Речевой сигнал можно представить как функцию времени р(/), которая выражает изменение во времени какой-либо физической величины, например звукового давления или напряжения.
Согласно общепринятой акустической теории ре-чеобразования [1], преобразование Лапласа от функции р(/) звукового давления можно представить в виде произведения
Р(8) = Б(5) • Т(5),
где Б(5) есть преобразование Лапласа для функции источника звука, а Т(5) - передаточная функция речевого тракта. Передаточная функция Т(5) является меро-морфной функцией и полностью определяется бесконечным числом комплексно-сопряженных полюсов и нулей.
Функция 5(/) источника звука представляет собой изменение во времени звукового давления, вызванное источником звука. Источник звука может быть либо голосовым, либо шумовым. Голосовой источник обусловлен колебанием голосовых связок, функция 5(/) в этом случае является квазипериодической. Среднее значение Т0 периода колебаний голосовых связок называется периодом основного тона речи, а обратная ей величина ^0=1/Т0 представляет собой частоту основного тона.
Отклонения периода основного тона от среднего значения можно классифицировать следующим образом:
1) случайные отклонения (джиттер-эффект [2]). Эти отклонения невелики и возможны как в сторону увеличения, так и в сторону уменьшения периода основного тона. Случайные отклонения не несут какой-либо смысловой нагрузки;
2) вибрато. Вибрато есть частотная модуляция колебаний голосовых связок и проявляется в периодическом или квазипериодическом изменении периода
основного тона. Вибрато характерно для пения, и целью его использования является придание голосу красоты. Как правило, вибрато не несет какой-либо смысловой нагрузки;
3) плавное увеличение или уменьшение периода основного тона. Этот тип отклонений является частью интонационного рисунка речи. Обычно увеличение частоты основного тона вызывается увеличением подсвязочного давления, вследствие чего увеличивается также и громкость речи, и наоборот, хотя увеличение частоты основного тона может происходить и без увеличением подсвязочного давления. Отслеживание изменения частоты основного тона совместно с интенсивностью звука является важным при определении границ слогов, фонетических слов и фонетических синтагм и необходимым для определения словесных и синтагматических ударений и интонации, которая указывает на общую коммуникативную цель говорящего.
Отслеживание изменения частоты основного тона требует анализа речевого сигнала на достаточно небольших временных интервалах.
Небольшие интервалы анализа (порядка одного периода основного тона) требуются и при отслеживании сдвига формант, имеющего место во время переходных процессов в речевом тракте. Например, слова [mama] и [nana] различаются человеком только за счет сдвига формант в начале произнесения первого и второго [a], в то время как сами по себе [m] и [n] практически не различаются [1].
ДИСКРЕТИЗАЦИЯ РЕЧЕВЫХ СИГНАЛОВ
Если в процессе произнесения данного речевого сигнала ему предшествовали и за ним следовали достаточно большие паузы, функцию s(t) источника звука можно считать финитной. В этом случае изображение S(s) будет целой функцией, а P(s) - мероморфной, из чего следует нефинитность функции p(t).
В реальности наблюдение функции p(t) ограничивают каким-либо интервалом T1, превращая ее в финитную. Преобразование Лапласа P^(s) такой финитной функции prfj(t) будет целой функцией, т.е. информация о точном нахождении полюсов функции P(s), а следовательно, и передаточной функции T(s), будет потеряна.
Пусть отсчет времени для функции p^(t) ведется от момента времени t=0. Тогда функцию p^(t) можно разложить в ряд Фурье на интервале [0, T1)
2п
ГО i—nt
p<i>(t) = £ cteT .
n=-ro
Здесь предполагается, что это разложение существует.
Перед обработкой на компьютере речевой сигнал p^(t) дискретизуется с шагом At. Соответствующее
сигналу рф(0 дискретное колебание рд(0 можно представить выражением [3]
N-1
Рд (/) = А/ X Рф (кА )8(/ - кА1),
к=0
где ^Т^А/ - количество отсчетов сигнала Рф(/), 5(/) -дельта-функция Дирака.
Коэффициенты сд ряда Фурье дискретного сигнала рд(/) выражаются через коэффициенты сф ряда Фурье сигнала Рф(/) следующим образом:
= У сф
/ -- п
к=-го
-Ж
(1)
а само значение коэффициентов сп дается дискретным преобразованием Фурье (ДПФ)
1 N —1
=N У рф(кы)е
2п ,
-г—пк
N
(2)
к=0
(3')
с ? = сф
г-, г-.
п = 0, ±1,
±'т—1
с—
—N/2 = %/2
Пусть сф = Яе{ сф } + г 1т{сф }, сфп = Яе{сф } — г 1т{сф },
тогда, поскольку
- N/2
/2
= 2Яе{сф ,2}.
(5)
Из равенства (5) можно сделать вывод, что синусная составляющая гармоники с номером N/2 сигнала Рф(/) при дискретизации теряется, что естественно, так как функцию 8ш(п//А/) нельзя корректно дискретизовать с шагом А/. С другой стороны, косинусная составляющая корректно дискретизуется, но из-за эффекта наложения спектров происходит ее удвоение.
Из (5) следует, что в случае четного N идеальный синтезирующий фильтр на выходе цифроаналогового преобразователя должен уменьшать в 2 раза амплитуды гармоник с номерами ±N/2.
В связи с вышеизложенным удобно ввести коэффициенты сПд, связанные с коэффициентами сд ряда Фурье дискретного сигнала Рд(/) следующим образом:
0 ±1 ■ " N-1"
для п = 0, ±1,
с? =
Интервал разложения для Рд(/) равен [0, Т1).
Предполагается, что в дискретном колебании Рд(/) на выходе цифроаналогового преобразователя будут отфильтрованы гармоники с круговыми частотами выше п/А/, и на выходе синтезирующего фильтра будет восстановленный сигнал Рв(/) [4]. В случае идеального синтезирующего фильтра спектр сигнала Рв(/) будет состоять из коэффициентов с^, п = 0, ±1, ..., ±[N/2], где квадратные скобки означают взятие целой части числа.
Если для коэффициентов сф ряда Фурье сигнала Рф(/) выполняется условие
сф = 0 для | п | > N/ 2 , (3)
то выражение (1) преобразуется к виду
сд = сф, п = 0, ±1, ..., ±[N/2], (4)
т.е. сигналы Рф(/) и Рв(/) будут совпадать.
Если условие (3) не выполняется, то происходит эффект наложения спектров [5]. Чтобы избежать этого эффекта, перед дискретизацией сигнала Рф(/) производится фильтрация гармоник с круговыми частотами, большими п/А/. Однако, поскольку аналоговые фильтры нижних частот не являются идеальными, то эффект наложения спектров может иметь место и в этом случае.
В случае четного N из условия (3) следует с-N /2 = cN /2 = 0 . Однако на практике при выполнении ДПФ часто оказывается, что cд_N /2 Ф 0 и сN /2 Ф 0. Поскольку в этом случае условие (3) уже не выполняется, введем более слабое условие
сф = 0 для|п |> N/2.
Из (1) и (3') в случае четного N следует
сд ^
— для п = ±— при четном N. 2 2
(6)
Таким образом, из выполнения условия (3') будет следовать
с? =
п = 0, ±1,
N — 1
N
(7)
Яе{сф} для п = ± — при четном N.
Поскольку синусная составляющая гармоники с номером N/2 сигнала Рф(/) потеряна, то можно считать, что она была равна нулю. Тогда из (3') и (7) получаем
сд = сф, п = 0, ±1, ..., ±[N/2], (4')
и с учетом вышесказанного сигналы Рф(/) и Рв(/) будут совпадать и в этом случае.
Если абсолютные значения коэффициентов сф для |п|>М2 достаточно малы, так что выполняется приближенное равенство
сд * сф, п = 0, ±1, ..., ±[N/2], то удобно исходным считать сигнал
[ N/2]
рф (1) = У
п=—[ N / 2]
2п
г—п
сде Т
(8)
Для сигнала (8) будут автоматически выполняться условия (3') и (4'). В дальнейших рассуждениях будем считать, что для сигнала Рф(/) условия (3') и (4') выполняются.
ПРОБЛЕМА ВЫБОРА ИНТЕРВАЛОВ АНАЛИЗА ДЛЯ РЕЧЕВЫХ СИГНАЛОВ
При распознавании речи большой интерес представляют спектральные свойства речевого сигнала Рф(/). Для нахождения спектра дискретного сигнала Рд(/) используется ДПФ. В этом случае мы в соответствии с равенством (4') одновременно получаем и спектр сигнала Рф(/).
Но поскольку сигнал Рф(/) в общем случае имеет достаточно большую длительность и представляет собой реализацию некоторого количества фонем, то
интересны, прежде всего, спектральные свойства не всего сигнала Рф(/), а его локальных участков.
В связи с этим возникает проблема сегментации речи: перед тем как определять, к какому классу относится тот или иной участок речевого сигнала, речевой сигнал необходимо сегментировать на эти участки.
Сегментация речевого сигнала неоднозначна и зависит от используемых критериев. Общепринятым подходом является разбиение сигнала на пересекающиеся или непересекающиеся интервалы времени фиксированной длины в диапазоне от 2 до 50 мс, исходя из предположения, что на интервале данной длины изменение характеристик речевого тракта достаточно незначительно.
Перед выполнением спектрального анализа из сигнала Рд(/) с помощью временного окна w(m) выделяется кратковременный сигнал
Фурье в тригонометрической форме:
N-1
pдт (t) = At У w(k - m)pi (kAt)S(t - kAt),
(9)
k=0
M
M-1
nk
У w(к)pф ((m + k)At)e M ,(10)
k=0
wH(m) = ■
;-(i-a)cos
2nm
M
0 при других m.
0 < m < M -1,
p^«(t) = p4> (t)
a-(i-a)cos| t
(11)
[ N/2]
pф (t) = ^ + У An c0s| + 0n I . (12)
2 n=1 V T2
Из (11) и (12) видно, что преобразование сигнала Pф(t) с помощью обобщенного окна Хэмминга эквивалентно амплитудной модуляции с круговой частотой 2п/Т2 каждой гармоники сигнала Pф(t), т.е. у каждой гармоники с круговой частотой 2пп/Т2 появляются копии с круговыми частотами 2n(n - 1)/T2 и 2п(п+1)/Т2 и амплитудой (1 - a)An/2. Эти копии накладываются на соседние гармоники, искажая их.
Это явление можно проинтерпретировать и другим образом, а именно: при преобразовании сигнала с помощью обобщенного окна Хэмминга происходит размытие спектра по частотной шкале. Из приведенных рассуждений видно, что при использовании временных окон, отличных от прямоугольного, будет происходить то или иное искажение спектра, которое можно представлять себе в виде его некоторого размытия.
Размытие спектра является недостатком, которого желательно избегать, так как при этом могут потеряться или притупиться формантные максимумы, что сделает задачу классификации данного участка речевого сигнала более сложной.
равный нулю вне интервала [/0, /0+Т2), /0=тА/, Т2=МА/. ДПФ этого сигнала на интервале разложения
[/0, /0+Т2) есть
где стп - п-й коэффициент ряда Фурье сигнала Рдт(/).
2п
-г—пт
Множитель е м в (10) представляет собой линейное изменение фазы, связанное с задержкой сигнала на т отсчетов. Обычно этот множитель не учитывается, поскольку положение сигнала Рдт(/) на временной оси, как правило, неважно.
Наиболее распространенными типами окон являются прямоугольное и обобщенное окна Хэмминга [6, 7].
Окна, отличные от прямоугольного, дополнительно искажают сигнал Рдт(/), поэтому искажается и ряд Фурье (10) в сравнении с рядом для прямоугольного окна.
Рассмотрим для примера искажение, вносимое м-точечным обобщенным окном Хэмминга
Пусть на интервал длиной Т2=МА/, выделяемый окном ^н(т), приходится ровно один период колебания голосовых связок, взятый из некоторого речевого сигнала. Для выделенного таким образом исходного сигнала Рф(/) длительностью Т2 ряд Фурье (2) совпадает с рядом Фурье (10) в случае использования м-точечного прямоугольного окна.
В случае же использования обобщенного окна Хэмминга ряд Фурье (10) будет совпадать с рядом Фурье сигнала
8000
6000
4000
2000
0
8000
6000
4000
2000
0
I iI II I i II ■
16
21
26
I ‘ ‘ ‘ I
31 36
Jjjj
16
21
26
I ■ I
31 36
Чтобы отследить отличия спектров у сигналов Рф(/) и Рфи(/), удобно представить сигнал Рф(/) рядом
Рис. 1. Сравнение ДПФ с использованием прямоугольного окна и окна Ханна: а) дискретный сигнал, для которого выполняется ДПФ, по оси абсцисс отложены номера отсчетов; б) амплитудный спектр для прямоугольного окна; в) амплитудный спектр для окна Ханна. Спектральные коэффициенты пронумерованы начиная с единицы
На рис. 1, б и в представлены для сравнения спектры, полученные с использованием прямоугольного окна и окна Ханна (обобщенного окна Хэмминга с параметром а = 0,5) на примере дискретного сигнала, представляющего собой один период речевого сигнала, взятый из реализации фонемы /а/, произнесенной
б
6
в
6
женским голосом (показаны только первые 40 спектральных коэффициентов, модули остальных коэффициентов пренебрежимо малы). На рис. 1,б отчетливо видны формантные максимумы, соответствующие первым двум формантам: ^ находится около гармоники п=3 (на рисунке она имеет номер 4) с частотой 723 Гц, Е2 находится около гармоники п=5 с частотой 1205 Гц. На рис. 1,в видно, что из-за размытия спектра формантный максимум, соответствующий первой форманте, исчез.
В [6] отмечается, что частотное разрешение у прямоугольного окна больше, чем у окна Хэмминга, однако прямоугольное окно редко используется при спектральном анализе речи из-за того, что преобразование Фурье
хт (е™) = X ™(к - т)Рф (кА/)е
-г'юк
к=-ж
дискретного сигнала Рдт(/) в случае прямоугольного окна имеет более изрезанный характер между гармониками, чем, например, преобразование Фурье с окном Хэмминга [6].
Но если заданы значения Хт(ет) в точках ю=2пп/М, п=0, ..., [М/2], т.е. гармоники, то значения Хт(егю) для других значений ю не несут информации, так как могут быть однозначно вычислены. Поэтому данное свойство прямоугольного окна не должно препятствовать его использованию при спектральном анализе речи.
Рассмотрим еще один недостаток прямоугольного окна, который не упоминается в литературе.
В случае прямоугольного окна выражение (9) приводится к виду
т +М-1
Рдт (/) = А/ X Рф (к А/ )§(/ - к А/),
а выражение (10) - к виду
,2п
-г—пт 1
сд = е М
тп
М X Рф((т+к )А/ )е
2п
г—п М
Для нас интересен случай, когда М<^ где N - количество отсчетов сигнала Рф(/). Найдем соответствие между спектрами дискретного сигнала Рдт(/) и сигнала
Р (/) = {Рф(/), /о - * </о + T2,
фт [0 в противном случае,
которому соответствует дискретное колебание Рдт(/).
Для этого разложим сигнал Рфт(/) в ряд Фурье на интервале [/0, /0+Т2):
.(/) = X
2пп
г------/
Искомое соответствие между спектрами Рдт(/) и Рфт(/) определяется выражением (1), которое можно переписать в следующем виде:
= X
(13)
к=-ж
В выражении (13) для нас важны только коэффициенты с номерами п=0, ±1, ..., ±[М/2], так как остальные коэффициенты получаются периодическим
повторением с шагом М коэффициентов из данного диапазона.
В дальнейшем при анализе скорости убывания коэффициентов т на производную функции Рфт(/) накладывается достаточно слабое ограничение абсолютной интегрируемости, которое следует хотя бы из условия (3').
Если периодически продолженная функция Рфт(/) терпит разрывы первого рода в точках /0+пТ2, пе.Ъ, т.е. Рфт(/0) ФРфт(/0+Т2-0), то ряд Фурье сигнала Рфт(/) будет иметь бесконечное число ненулевых коэффициентов ^п , которые при п ^ ж будут бесконечно малыми первого порядка [8], т.е. будут сравнительно медленно убывать.
Анализируя выражение (13), приходим к выводу, что в случае разрыва функции Рфт(/) на краях интервала [/0, /0+Т2) спектр дискретного сигнала Рдт(/) будет искажен или зашумлен медленно убывающими хвостами спектра континуального сигнала Рфт(/).
Абсолютная величина разрыва
с =| Рфт (/0 + Т2 - 0) - Рфт (/0)| также непосредственно влияет на величину коэффициентов стфп . Если представить стфп в виде
сфп = 2 (п - гЬп ),
то абсолютные значения коэффициентов Ьп будут асимптотически эквивалентны с/пп при п ^ ж [8], из чего следует, что чем больший разрыв терпит функция Рфт(/) на краях интервала [/0, /0+Т2), тем более зашумленным окажется спектр сигнала Рдт(/).
Описанный эффект зашумления спектра дискретного сигнала является недостатком, которого необходимо избегать. Из вышеизложенного виден и способ уменьшения зашумления спектра сигнала Рдт(/): необходимо уменьшить величину разрыва функции Рфт(/) на краях интервала анализа, доведя ее в идеале до нуля. Этого можно достичь путем выбора подходящего интервала анализа.
На скорость убывания хвостов спектра функции Рфт( ) влияют также и разывы ее производных, а именно: если у функции Рфт(/) к-1 производных непрерывны, а к-я производная абсолютно интегрируема, то для модулей коэффициентов стфп ряда Фурье функции Рфт( ) имеет место соотношение
Иш п
= 0.
Выбирая интервал анализа таким образом, чтобы выполнялись условия
Р$1«0 + 0) = Р(й(/0 + Т2 -0) (14)
для к<к0, и предполагая, что к0-1 производных у функции Рфт(/) непрерывны на интервале [/0, /0+Т2), а к0-я ее производная абсолютно интегрируема на этом же интервале, что будет выполняться, например, при выполнении условия (3'), мы получаем, что последовательность |сфп | при п ^ ж будет бесконечно малой
более высокого порядка малости, чем 1/ пк°.
п=-ж
Таким образом, обобщая, можно сказать, что для минимизации эффекта зашумления спектра дискретного сигнала рдт(0 интервал анализа надо выбирать таким образом, чтобы выполнялись условия (14) для как можно большего к. В пределе, если условия (14) выполняются для любого к, мы получаем скорость убывания хвостов спектра функции рфт(() быстрее любой полиномиальной функции.
Чтобы определить абсолютную величину с(М) разрыва функции рфт(() на краях интервала [0, ¿0+Т2), имея в наличии только дискретную функцию рд(0, предположим непрерывность функции Рф(0. Отсюда следует
Рфт (¿0 + Т2 - 0) = Рф (¿0 + Т2) = Рд (Ч + Т2 ) и с(М) = |рд ((т + М)Д/) - рд (тД?)|. (15)
МЕТОД ВЫБОРА ИНТЕРВАЛА АНАЛИЗА, ДАЮЩИЙ ТОЧНОЕ РЕШЕНИЕ
Будем считать, что начало ¿0=тД интервала анализа задано. Требуется найти его длину Т2=МД, исходя из критерия минимума зашумления спектра. На длину накладываются ограничения Тт1П<Т2<Ттах.
Ограничение Тт1П<Т2 предполагает, что отрезок речевого сигнала с длительностью, меньшей чем Ттш, недостаточно хорошо отражает свойства кратковременно установившегося процесса в речевом тракте. В свою очередь, ограничение Т2<Ттах предполагает, что на отрезке речевого сигнала с длительностью, большей чем Ттах, процесс, происходящий в речевом тракте, успевает достаточно сильно измениться, а для распознавания речи наибольший интерес предсталяют именно кратковременно установившиеся процессы.
Минимизация зашумления спектра заключается в
минимизации модулей |сфп | спектральных коэффициентов сигнала рфт(/) для п>М/2. В качестве минимизируемой величины естественно выбрать среднюю на интервале Т2 мощность
ж 2
Р (М) = 2 X \сЦ
п=[М/2]+1
сигнала рфт( ) вне основной полосы с граничной круговой частотой п/Д. Средняя мощность Р(М) зависит и от начала 0 интервала анализа, который мы полагаем фиксированным.
При выполнении условия (3') можно найти сколь угодно точное решение этой задачи. Действительно, поскольку функция рф( ) представима конечной суммой известных гармоник, мы можем восстановить ее значения в любой точке, а следовательно, разложить ее в ряд Фурье на интервале [/0, ¿0+Т2), после чего вычислить Р (М) на основании равенства
1 ¿0 +Т2 [М /2] 2
Р (М) = — I р2фт (0* - X \сФфп\ , (16)
Т
-[ M/2]
ется вычислять значения подынтегральных функций в различных точках. Вычисление значений подынтегральных функций в каждой точке потребует суммирования N синусоидальных и косинусоидальных составляющих ряда Фурье функции рф(/). Количество элементарных операций, которое нужно выполнить для вычисления Р (М), зависит от метода и точности численного интегрирования. Пусть для приближенного вычисления интеграла на интервале Д необходимо вычислить значения рф(/) в С точках. Тогда на интервале Т2 необходимо вычислить значения рф(/) в СМ точках, и общее число элементарных операций, нужное для вычисления Р(М), составит СЩМ2+2М).
Для нахождения наилучшего интервала анализа понадобится вычислить значения Р(М) для М=Мтш,
• • •, ^^тах, где ^Мш1п=Тт1п/Д ¿, ^Мтах=Ттах/Д ¿.
Пусть Q = Мтах - Мт1п. Тогда трудоемкость нахождения наилучшего интервала анализа при численном вычислении интегралов составит
T (Q) = £ CN (M 2 + 2M) =
где
= CN (aQ3 + bQ2 + cQ + d),
а = 1/3, b = Mmin + 3/2,
(17)
C = Mmin + 3Mmin + 7/6 ,
d = Mmin + 2Mmin-Трудоемкость T(Q), даваемая выражением (17), является кубической относительно Q, квадратичной относительно Mmin и линейной относительно N. Также необходимо учитывать, что для разложения функции Рф(0 в ряд Фурье на интервале [0, T1) потребуется выполнить ДПФ, которое имеет в общем случае квадратичную трудоемкость относительно N.
Эксперимент показывает, что скорость вычисления значений P (M) с помощью численного интегрирования является крайне низкой, что делает данный способ нахождения наилучшего интервала анализа непригодным для практического применения. Кроме того, из-за вычислительных погрешностей, возникающих при численном интегрировании, полученный результат может отличаться от точного решения.
Интегралы в выражении (16) можно вычислить и аналитически. Для этого представим функцию Рф(0 в виде ряда Фурье в тригонометрической форме:
а0 [ х-?] . . (2пп
Рф(t) = — + £ An sml------1 + {
n=1
T1
Тогда средняя на интервале [/0, /0+Т2) мощность сигнала рфт( ) будет равна
¿0 +Т2
j р2Фш(t)dt =
1 [N/2]
+2 £ A
2 п=1
2
N
nM
которое выполняется ввиду полноты тригонометрической системы.
Для вычисления Р (М) по формуле (16) нам понадобится вычислить значения М+2 интегралов. При использовании численного интегрирования потребу-
а [ N/2] 1
т £ а. -
2 п=1 п
cos
, „ ш + M . - cos l 2п------п + 9
N
2п—п + 9п |-
N
1[ N^,2 1 8 £ а. ->
8 п=1 п
0
X
. „ т _ ) . т + М
Б1п21 2л—п + 9п |-81п2| 2л---------------п + (
N ) \ N
1 [N/2] N/2]
+4 X ^п X А
4 п=1 к=1
к ф п
п - к
. . „ т + М _ , ч _ _ .
2л (п - к) + 9п-9к |-
т
- б1п| 2л ^ (п - к) + 9п-9к
1 I • Г ~ т + М _ ,ч _ _ ,
Б1п| 2л——— (п + к) + 9п +9к |-
п + к
N
-Н 2лN(п + к) + 9п +9к
(18)
Re{cфn } = Т- ] рфт () 008 п^ , (19)
где
2 ¿0 ¿0 +Т2
1 Г /ч- 2л
1т{сфп } = — I рфт 0)81п — М .
гр I •* ууш ' ^ /Т-Т
Т2 г Т 2
¿0
(20)
Интеграл (19) для пф0 можно представить в виде
¿0 +Т2
N° л • О N
X Л 81п9к + —>
к=1 2л
кМ=nN
[ N/2]
к=1 кМ ф п^
-0081 2л 1
кМ - nN т + М
ж-2 ] . 1 | Г _ т тп
X А------------^ 0081 2л—к - 2л-------------------------+(
N
~\г--к - 2^ТТ + 9к ||+ X Ак:
N М Л к=
М
[ N/2]
ооб| 2лтк + 2л— + 9к
кМ + nN I I N М
-00Б | 2л
т + М
к + 2л
тп
N М
Интеграл (20) для пф0 можно представить в виде
(21)
¿0 +Т2
[ Х"°] л 9 N
X Ак008 9к+—> к=1 2л
кМ=пы
[ N/2] ,
X А '
к кМ - п^
кМ фпЫ
бш| 2л
т + М
N
к- 2л
тп
М
.1 т , „ тп „п V
-Н 2лТтк-2^ТТ + 9к 1г- X Ак-
N
М
[ N/2]
=1
1
кМ + nN
. , „ т + М , „ тп ,
бш| 2л----------к + 2л-------+9
N
М
. , - т , „ тп .
-бш| 2л—к + 2л-----------+(
N
М
В случае п=0 коэффициент стф,0 равен
1 ¿0 +Т2 а лг [ N/2] 1
ф 1 г - ч, -0 N 1 ^ -1 ,1
Ст,0 = ~ | Рфт (0 ^ = -2" + 2лМ ^ кЛ
Для вычисления средней мощности сигнала Рфт(() в основной полосе с граничной круговой частотой п/Д необходимо вычислить квадраты модулей коэффициентов стфп , которые равны
1сфп Г = Ке2{сфп } + 1т2{сфп },
т
х ^ ооб| 2л^к + 9к | -ооб| 2л
т + М
N
к + 9к |^. (23)
Трудоемкость вычисления Р^М) по формулам (18), (21) - (23) равна N2/4+2N+NM, и трудоемкость нахождения наилучшего интервала анализа при аналитическом вычислении интегралов составит
Мтах
Т (О) = X (N2/4 + 2 N + NM) =
М =М т1п
= -Q2 + bQ + с, (24)
где Q = Мтах — Мт1п
и - = N/2,
Ь = N2 / 4 + N(Мт1п + 5/4),
с = N2/4 + N (Мт1п + 2).
Трудоемкость (24) является квадратичной относительно Q и линейной относительно Мт1п, т.е. по этим параметрам лучше, чем трудоемкость (17). Относительно N трудоемкость (24) является квадратичной, т.е. теоретически хуже, чем трудоемкость (17), однако необходимо учитывать, что и при численном, и при аналитическом вычислении интегралов потребуется нахождение ДПФ функции Рф(0, а оно имеет квадратичную относительно N трудоемкость.
Экспериментально было установлено, что скорость вычисления значений Р(М) с помощью формул (18), (21) - (23) существенно выше, чем с помощью численного интегрирования. Тем не менее, эта скорость остается достаточно низкой, и данный способ нахождения наилучшего интервала анализа также является малопригодным для практического применения.
Кроме этого, вычисления по формулам (18), (21) -(23) даже при использовании чисел расширенной точности стандарта 1ЕЕЕ-854 дают неточный результат. Например, были получены отрицательные значения Р(М), составляющие около 2,8% от вычисленной средней мощности сигнала рфт( ) на интервале [¿0, ¿0+Т2). В таких случаях вычисление Р (М) с помощью численного интегрирования при достаточно малом шаге интегрирования дает более точный резуль-
X
тат. Это не позволяет говорить, что вычисления по формулам (18), (21) - (23) в общем случае дают более точный результат, чем результат, полученный с помощью численного интегрирования.
Таким образом, несмотря на то, что метод минимизации Р(М) теоретически дает точное решение задачи выбора интервала анализа по данному критерию, на практике получить точное решение может оказаться затруднительным.
На рис. 2,6 для иллюстрации эффекта зашумления спектра приведена функция Р (М) зашумления спектра дискретного сигнала, изображенного на рис. 2,а. Как видно из примера на рис. 2, зашумление спектра дискретного сигнала увеличивается по мере увеличения разрыва функции на краях интервала анализа. Также видно, что чем больше длина интервала анализа, тем меньше зашумление спектра при той же величине разрыва функции на краях интервала. Следствием этого является то, что точки максимумов разрыва функции и зашумления спектра не совпадают. Однако точки минимумов разрыва функции и зашумления спектра совпадают, функция Р (М) на области определения М=1,...,100 достигает минимума в трех точках: М=1, М=50 и М=100 (флуктуации при больших значениях Р(М) не учитываются). Функция Р (М) достигает наименьшего значения в точке М=100, Р^100)и0, поскольку спектр исходного континуального сигнала для М=100 состоит из одной гармоники.
5,0-10°
3,8-10“
2,510°
1,3-10°
|||......... ..........и
ІШд,
интервале [t0, t0+T2)
і -0 ' -*2
Рф (M) = Т i рфm {t)dt •
llw.
1 11 21 31 41 51 61 71 81 91
0,018
0,014
0,009
0,005
1 11 21 31 41 51 61 71 81 91
Рис. 2. Пример зашумления спектра дискретного сигнала:
а) дискретный сигнал - один период дискретизованной синусоиды, N=100; 6) функция Р/М) для т=0 иМ=1,...,100
Минимизация Р(М) есть минимизация абсолютного зашумления спектра, но можно минимизировать и относительное зашумление спектра, т.е. минимизировать величину
^ (м )=,
‘ Рф (М)
где Рф(М) есть средняя мощность сигнала Рфт(ґ) на
t0 +T2
1 11 21 31 41 51 61 71 81 91
Рис. 3. Сравнение абсолютного и относительного зашумления спектра дискретного сигнала: а) дискретный сигнал -один период дискретизованной функции (25), jV=100;
б) функция P (M) для т=0 и M=1,...,100; в) функция R(M) для т=0 и M=1,...,100
Сравнение функций Pt(M) и Rt(M) приведено на рис. 3, б и в на примере дискретной функции (рис. 3,а), полученной при дискретизации одного периода континуальной функции
Рф (t) = 16383(sin rajt + cos2ra1t), (25)
где ю1=2п/7’1. В примере на рис. 3 функции Pt(M) и Rt(M очень похожи с точностью до масштабного множителя. Точки минимумов у функций P(M) и Rt(M на области определения M=1,...,100 совпадают. Этими точками являются точки M=1, M=8, M=42, M=50 и M=100. Эти же точки являются точками минимумов разрыва функции на краях интервала анализа.
Сходство функций Pt(M) и Rt(M) будет наблюдаться в тех случаях, когда функция Рф(М) является слабо осциллирующей. Для речевых сигналов это обычно выполняется. Исключением является случай, когда функция Рфт(0 близка к нулю в точке t0. В этом случае функция Рф(М) будет близка к нулю при небольших M с последующей стабилизацией на некотором среднем уровне при увеличении M. При таком поведении функции Рф(М) функции P(M) и Rt(M) будут похожи с точностью до фиксированного масштабного множителя, за исключением небольших значений M, при которых значения функции Rt(M) будут резко завышены. Однако, учитывая ограничение Mmm<M и предполагая, что на интервале длительностью Mmin функция Рф(М) успевает стабилизироваться, можно в первом приближении считать, что минимизация абсолютного и относительного зашумления спектра будет давать одинаковый результат.
3,510
б
2,610
1,810
8,810
0
в
б
0
МЕТОД МИНИМИЗАЦИИ ВЫСОКОЧАСТОТНЫХ ГАРМОНИК ДИСКРЕТНОГО СИГНАЛА
При спектральном анализе речи рассматриваемый диапазон частот обычно ограничивают сверху значением 8 - 10 кГц [2], поскольку в этом диапазоне заключается основная доля энергии речевого сигнала, а энергия более высокочастотных гармоник пренебрежимо мала. При выполнении ДПФ на высокочастотные гармоники дискретного сигнала рдт(ґ) будут накладываться хвосты спектра исходного континуального сигнала Рфт(ґ). Можно предположить, что чем больше средняя мощность Р(М хвостов спектра сигнала Рфт(0, тем большей окажется средняя мощность Ри(М высокочастотных гармоник соответствующего ему дискретного сигнала рдт(ґ). В общем случае это предположение выполняться не будет, так как при наложении хвостов спектра континуального сигнала на высокочастотные гармоники дискретного сигнала будет происходить как их увеличение, так и их уменьшение, в зависимости от фаз суммируемых гармоник. Однако в среднем поведение функций Р(М и Рн(М будет схожим, и минимизацию функции Ри(М можно использовать для приближенного поиска интервалов анализа, на которых зашумление спектра дискретного сигнала минимально.
Пусть / - частота, ниже которой заключена основная доля энергии речевого сигнала, тогда средняя мощность высокочастотных гармоник дискретного сигнала рдт(і) будет равна
функция Ри(М) достигает наименьшего значения Ри(100)«0,053.
пшах
Ри (М) = 2 £
(26)
где пш
д
тп
п_пшт
= [ /¡М Аґ ], пшах=[М/2], а коэффициенты с
задаются аналогично (6).
Пример функции Ри(М) для /,= 15000 Гц, Аґ= 1/44100 с и дискретного сигнала, изображенного на рис. 2,а, приведен на рис. 4.
5,0-106 3,8-106 2,5106 1,3106
ЇЛИ
1 11 21 31 41 51 61 71 81 91
Рис. 4. Пример функции Р*(М) для дискретного сигнала, изображенного на рис. 2,а, у^=15000 Гц, Дг=1/44100 с, т=0, М=1,...,100
Сравнивая функции на рис. 2,б и 4, можно обнаружить их сходство, за возможным исключением небольших значений М (М=1,...,30), при которых функция Ри(М) ведет себя достаточно неустойчиво. Как и следовало ожидать, значения функции Рь(М) оказываются несколько заниженными в сравнении со значениями функции Р(М). На области определения М=40,...,100 функция Ри(М) достигает минимума в двух точках: М=51 и М=100. Второе решение (М=100) совпадает с решением, полученным с помощью точного метода выбора интервала анализа. В этом случае
105
104
103
102
10
1
1 11 21 31 41 51 61 71 81 91
Рис. 5. Сравнение спектров дискретного, изображенного на рис. 2,а, и исходного континуального сигналов для интервала анализа длинойМ=50 отсчетов: а) амплитудный спектр дискретного сигнала; 6) амплитудный спектр исходного континуального сигнала, показаны только первые 100 коэффициентов
Решение М=51 очень близко к точному решению М=50, но все же отличается от него. Несмотря на то, что средняя мощность хвостов спектра сигнала Рфт(ґ) в случае М=51 больше (рис. 6,6), чем в случае М=50 (рис. 5,6), их суммирование с высокочастотными гармониками сигнала Рдт(ґ) приводит к тому, что в случае М=51 высокочастотные гармоники оказываются заниженными (рис. 6,а) в сравнении со случаем М=50 (рис. 5,а). Для /¡=15000 Гц, Аґ=1/44100 с и М=50,51 номер нижней высокочастотной гармоники равен
пшт=18.
105
104
103
102
10
1
1 11 21 31 41 51 61 71 81 91
105
104
103
102
10
1
1 11 21 31 41 51 61 71 81 91
Рис. 6. Сравнение спектров дискретного, изображенного на рис. 2,а, и исходного континуального сигналов для интервала анализа длиной М=51 отсчет: а) амплитудный спектр дискретного сигнала; 6) амплитудный спектр исходного континуального сигнала, показаны только первые 100 коэффициентов
Трудоемкость вычисления Ри(М по формуле (26) составляет М(пшах - пшт + 1), ее можно оценить сверху величиной
а
2
а
0
6
010001003101313201310230530131010131013131010001533101310201
Тм = (1/2 - / Аґ)М2 + 2М,
отсюда трудоемкость нахождения наилучшего интервала анализа методом минимизации высокочастотных гармоник дискретного сигнала составит
Мшах
Т(0 = £ Тм = + 602 + с0 + й , (27)
где Q Mmax Mmin ;
a = z/3 ; b = zMmin + z /2 +1;
С = zMmm + (z + 2)Mmin + z / 6 + 1 ;
d = zMmin + 2Mmin; z = 1/2- fhAt.
Трудоемкость (27) похожа на трудоемкость (17) за исключением того, что в формуле (27) отсутствует множитель CN. На первый взгляд, трудоемкость (24) лучше, чем (27), но если учесть, что N>Q и N>Mmin, то трудоемкость (27) оказывается лучше трудоемкости (24). Кроме того, для вычисления Ph(M) не надо выполнять ДПФ для всего сигнала Рф(0.
В эксперименте скорость вычисления Ph(M) оказалась весьма высокой, что делает метод минимизации высокочастотных гармоник дискретного сигнала пригодным для практического применения.
Вместо минимизации функции Ph(M) можно минимизировать величину
Ph (M)
Rh (M) =
Pö (м)
СРАВНЕНИЕ ТОЧНОГО И ПРИБЛИЖЕННЫХ МЕТОДОВ ВЫБОРА ИНТЕРВАЛА АНАЛИЗА
На рис. 7, б - г приведены функции Р(М), Рь(М) и с(М), вычисленные для речевого сигнала, представляющего собой реализацию фонемы /и/ длительностью 298 мс. Вычисления проведены для т=4400 и М=71,...,491. В точке т=4400 значение сигнала равно нулю, и легко видно, что чем больше значение сигнала на правой границе интервала анализа отклоняется от нуля, тем больше зашумление спектра дискретного сигнала. В целом функции Р(М), Рь(М) и с(М) очень похожи, но наибольший интерес представляют точки минимумов этих функций.
- аналог относительного зашумления спектра, где Рд(М) есть средняя мощность дискретного сигнала Рдт( ) на интервале [ 0, 0+Т2)
[М/2] 2
п=-[М/2]
МЕТОД МИНИМИЗАЦИИ РАЗРЫВА ФУНКЦИИ НА КРАЯХ ИНТЕРВАЛА АНАЛИЗА
Для приближенного поиска интервалов анализа, на которых зашумление спектра дискретного сигнала минимально, очевидным является использование метода минимизации разрыва функции на краях интервала анализа. Минимизируемой величиной в этом случае будет функция с(М), даваемая выражением (15).
В примерах на рис. 2 и 3 функция с(М) достигает минимума в тех же точках, что и функция Р (М), т.е. в этих случаях решение, полученное с помощью приближенного метода минимизации разрыва функции на краях интервала анализа, совпадает с точным решением.
Трудоемкость приближенного поиска наилучшего интервала анализа методом минимизации разрыва функции на краях интервала является линейной и равна
Т ©) = Q .
Рис. 7. Сравнение функций Р(М), Р*(М) и с(М) на примере реального речевого сигнала: а) фрагмент речевого сигнала, номер левого граничного отсчета равен 4470, ему соответствует номер 1 на оси абсцисс, длина фрагмента - 421 отсчет, N=13149; б) функция Р(М); в) функция Р*(М), /¡,= 15000 Гц, Д¿=1/44100 с; г) функция с(М). Для всех трех функций левая граница интервалов анализа равна т=4400, длина интервалов меняется от Мт;п=71 до Мтах=491, по оси абсцисс отложены номера интервалов
Первые десять минимумов с наименьшими значениями функции Р(М) приведены в табл. 1. Точки первых двух минимумов очень близки друг к другу, так что в практических приложениях их можно объединить в одну, поскольку за время двух отсчетов речевой сигнал не может существенно измениться.
Т аблица 1
Длина интервала, М Номер интервала Р(М)
173 103 73,4
171 101 134,0
257 187 167,8
449 379 176,5
460 390 266,2
357 287 478,0
349 279 506,8
91 21 920,3
432 362 1183,1
345 275 1444,2
Первые десять минимумов с наименьшими значениями функции Рь(М) приведены в табл. 2. Здесь точки первых двух минимумов также очень близки друг к другу, и их можно объединить в одну. В первую десятку минимумов функции Рь(М) вошли две точки М=103 и М=161, для которых зашумление спектра достаточно велико (Р^103)=2704,3; Р^161)=2697,5), что объясняется приближенным характером метода минимизации высокочастотных гармоник дискретного сигнала. Точки оставшихся восьми минимумов функции Рь(М) очень близки к соответствующим точкам минимумов функции Р (М) из табл. 1, максимальное расхождение между точками минимумов составляет 3 отсчета. Точное совпадение минимумов функций Р (М и Рн(М) происходит в единственной точке М=357. Длина интервала М=357 отсчетов замечательна еще тем, что она равна периоду основного тона данного фрагмента речевого сигнала. Также необходимо обратить внимание на то, что точки минимумов из табл. 2 расположены в совершенно ином порядке, чем точки из табл. 1.
Первые десять минимумов с наименьшими значениями функции с(М) приведены в табл. 3. Как и в случае функции Ри(М), в первую десятку минимумов функции с(М) вошли точки М=103 и М=161, для которых зашумление спектра достаточно велико. Точки оставшихся восьми минимумов функции с(М) полно-
стью совпадают с соответствующими точками минимумов функции Р(М) из табл. 1, т.е. в данном случае метод минимизации разрыва функции на краях интервала анализа показывает лучший результат, чем метод минимизации высокочастотных гармоник дискретного сигнала. Этот факт вместе с линейной трудоемкостью метода минимизации разрыва функции делает его предпочтительным в сравнении с методом минимизации высокочастотных гармоник. Порядок расположения точек в табл. 3 отличается от того же порядка в табл. 1. Это говорит о том, что по величине с(М) разрыва функции на краях интервала анализа можно лишь очень грубо судить о величине Р(М) зашумле-ния спектра дискретного сигнала.
Таблица 2
Длина интервала, М Номер интервала Рн(М)
459 389 35,3
461 391 90,0
260 190 117,8
450 380 123,8
172 102 129,3
103 33 152,6
161 91 252,7
93 23 475,9
350 280 491,8
357 287 580,9
Таблица 3
Длина интервала, М Номер интервала с(М)
173 103 3
257 187 47
349 279 164
449 379 347
91 21 555
460 390 995
357 287 1274
103 33 1588
161 91 2013
345 275 2172
ЛИТЕРАТУРА
1. Фант Г. Акустическая теория речеобразования: Пер. с англ. М.: Наука, 1964. 284 с.
2. Кодзасов С.В., Кривнова О.Ф. Общая фонетика. М.: Рос. гос. гуманит. ун-т, 2001. 592 с.
3. Голд Б., Рэйдер Ч. Цифровая обработка сигналов: Пер. с англ. М.: Сов. радио, 1973. 368 с.
4. Гоноровский И.С. Радиотехнические цепи и сигналы. М.: Радио и связь, 1986. 512 с.
5. ОппенгеймА.В., ШаферР.В. Цифровая обработка сигналов: Пер. с англ. М.: Связь, 1979. 416 с.
6. Рабинер Л.Р. Шафер Р.В. Цифровая обработка речевых сигналов: Пер. с англ. М.: Радио и связь, 1981. 496 с.
7. Рабинер Л. Гоулд Б. Теория и применение цифровой обработки сигналов: Пер. с англ. М.: Мир, 1978. 848 с.
8. Толстов Г.П. Ряды Фурье. М.: Наука, 1980. 384 с.
Статья представлена кафедрой прикладной информатики факультета информатики Томского государственного университета, поступила в научную редакцию 23 июля 2003 г.