Книга. Речевых сигналов
 Скачать 1.72 Mb.
|
|
2.1.6. Быстрое преобразование Фурье Пару ДПФ часто записывают в виде 1 , 0, , 1, 0 1 1 , 0, , 1, 0 N kn A X W k N k n n N kn X A W n N n k N k − ⎧ = = − ∑ ⎪ ⎪ = ⎨ − ⎪ − = = − ∑ ⎪ = ⎩ … … (2.25) где ( ) exp 2 , ( ) W j N X rep X n t n N = − = Δ π – отсчеты сигнала, ( ) A rep F k f t k N = Δ Δ – коэффициенты ДПФ. Для вычисления одного элемента последовательности , 0, , 1, k A k N = − … необходимо примерно 2 N операций комплексных ум- ножений и сложений. Число операций возрастает пропорционально квад- рату размерности ДПФ [43]. Однако если N не является простым числом и может быть разложено на множители, процесс вычислений можно уско- рить, разделив анализируемый набор отсчетов на части, вычислив их ДПФ и объединив результаты. Такие способы вычисления ДПФ называются бы- стрым преобразованием Фурье (БПФ; Fast Fourier Transform, FFT). Пре- 10 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 0 8 X r ( ) r а) 0 0 r = 10 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 0 4.09 X r ( ) r б) 1 0 r = 10 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 0 4 X r ( ) r 42 имущества алгоритма БПФ особенно проявляются с ростом N , что суще- ственно при обработке массивов большой размерности [38]. Существует несколько разновидностей алгоритма БПФ. Изложим две модификации: с прореживанием по времени и по частоте. Прореживание по времени Разделим последовательность Xn , состоящую из N отсчетов, на две подпоследовательности Yn и Zn, каждая из которых включает 2 N отсче- тов (рис. 2.5). Рис. 2.5. График последовательностей Xn , Yn и Zn Отсчеты Yn образованы из четных отсчетов исходной последова- тельности Xn , а отсчеты Zn – из нечетных: 2 2 1 , , 0,1, 2, , 1 2 n n N Y X Z X n n n + = = = − … . (2.26) Поскольку подпоследовательности Yn и Zn состоят из 2 N отсчетов каждая, ДПФ для них имеет вид 2 2 ( 2) 1 , 0 0,1, , 1 ( 2) 1 , 0 k n k n N B Y W k n n k N N C Z W k n n − ⎫ = ⎪ ∑ ⎪ = = − ⎬ − ⎪ = ∑ ⎪ = ⎭ … . (2.27) 43 Необходимо найти последовательность Ak , которую можно предста- вить через четные и нечетные элементы исходной последовательности Xn: ( 2) 1 (2 1) 2 ( ) 0 , 0 1. 2 N k n kn A Y W Z W k n n n N k B W C при k k k − + = + = ∑ = = + ≤ ≤ − (2.28) Поскольку Bk и Ck периодичны с периодом 2 N , можем записать 2 п р и 0 1. 2 2 N k N k A B W C B W C k k N k k k k + = + = − ≤ ≤ − + (2.29) Таким образом, первые 2 N и последние 2 N отсчетов ДПФ от Xn могут быть получены комбинацией отсчетов ДПФ двух подпоследователь- ностей Yn и Zn. На рис. 2.6 показан сигнальный граф, наглядно представ- ляющий процедуру конструирования отсчетов Ak из отсчетов Bk и Ck для случая 8 N = . Рис. 2.6. Сигнальный граф, представляющий процедуру конструирования отсче- тов Ak из отсчетов Bk и Ck для случая 8 N = 6 3 x y = ДПФ (N=4) y x2 1 = y x4 2 = ДПФ (N=4) x z0 1 = x z1 3 = x z2 5 = y x7 3 = B0 B1 B2 B3 C0 C1 C2 C3 A0 A1 A2 A3 A4 A5 A6 A7 0 W 1 W 2 W 3 W 4 W 5 W 6 W 7 W y x2 1 = 44 Поскольку удалось задачу вычисления N-точечного ДПФ редуциро- вать к задаче вычисления двух 2 N -точечных ДПФ, естественно попы- таться «развить» успех в данном направлении. На рис. 2.7 показаны два следующих аналогичных шага, после которых отсчеты сигнала Xn оказы- ваются связанными с коэффициентами ДПФ Ak своеобразными нитями- операциями, похожими на «бабочку». Итак, для случая 3 2 8 N = = вычисления совершаются в три этапа. На первых двух этапах находят некие «промежуточные» массивы из вось- ми точек каждый. На третьем этапе вычисляется «окончательный» восьми- точечный массив. Для определения каждого элемента из этих трех масси- вов необходимо выполнить одно комплексное умножение и одно ком- плексное сложение. Итого получается 3 8 24 ⋅ = комплексных умножений и сложений вместо 8 8 64 ⋅ = при «лобовых» вычислениях. Обобщая рассуждения на случай N -точечных массивов, заключаем, что для вычислений в соответствии с алгоритмом БПФ необходимо log2 N N комплексных умножений и сложений, тогда как при прямых вы- числениях требуется 2 N операций. Прореживание по частоте Как и прежде, разделим последовательность Xn на две подпосле- довательности Yn и Zn, каждая из которых включает 2 N отсчетов. Од- нако теперь Yn будет состоять из первых 2 N отсчетов, а Zn – из по- следних 2 N отсчетов , , 0, 1, 2, , 1 2 2 N Y X Z X n n n n n N … = = = − + . (2.30) Тогда ( ) ( 2) 1 2 0 N k n N k n A Y W Z W k n n n ⎛ ⎞ + − ⎜ ⎟ = + = ∑ ⎜ ⎟ = ⎜ ⎟ ⎝ ⎠ ( 2) 1 2 , 0 1 0 N k N k n Y Z W W k N n n n ⎛ ⎞ − ⎜ ⎟ = + ≤ ≤ − ∑ ⎜ ⎟ = ⎜ ⎟ ⎝ ⎠ . (2.31) 45 Рис. 2.7. Связь отсчетов сигнала Xn с коэффициентами ДПФ Ak X0 X4 X2 X6 X1 X5 X3 X7 ДПФ ( N=2) ДПФ ( N=2) ДПФ ( N=2) ДПФ ( N=2) B0 B1 B2 B3 0 W 2 W 4 W 6 W A0 A1 A2 A3 A4 A5 A6 A7 0 W 1 W 2 W 3 W 4 W 5 W 6 W 7 W C0 C1 C2 C3 0 W 2 W 4 W 6 W X0 X4 X2 X6 X1 X5 X3 X7 0 W 4 W 0 W 4 W 0 W 4 W 0 W 4 W 0 W 2 W 4 W 6 W 0 W 2 W 4 W 6 W 0 W A0 1 W 2 W 3 W 4 W 5 W 6 W 7 W A1 A2 A3 A4 A5 A6 A7 46 Теперь рассмотрим четные и нечетные отсчеты массива Ak , т.е. осуществим прореживание по частоте: , , 0, 1, 2, , 1 2 2 1 2 N R A S A k k k k k = = = − + … . (2.32) Тогда для четных отсчетов ( ) ( 2) 1 2 , 0 1 2 2 0 N N k n R A Y Z W k k k n n n − = = + ≤ ≤ − ∑ = . (2.33) Соотношение (2.33) есть 2 N -точечное ДПФ массива Y Z n n + , пред- ставляющего собой сумму первых 2 N и последних 2 N отсчетов исход- ного временного массива Xn . Аналогично для нечетных спектральных от- счетов ( ) ( 2) 1 2 , 0 1 . 2 1 2 0 N N n k n S A Y Z W W k k k n n n − ⎡ ⎤ = = − ≤ ≤ − ∑ + ⎢ ⎥ ⎣ ⎦ = (2.34) Соотношение (2.34) есть 2 N -точечное ДПФ массива ( ) n Y Z W n n − , представляющего собой взвешенную разницу первых 2 N и последних 2 N отсчетов исходного временного массива Xn . Итак, задача вычисления N -точечного ДПФ снова была редуцирова- на, хотя и несколько иным способом, к задаче вычисления двух 2 N - точечных ДПФ. Сигнальный граф для этого случая показан на рис. 2.8. Таким образом, при прореживании по частоте и по времени проце- дуру вычислений делят на log2 N этапов. При этом на каждом этапе на определение элементов N-точечного массива затрачивается N комплекс- ных сложений и умножений. В результате вычисления производятся при- мерно за log2 N N комплексных сложений и умножений против 2 N – при «лобовых» вычислениях. 47 Рис. 2.8. Сигнальный граф редуцирования задачи вычисления N -точечного ДПФ к задаче вычисления двух 2 N -точечных ДПФ 2.1.7. Вычисление спектрограммы Спектрограммой сигнала называется его мгновенный спектр, зави- сящий от времени: ( , ) ( ) t j t F t X t e d t t T ω ω − = ∫ − . (2.35) Для вычисления спектрограммы дискретного сигнала его разбивают на сегменты (возможно, с перекрытием). Для каждого сегмента находят его спектр в виде коэффициентов ДПФ. Набор спектров и образует спек- трограмму (рис. 2.9) [38]. Разрешающая способность по частоте такого спектрального анализа определяется величиной 1 1 f T Δ = , а разрешающая способность по време- ни – величиной 1 T , если сегменты анализируемого процесса не перекры- ваются. Если же сегменты перекрываются, то разрешение по времени мо- жет быть равным даже 1 1 t T N Δ = , где 1 N – число отсчетов сегмента, под- вергаемого преобразованию Фурье (поэтому число 1 N часто называют па- раметром БПФ и принимают равным степени «2»). Однако на практике X0 X1 X2 X3 X4 X5 X6 X7 0 W 0 W − 1 W 1 W − 2 W 2 W − 3 W 3 W − A0 A1 A2 A3 A4 A5 A6 A7 ДПФ ( N=4) ДПФ ( N=4) 48 степень перекрытия выбирают из неких «разумных» предпосылок, по- скольку при чересчур высокой степени перекрытия объем вычислений мо- жет стать неприемлемо высоким. Рис. 2.9. Спектрограмма сигнала ( ) X t : A – амплитуда, f – частота, t – время Один и тот же термин «спектрограмма» применяют как к комплекс- ной функции частоты и времени, так и к ее модулю (набору амплитудных спектров). В программе MATLAB для получения комплексного массива B ис- пользуют функцию specgram c синтаксисом B = specgram(x, Nfft, Fs, window, numoverlap), где x – массив отсчетов исходного сигнала; Nfft – параметр ДПФ, вы- числяемого с помощью алгоритма БПФ; Fs – частоты дискретизации сиг- нала; window – окно для взвешивания сегментов сигнала; numoverlap – ко- личество перекрывающихся отсчетов сегментов. Таким образом вычисля- ется и строится модуль массива B, т.е. набор амплитудных спектров, уро- вень значений которых кодируется цветом. Рассмотрим пример вычислений спектрограммы из командного окна. После указания пути к папке по имени \toolbox\signal\signal, где на- ходятся необходимые для работы исходные данные и программы, откройте файл mtlb.mat, в рабочем пространстве при этом появится информация о считывании в среду MATLAB двух переменных – одномерного вещест- венного массива по имени mtlb из 4001 отсчета и числа Fs = 7418 Гц. A f 1 T 1 T 1 T 1 T t T t 49 Сигнал mtlb (слово “MATLAB”) можно прослушать с помощью ко- манды wavplay(mtlb, Fs, 'async'). Можно построить график сигнала mtlb (рис. 2.10) с помощью следующих команд: >> t=1:4001; >> plot(t,mtlb). Рис. 2.10. График речевого сигнала: слово «MATLAB» Теперь сформируем окно Бартлетта протяженностью 512 отсчетов: >> win = bartlett (512), график которого можно просмотреть с помощью команд >> x=1:512; plot(x,win). Наконец, выполним команду >> specgram(mtlb, 512, Fs, win, 256). При этом будем наблюдать спектрограмму, состоящую (при данных значениях параметров) из 14 спектров сегментов (рис. 2.11). Нетрудно подсчитать, что длительность анализируемого сигнала со- ставляет 4001/ 7418 0,54 T = = с, диапазон анализируемых частот равен [0, Fs/2], т.е. [0, 3709] Гц. Для Nfft=512 получаем, что длительность сегмента равна 512 / 7418 0,069 1 T = = − с, разрешающая способность по частоте 7418/ 512 14,5 1 f F N s Δ = = = Гц. Для степени перекрытия 256 отсчетов разрешение по времени составляет 256 / 7418 0,035 = с. 50 Рис. 2.11. Спектрограмма сигнала в окне Бартлетта 2.2. Корреляционная обработка сигналов. Процедура и примеры применения Задача выявления периодических колебаний на фоне шумов встреча- ется практически повсеместно. Периодически меняется интенсивность из- лучения такого космического объекта, как Солнце. Планетные системы, вращаясь вокруг единого центра, периодически заслоняют друг друга, что позволяет астрономам обнаруживать объекты и измерять их характеристи- ки. Наконец, все земные явления и процессы несут на себе «печать» пе- риодичности вследствие вращения Земли вокруг Солнца и собственной оси, а также вращения Луны вокруг Земли. Выявление периодичности присутствует в задачах прогнозирования погоды и климата, тенденций развития растительного и животного мира, социальной и экономической активности людей [38, 43]. Разумеется, периодичность свойственна не только явлениям косми- ческого и планетарного масштаба. Анализ шумов и вибраций механизмов позволяет осуществлять раннее выявление неисправностей при техниче- ской диагностике. Анализ шумов сердца и легких – хорошо известные примеры медицинской диагностики. А задачи шумопеленгования либо ак- тивной локации объектов являются базовыми при обнаружении, слежении и классификации объектов-целей в военном деле. Зачастую периодический характер явлений и процессов замаскирован шумами, порой весьма интенсивными. В этой связи становится понятной ак- туальность задачи выявления «скрытых периодичностей». Одним из эффек- тивных методов решения такой задачи является корреляционный анализ. 51 Для аддитивной смеси ( ) ( ) ( ) Y t S t t ξ = + (2.36) сигнала 0 ( ) cos(2 ) S t A f t π ϕ = + и шума ( ) t ξ отношение сигнал-шум (от- ношение средней мощности 2 2 A сигнала к средней мощности D ξ шума) 2 2 вх A D = ρ ξ может быть разным: как очень малым 1 вх ρ << , так и большим. На рис. 2.12 показаны графики аддитивной смеси для 10 lg 17 дБ, вх дБ вх = ⋅ = − ρ ρ построенные для соотношения (2.36), т.е. «сигнал + шум», где t i t = Δ , j t = Δ τ , B – верхняя граница частоты белого шума на интервале [ ] 0, B Гц. Как видим, на фоне интенсивного шума сиг- нал не наблюдается. % Вычисление амплитуды A сигнала ( ) Y t >> Dksi = 1; % дисперсия шума >> Rvh = -17; >> A = 10^((10*log10(2)+Rvh)/20);% Амплитуда сигнала >> %== сигнал плюс шум (3 периода гарм.сигнала)==== >> B = 5000; % верхн.границна частоты шума >> fi = rand(1)*2*pi; % фаза >> f0 = 220; >> f0B = f0/B; % входная частота сигнала >> Ngraf = ceil(3/f0B)*2+1; % (Ngraf = 139) >> i = 1: Ngraf; >> S = A*cos(pi*f0B*i+fi); % функция вх. сигнала >> ksi = randn(1, Ngraf); % функция шума >> Sksi = S + ksi; % сигнал + шум >> subplot(3,1,1); plot(i,S); % график сигнала >> title('Signal'); >> subplot(3,1,2); plot(i,ksi); % график шумa >> title('Noise'); >> subplot(3,1,3); plot(i,Sksi,'r'); % график "сигнал + шум" >> title('Signal+Noise'); Покажем, что корреляционный анализ случайного процесса помогает решить задачу выявления периодического сигнала на фоне шума. Поскольку составные части процесса ( ) Y t статистически независимы, ( ) ( ) ( ) K K K Y S τ τ τ ξ = + , (2.37) 52 0 2 s i n 2 ( ) c o s ; ( ) 2 2 A B K K D S B π τ τ ω τ τ ξ ξ π τ = = , (2.38) где B – верхняя граничная частота шума ( ) t ξ Рис. 2.12. Графики аддитивной смеси Для построения графиков корреляционных функций (2.37) – (2.38) необходимо дискретизовать функции (2.38) с шагом 1 2 t B Δ = , в результа- те чего получим 2 s i n 0 ( ) c o s ; ( ) 2 f A i K i K D S B i π π τ τ ξ ξ π ⎛ ⎞ = = ⎜ ⎟ ⎝ ⎠ , где 100 200, 5000 0 f B = = … Как следует из соотношений (2.37) – (2.38), форма корреляционной функции процесса ( ) Y t для 1 2B τ ≥ практически не отличается от формы гармонического сигнала ( ) S t . Этот факт позволяет по частоте нуль- пересечений оценить частоту 0 f . Понятно, что действенность такого спо- соба оценивания 0 f особенно эффективна при малых отношениях сигнал- шум вх ρ 0 20 40 60 80 100 120 140 -0.2 0 0.2 Signal 0 20 40 60 80 100 120 140 -5 0 5 Noise 0 20 40 60 80 100 120 140 -5 0 5 Signal+Noise 53 2.2.1. Выявление периодического колебания на фоне шума На практике можно лишь оценить функцию корреляции, поэтому, естественно, результаты выявления гармонического сигнала на фоне шума будут не такими «красивыми». Структурная схема коррелометра показана на рис. 2.13 [38, 43]. Рис. 2.13. Структурная схема коррелометра Для отрезка ( ) Y t T процесса ( ) Y t несмещенная оценка функции кор- реляции имеет вид 1 * ( ) ( ) ( ) 0 T K Y t Y t dt Y T T T τ τ τ τ − = + ∫ − . (2.39) Покажем, что коррелометр действительно способен повышать отно- шение сигнал-шум: тем больше, чем больше некоррелированных отсчетов шума N содержится в реализации анализируемого процесса. Для аддитивной смеси сигнала и шума (2.36) имеем ( ) ( ) ( ) ( ) ( ) K K K K K Y S S S τ τ τ τ τ ξ ξ ξ ∗ ∗ ∗ ∗ ∗ = + + + . (2.40) Для статистически независимых сигнала и шума математическое ожидание оценки (2.40) будет иметь вид * ( ) ( ) ( ) ( ) K K K K Y Y S τ τ τ τ ξ = = + , (2.41) где ∗ – символ математического ожидания. При T к τ τ < << имеем 2 0 * ( ) ( ) ( ) c o s 2 A K K K Y Y S τ τ τ ω τ = ≈ = . (2.42) Рассмотрим теперь дисперсию оценки выражения (2.40): * * * * [ ( ) ] [ ( ) ] [ ( ) ] [ ( ) ] * [ ( ) ] 2 2 * * * * D K D K D K D K Y S S D K K K S K K K K S S S τ τ τ τ ξ ξ τ ξ ξ ξ = + + + + + + + Y(t) ) t ( Y ) t ( Y τ + ) t ( Y τ + ∫ − − x T 0 T 1 τ 1 * ( ) ( ) ( ) 0 T K Y t Y t dt Y T T T τ τ τ τ − = + ∫ − 54 2 2 2 * * * * * * 2 * * K K K K K K K K K S S S S K K K S S + + + + + ξ ξ ξ ξ ξ ξ ξ (2.43) После громоздких, хотя и несложных, выкладок получаем ( ) 0 0 0 0 4 * [ ( ) ] 1 c o s 2 4 0 0 2 2 2 c o s 2 2 2 к . к в к T A D K d Y T T D D A T T ⎛ ⎞ ≈ − + ∫ ⎜ ⎟ ⎝ ⎠ + + + τ τ ω τ τ ξ ξ ω τ τ τ (2.44) Ненулевой вклад обеспечивают слагаемые ( ) , ( ) D K D KS τ τ ξ ∗ ∗ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ и K K K S S ξ ξ ⎡ ⎤ ⎢ ⎥ ∗ ∗ ⎢ ⎥ ⎣ ⎦ Подставляя 2 в х 2 A D ξ ρ в (2.43), имеем 2 2 * [ ( ) ] 1 c o s 2 0 0 0 T D в х D K d Y T T ρ τ ξ τ ω τ τ ⎛ ⎞ ≈ − ∫ ⎜ ⎟ ⎝ ⎠ + ( ) 2 4 1 4 2 c o s 2 0 0 2 2 2 к . к в к в х D A в х T T + ⋅ + ⋅ + ξ ρ ω τ τ τ ρ . (2.45) Нетрудно увидеть, что для малых значений входного отношения сиг- нал-шум, наиболее интересных для исследования ( 1 вх ρ << ), наибольший вклад в дисперсию оценки функции корреляции делает слагаемое ( ) D K τ ξ ∗ ⎡ ⎤ ⎢ ⎥ ⎣ ⎦ , т.е. 2 4 2 1 4 * * [ ( ) ] [ ( ) ] к . к в 2 2 к.кв D A D K D K Y T T в х ξ τ τ τ ξ τ ρ ≈ ≈ = ⋅ . (2.46) 55 Таким образом, если трактовать выражение (2.42) как сигнал на вы- ходе коррелометра, а выражение (2.46) – как мощность шума на выходе коррелометра, тогда отношение сигнал-шум на выходе коррелометра будет иметь вид 4 4 ( 2 ) 8 к . к в в ы х * 2 [ ( ) ] 8 ( 2 ) 2 2 к . к в в х в х 2 A T A D K D T T B τ ρ τ ξ ξ τ ρ ρ = = = = = (2.47) Обозначая 2 к кв T N τ = , получаем 2 вых вх 2 N = ρ ρ (2.48а) или, если отношение сигнал-шум измеряется в децибелах, 2 1 0 l g в ы х . д Б в х 2 2 1 0 l g 1 0 l g 2 . в х . д Б N N ⎛ ⎞ = ⋅ = ⎜ ⎟ ⎝ ⎠ = + ⋅ − ⋅ ρ ρ ρ (2.48б) Соотношения (2.48а), (2.48б) показывают, что применение корреля- тора позволяет существенно повысить отношение сигнал-шум за счет ус- реднения некоррелированных отсчетов шума. Используя соотношения (2.48а) и (2.48б), можно оценить объем N экспериментальной выборки отсчетов смеси сигнала и помехи, который необходим для обеспечения заданного отношения сигнал-шум на выходе коррелятора: 2 вых 2 вх P N ρ = , (2.49а) 0,1( 2 10 lg 2) вых.дБ вх.дБ 10 N − + ⋅ = ρ ρ . (2.49б) Формулы (2.49а), (2,49б) чрезвычайно полезны при планировании эксперимента. Они позволяют определить объем выборки отсчетов про- цесса ( ) Y t при заданных входном и выходном отношениях сигнал-шум. Вычислим N, которое необходимо для обеспечения выходного отно- шения сигнал-шум, равного 10 дБ, и проведем моделирование соответст- вующего процесса на выходе коррелятора. % вычисление N , которое обеспечивает Rvyh=10 дБ % Rvh – отношение сигнал-шум на входе коррелометра (Rvh =-17 - 10). Rvh = - 17; Rvyh=10; N = 10^(0.1*(Rvyh-2*Rvh+10*log10(2))); 56 В результате вычислений получаем N = 50238. Наконец, смоделиру- ем процедуру обработки сигнала ( ) Y t T коррелометром. На рис. 2.14 приведены графики оценки функции корреляции (т.е. графики сигнала на выходе коррелометра). Рис. 2.14. Графики сигнала на выходе коррелометра % моделирование коррел. обработки смеси сигнала с шумом % =ceil(N); >>figure; i = 1:N; fi = rand(1)*2*pi; >> f0 = 200; B = 5000; >> f0B = f0/B; >> A = 10^((10*log10(2)+Rvh)/20); >> Ngraf = ceil(3/f0B)*2+1; >>S = A*cos(pi*f0B*i+fi); % N отсчетов сигнала >>ksi = randn(1,N); % N отсчетов шума >>Sksi = S + ksi; % N отсчетов смеси >>[K_y,tau] = xcorr(Sksi,Ngraf,'unbiased'); % Ngraf*2+1 значений функц.коррел. >>subplot(2,1,1); plot(tau,K_y); % график оценки функции коррел. >> title('otsenka corr sig+nois'); >>subplot(2,1,2); plot(tau(Ngraf+2: 2*Ngraf+1),K_y(Ngraf+2: 2*Ngraf+1)); % фрагмент графика >> title('frogment otsenka corr sig+nois') 57 2.2.2. Оценивание длины периода основного тона речевого сигнала. Экспериментальное измерение в системе MATLAB Высота тона – основное различие между мужскими и женскими спи- керами. Высота тона человека определяется в голосовых связках, и темп, с которым вибрируют вокальные складки, представляет собой частоту высо- ты тона. Когда воздух проходит через вибрирующие вокальные складки, создаются также гармоники. Оценивание периода (или частоты) основного тона – одна из наибо- лее важных задач в обработке речевых сигналов [39]. Выделители основ- ного тона (ОТ) используются в вокодерах, системах распознавания и ве- рификации дикторов, в устройствах, предназначенных для глухих. Для решения задачи предложен ряд способов. Например, для оценивания дли- ны периода основного тона речевого сигнала можно использовать кратко- временную автокорреляционную функцию [34]. Для ее вычисления разра- ботаны эффективные алгоритмы. При длине окна N > 256 отсчетов лучше всего использовать БПФ, которое предполагает выполнение 8(2·log2N+3)N умножений. С помощью системы MATLAB можно провести экспериментальное измерение частоты основного тона речевого сигнала [38], которое базиру- ется на следующей априорной информации: – частота основного тона голоса человека находится в пределах 60 – 200 Гц; – для качественной передачи речевого сигнала достаточно использо- вать полосу частот 20 Гц – 5 кГц. Вначале проводят проверку работоспособности системы ввода рече- вого сигнала в компьютер. Для этого используют телефонную гарнитуру и стандартную программу записи звука, активизируемую в Windows с по- мощью меню Пуск > Программы > Стандартные > Развлечения > Звукоза- пись. Далее, используя MATLAB, вводим речевой сигнал в компьютер (на- пример слово «барбарис»): >>Y = wavrecord(15000); % частота дискретизации 11025 Гц. Проконтролировать на слух результаты записи можно, выполнив следующую команду: >>wavplay(Y). Далее нужно выделить из речевого сигнала фрагмент с гласным зву- ком. С этой целью активизируем программу sptool. Анализируем записанный сигнал и находим границы самого протя- женного гласного звука (рис. 2.15). Записываем координаты маркеров, со- ответствующих этим границам: 1 t = 0,816 и 2 t = 1,018. 58 Рис. 2.15. График речевого сигнала (слово «барбарис») Формируем массив чисел 1 Y из выделенного фрагмента: >>Fs = 11025; t1 = 0.816; t2 = 1.018; % границы фрагмента >>j1 =ceil(t1*Fs); j2 = ceil(t2*Fs); % номера отсчетов границ фрагме- нта >>Y1 = Y(j1:j2); % массив отсчетов фрагмента. Импортируя массив 1 Y в рабочее пространство программы sptool, можно визуально и на слух проконтролировать полученный фрагмент (рис. 2.16). Рис. 2.16. График сформированного фрагмента (гласный звук) 59 Наконец, выполняем автокорреляционный анализ выделенного фрагмента и оцениваем частоту основного тона: >> [Kzvuk,lags] = xcorr(Y1,400). Результат вычисления корреляционной функции вновь целесообраз- но импортировать в среду sptool, где с помощью вертикальных маркеров (рис. 2.17) провести измерения периода колебаний. Рис. 2.17. График корреляционного анализа гласного звука Частота основного тона определяется по формуле 0 0 1 f T = . Изме- рить период 0 T очень просто: в программе sptool установите маркеры примерно так, как показано на рис. 2.17. Между маркерами в данном слу- чае уложилось восемь периодов 0 T . Величина периода вычисляется следующим образом: ( ) ( ) 0 2 1 8 0,06549 0,00698 8 0,007314 T t t = − = − = В результате таких измерений получаем оценку частоты основного тона 0 137 f = Гц. Теперь рассмотрим проблемы, обусловленные наличием нескольких гармонических компонентов. Анализ рис. 2.16 показывает, что измерение частоты основного тона произвести весьма непросто несмотря на то, что шумовой сигнал действи- тельно практически полностью подавлен. Дело в том, что помимо основ- ного тона в гласном звуке присутствуют обертоны – гармоники с кратны- ми частотами. Более того, как показывает практика измерений, мощность 60 обертонов может быть выше мощности основного тона, поэтому, измеряя частоту основного тона, легко ошибиться, приняв период колебания обер- тона за период основного тона. Смоделируем рассмотренную ситуацию. Предположим, что сигнал представляет собой сумму двух гармоник одинаковой амплитуды, но частота второй гармоники вдвое выше частоты первой: 1 2 1 0 1 2 0 2 ( ) ( ) ( ) ; ( ) c o s (2 ) ; ( ) c o s (4 ) . S t S t S t S t A f t S t A f t = + = + = + π ϕ π ϕ (2.50) В силу статистической независимости случайных величин 1 ϕ и 2 ϕ корреляционная функция такого сигнала равна сумме корреляционных функций слагаемых: ( ) 1 2 2 0 0 ( ) ( ) ( ) c o s c o s 2 2 S S S A K K K τ τ τ ω τ ω τ = + = + . (2.51) Обобщая рассмотренную модель на случай трех гармоник (частота третьей гармоники в три раза больше частоты первой), получим ( ) 1 2 3 0 0 0 ( ) ( ) ( ) ( ) 2 c o s c o s 2 c o s 3 2 S K K K K S S S A = + + = = + + τ τ τ τ ω τ ω τ ω τ (2.52) Построим графики функций корреляции (2.51) и (2.52) (с точностью до множителя 2 2 A ) (рис. 2.18). % функция коррел. полигармонич. сигнала >>figure; >> f0 = 137; Fs = 11025; itau = 0:400; >>Ks2 = cos(2*pi*f0*itau/(2*Fs))+ cos(4*pi*f0*itau/(2*Fs)); >>Ks3=cos(2*pi*f0*itau/(2*Fs))+cos(4*pi*f0*itau/(2*Fs))+… cos(6*pi*f0*itau/(2*Fs)); >>subplot(2,1,1); plot(itau,Ks2); >>title('dve gormonic'); grid on; >>subplot(2,1,2); plot(itau,Ks3); >>title('tri gormonic');>>grid on; Период основного тона на рис. 2.18 виден неплохо: в моменты вре- мени, кратные этому периоду, все гармонические компоненты функции корреляции суммируются синфазно, в результате чего обеспечивается «дружный» периодический «всплеск». Вместе с тем очевидно, что вклад обертонов настолько существен, что назвать ситуацию удобной для прак- тического анализа нельзя. 61 Рис. 2.18. Графики функций корреляции Особенно сложная ситуация, когда мощность обертонов выше мощ- ности основного тона. Действительно, для простой двухкомпонентной модели (2.52) в предположении, что мощность обертона в два раза (т.е. на 6 дБ) выше мощности основного тона ( ) 1 2 2 0 0 ( ) ( ) ( ) cos 2 cos 2 2 S S S A K K K τ τ τ ω τ ω τ = + = + ⋅ . (2.53) Результат показан на рис. 2.19. Выполним моделирование. % функц. коррел., обертон мощнее в 2 раза >>figure; >>Ks2 = cos(2*pi*f0*itau/(2*Fs))+ 2*cos(4*pi*f0*itau/(2*Fs)); >>plot(itau,Ks2); % две гармоники >>title(‘Две гармоники, обертон мощнее в 2 раза’); >>grid on. Рассмотренный пример наглядно демонстрирует мешающий харак- тер мощного обертона. Сравнивая рис. 2.17 и 2.18, видим, что результаты 62 моделирования объясняют природу отмеченных ранее трудностей измере- ния периода (частоты) основного тона. Рис. 2.19. График функции корреляции (обертон мощнее в два раза) 2.2.3. Преобразование Фурье функции корреляции как способ выявления периодического колебания Решение указанной выше проблемы можно найти путем вычисления преобразования Фурье от корреляционной функции. Гармоники, из кото- рых состоит функция корреляции, превратятся в спектральные пики, раз- несенные по частоте. Тем самым будет решена проблема разделения ос- новного тона и обертонов. В сущности, мы вплотную подошли к идее вычисления спектра мощности стационарного случайного процесса (ССП), провозглашенной в свое время независимо друг от друга советским математиком А. Хинчи- ным и американским кибернетиком Н. Винером. Пара преобразований Ви- нера – Хинчина – это преобразования Фурье, связывающие между собой функцию корреляции и спектр мощности ССП [38]: ( ) ( )exp( 2 ) , ( ) ( )exp( 2 ) . P f K j f d K P f j f df ∞ −∞ ∞ −∞ ⎧ = − ⎪ ⎪ ⎨ ⎪ = ⎪ ⎩ ∫ ∫ τ π τ τ τ π τ (2.54) Применяя преобразование Фурье к функции (2.53), получим 63 2 0 0 0 0 ( ) [ ( ) ( ) ] 4 2 [ ( 2 ) ( 2 ) ] . 2 A P f f f f f A f f f f = + + − + + + + − δ δ δ δ (2.55) График функции (2.55) показан на рис. 2.20. В области положитель- ных частот мы имеем два идеально разделяемых спектральных пика. Разумеется, результат изме- рений спектра мощности, именуе- мый «оценкой спектра мощности», будет выглядеть несколько хуже, как и «положено» всякой оценке. Оценка спектра будет отличаться от истинного спектра на величину ошибки измерений, содержащей систематическую и случайную со- ставляющие. Оценку спектра мощности с учетом соотношений (2.54) естест- венно сформировать в виде ( ) ( )exp( 2 ) T T P f K j f d τ π τ τ − = − ∫ , (2.56) где ( ) K τ – оценка корреляционной функции. Систематическая составляющая погрешности оценки (2.56) обуслов- лена тем, что длительность T отрезка наблюдаемого процесса конечна. Можно показать, что математическое ожидание оценки спектра мощности двухкомпонентного ССП имеет вид { } { } 2 2 2 ( ) ( ) ( ) 0 0 4 2 2 2 ( 2 ) ( 2 ) 0 0 2 A T P f Sa f f T Sa f f T A T Sa f f T Sa f f T ⎡ ⎤ ⎡ ⎤ = − + + + ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤ + − + + ⎣ ⎦ ⎣ ⎦ π π π π График этой функции для двух значений параметра T показан на рис. 2.21. На нижнем графике параметр T вдвое больше, чем на верхнем. При построении графиков приняты следующие обозначения: ' ' ' 2 r f N r f T r f N t B f f Δ = Δ Δ = = Δ Δ , 0 0 f T r = %===== функц. коррел. полигармонич. сигнала ====== figure; r = 0:200; fT1 = r/10; fT2 = r/5; ( ) P f 2 4 A 2 2 A 0 0 f 2 0 f f Рис. 2.20. График функции (2.55) 64 f0T1 = 6; f0T2 = 12; P1 = (sinc(fT1 – f0T1)).^2 + 2*(sinc(fT1 – 2*f0T1)).^2; P2 = (sinc(fT2 – f0T2)).^2 + 2*(sinc(fT2 – 2*f0T2)).^2; subplot(2,1,1); plot(r,P1); % fT = r/10 title(‘fT = r/10’); grid on; subplot(2,1,2); plot(r,P2); % fT = r/5 title(‘fT = r/5’); grid on; Рис. 2.21. График математического ожидания оценки спектра мощности двух- компонентного ССП Верхний график (рис. 2.21) соответствует случаю 0 , 6; 10 r f T r = = нижний – 0 ; 12 5 r f T r = = . Контрольные вопросы 1. Что такое ряд Фурье и каковы условия Дирихле? 2. Как определяются ППФ и ОПФ? 3. Какая связь между коэффициентами Фурье и спектром сигнала? 4. Каковы свойства ДПФ? 5. Что такое спектрограмма сигнала и как она вычисляется в MATLAB? 65 6. Что такое аддитивная смесь сигнала и как происходит процесс ее построения в MATLAB? 7. Какова функция корреляции? 8. Какова оценка спектра мощности сигнала с учетом преобразова- ния Фурье? |