Диссертация (1149280), страница 16
Текст из файла (страница 16)
Подстановка в формулу для корреляционной функции rP (τ ) приτ = qP даётN−1∞XX2πi1|k| 1Sn∗ Gn,k e 2N τ n .αFrP (τ ) =21 − αF2N n=−Nk=−∞Упростим полученное выражение. Очевидно, Gn,0 = Sn . При k ≥ 1 из определенийследует, чтоS−n = Sn∗ ,ПоэтомуN−1Xn=−NSen∗ Gn,−k e2πiτn2NG∗−n,−k = Gn,k .!∗=N−1Xn=−N2πiSen∗ Gn,k e 2N τ n .Отсюда в формуле для функции rP (τ ) можно обойтись только положительными индексами k:1rP (τ ) =1 − αF2"#∞N −1N−1XX2πi11 X e 2 2πi τ nαFk|Sn | e 2N + 2 ReSen∗ Gn,k e 2N τ n .2N n=−N2Nk=1n=−N111Разложение экспонентыВведём следующие обозначения:N/2−1XFn,j =t=−N/2ste2tNj2πie− 2N tn ,−N ≤ n ≤ N − 1,N −12πi1 X e∗Sn Fn+ℓ,j e 2N τ n ,2N n=−Nρℓ,j (τ ) =j ≥ 0,ℓ ≥ 0,j ≥ 0,|τ | < N,причём Fn,j продолжается периодически по n, так что последняя свёртка циклическая.
Изопределения следует, что Fn,0 = Sen есть ДПФ от set .Теорема 8. При |τ | < N1rP (τ ) =1 − αF2"#∞∞N −1jXX(−πix/2)1 X e 2 2πi τ n2FkαFk|Sn | e 2N + 2 Reρℓ2F k ,j (τ ) .2N n=−Nj!j=0k=1Доказательство. Пусть k ≥ 1. В показателе функции2πizPkt = e− 2N 2F ktвыделим целую и дробную части относительно ДПФ на 2N отсчётов:2F k = ℓ2F k + x2F k ,|x2F k | ≤ 1/2.Сомножитель с дробной частью разложим в ряд Тейлора в нуле, что даётN/2−1XGn,k =t=−N/2=N/2−12πist zPkt e− 2N tne=Xt=−N/2N/2−1∞X(−πix2F k /2)j Xj=0Таким образом,1rP (τ ) =1 − αF2"j!t(n+ℓ2F k )− 2πi2Nt=−N/2setset e2tNjj∞Xπit1− x2F kj!Nj=02πie− 2N t(n+ℓ2F k ) .#N −1∞∞jXX1 X e 2 2πi τ n(−πix/2)2F k|Sn | e 2N + 2 Reρℓ2F k ,j (τ ) .αFk2N n=−Nj!j=0k=1Отметим, что |πix2F k /2| ≤ π/4 < 0.8, поэтому слагаемые быстро убывают по j.112Окончательное приближениеПриближение функции rP основано на выборе конечного множества пар неотрицательных индексов M = {(ℓ, j)}, для которых вычисляются величины ρℓ,j (τ ).Следствие 3.
Общее количество преобразований Фурье, необходимых при вычислении приближенияrbP,M (τ ) =1 ρ0,0 (τ ) + 2 Re1 − αF2для −N ≤ τ < N, равноXαFk(k,j):(ℓ2F k ,j)∈Mj(−πix2F k /2)ρℓ2F k ,j (τ ) .j!Nfft = |M| + Jmax + 1.Доказательство. Выражение для rbP,M (τ ) следует из теоремы 8. ДПФ на 2N отсчётоввыполняется для каждого элементы (ℓ, j) множества M для расчёта ρℓ,j (·), а также длякаждого j ≤ Jmax при вычислении F·,j , где Jmax — наибольшее значение j в парах из M.Поэтому общее количество ДПФ равноNfft = |M| + Jmax + 1,где |M| — количество элементов в множестве M.При фиксированной первой компоненте ℓ для повышения качества аппроксимации естьсмысл включать в множество M все пары (ℓ, j) при 0 ≤ j ≤ J(ℓ) − 1, где J(ℓ) — количествокоэффициентов многочлена Тейлора в разложении экспоненты.
Далее будем считать, чтоэто условие выполнено.В частности, условие J(ℓ) = 0 означает, что ни одной пары вида (ℓ, j) не входит в M.Очевидно, что|M| =∞XJ(ℓ).ℓ=0Исходная задача поиска целых периодов основного тона P была сведена по лемме 12к расчёту значений φ0 (P ) по корреляционным функциям rP (τ ). Полученные приближенияфункций rP дают возможность рассчитать следующие приближения целевых функций φ0 (P ).113Пусть P — натуральное число и F = N/P > 1. Тогда⌊F ⌋Xφb0 (P ) = rbP,M (0) + 2q=1rbP (qP )⌊F ⌋X1ρ0,0 (0) + 2ρ0,0 (qP )=1 − αF2q=1X(−πix2F k /2)jk+αF 2 Reρℓ2F k ,j (0)j!(k,j):(ℓ2F k ,j)∈M⌊F ⌋j XX(−πix2F k /2)k+αF 4 Reρℓ2F k ,j (qP ).j!q=1(k,j):(ℓ2F k ,j)∈MОценка погрешностиМаксимальное число ℓ, для которого (ℓ, 0) ∈ M, обозначим L.Пусть P — целое число и F = N/P . Максимальное целое число k, для которого 2F k ≤L + 0.5, обозначим K(P ).
Очевидно,P (2L + 1)K(P ) =.4NТеорема 9. Погрешность приближения функции φ0 оценивается сверху следующим образом:N/2−1|φ0 (P ) − φb0 (P )| ≤гдеγP (t) =K(P )X2 1 − αF2 k=1Xt=−N/2γP (t)|est | |est | + 2⌊ Pt + F2 ⌋Xq=1|est−qP | ,J(ℓ2F k )1 + |αF | 1 πt x2F k + |αF |K(P )+1|αF |k,J(ℓ2F k )! N|1 + αF zPt |2Доказательство. Погрешность оценкиrP (τ ) − rbP,M (τ ) = δα (τ ) + δx (τ )складывается из погрешности усечения геометрической прогрессии и погрешности многочлена Тейлора:1δα (τ ) =2 Re1 − αF2∞Xk=K(P )+1K(P )αFkN −12πi1 X e∗Sn Gn,k e 2N τ n ,2N n=−NN−1XX1k 1τnbn,k )e 2πi2NαFSen∗ (Gn,k − G,2 Reδx (τ ) =21 − αF2Nn=−Nk=1114где при τ = qPN/2−1Gn,k =Xt=−N/2N/2−12πiset zPkt e− 2N tnX=t=−N/2J(ℓ2F k )−1N/2−1bn,k = =GXt=−N/2Xt(n+ℓ2F k )− 2πi2Nset e2πi2πiset e− 2N t(n+ℓ2F k ) e− 2N x2F k t ,j=01j!jπit− x2F k .NПогрешность усечения геометрической прогрессии снова сворачивается, как геометрическая прогрессия:N/2−1N −12πi1 X e∗ 2πi τ n XSn e 2Nset zPkt e− 2N tn2N n=−Nk=K(P )+1t=−N/2N/2−1N−1(K(P )+1)tK(P )+1XX2πi2πi1zαFSen∗2 Ree− 2N tn e 2N τ n .set P=2t1 − αF2N n=−N1 − αF zP∞X1δα (τ ) =2 Re1 − αF2αFkt=−N/2ДПФ сохраняет скалярное произведение с точностью до постоянного множителя, аN/2−12πiN −1есть ДПФ от вектора (st−τ )t=−N/2 .
Поэтому(Sn e− 2N τ n )n=−NN/2−1K(P )+1(K(P )+1)tXαFzPδα (τ ) =esse2Ret−τ t1 − αF21 − αF zPtt=−N/2Обозначим для краткости κ = K(P ) + 1. ТогдаαFκδα (τ ) =1 − αF2αFκ=1 − αF2Следовательно,|δα (τ )| ≤N/2−1Xt=−N/2N/2−1Xt=−N/2zP−κtzPκt+set−τ set1 − αF zPt1 − αF zP−t(κ−1)tset−τ est2 Re(zP− αF zPκt )|1 − αF zPt |2N/2−12|αF |K(P )+1 X |est−τ | |est |.1 − |αF ||1 − αF zPt |2t=−N/2Погрешность многочлена Тейлора для экспоненты с чисто мнимым показателем имеетвидn−1 ix X (ix)j |x|n.e −≤j! n!j=01152πiПри 1 ≤ k ≤ K(P ) применим это неравенство к функции e− 2N x2F k t .
Это приводит к верхнейграницеJ(ℓ2F k )N/2−1K(P )XX πt12k x2F k |α||es||es|.|δx (τ )| ≤Ft−τt1 − αF2 k=1J(ℓ2F k )! Nt=−N/2В итоге при q ≥ 0N/2−1|rP (qP ) − rbP,M (qP )| ≤гдеγP (t) =K(P )X2 1 − αF2 k=1Xt=−N/2|est−qP | |est |γP (t),J(ℓ2F k ) πt1 + |αF | 1 x2F k + |αF |K(P )+1,|αF |kJ(ℓ2F k )! N|1 + αF zPt |2и это выражение не зависит от q.Погрешность приближения функции φ0 оценивается сверху следующим образом:N/2−1|φ0 (P ) − φb0 (P )| ≤Xt=−N/2⌊ Pt + F2 ⌋γP (t)|est | |est | + 2Xq=1|est−qP | .Гарантированное оцениваниеN/2−1Точность оценивания функции φ0 сравнивается с нормой вектора sw = (sw,t )t=−N/2 . ПоN/2−1определению, set = sw,t wt = st wt2 , где (wt )t=−N/2 — окно Ханнинга.Найдём минимальное число λ > 0, для которогоN/2−1Xt=−N/2+1γP (t)|est | |est | + 2⌊ Pt + F2 ⌋Xq=1|est−qP | ≤ λksw k2для любого сигнала s. Левая часть является квадратичной формой, обозначим её f (sw ).Искомое число λ есть норма матрицы этой квадратичной формы.Форма f (sw ) обладает рядом особенностей, упрощающих её исследование.
Сумму можно переставить таким образом, чтобы в каждом слагаемом оказались произведения сигналовс индексами из одного класса эквивалентности по модулю P . Выполним это преобразование.Для каждого вычета k = 0,1, . . . , P − 1 минимальное число из промежутка −N/2 + 1 ≤t ≤ N/2−1, сравнимое с k, обозначим t0k , а количество чисел из этого промежутка, сравнимых116с k, обозначим Nk . Тогдаf (sw ) =P −1 Nk −1XXγP (t0k + nP )est0k +nPset0k +nP + 2k=0 n=0nXq=1set0k +(n−q)P!.Отсюда видно, что исходная задача гарантированного оценивания распадается на P независимых подзадач:Nk −1XγP (t0k+ nP )est0k +nPn=0st0k +nP + 2enXq=1st0k +(n−q)Pe!≤ λkNk −1Xs2w,t0 +nP ,0 ≤ k ≤ P − 1.kn=0Затем в качестве общего множителя λ можно гарантированно выбрать наибольшее из λk , которое затем максимизируетс и по P . В практических задачах можно выбирать и взвешенноесреднее значений λk с учётом интенсивностей сигнала.Каждая частная задача решается при помощи следующей леммы, в которую подставляется xn = set0k +nP , cn = wt0k +nP , dn = γP (t0k + nP ), N = Nk − 1.NЛемма 17.
Пусть заданы положительные числа (cn )Nn=0 и (dn )n=0 . Минимальное значениеλ, для которогоNXdn xnxn + 2n=0nXxn−qq=1!≤λNXx2nn=0c2n∀ x = (xn )Nn=0 ,,совпадает с нормой самосопряжённой матрицы A = (ai,j )Ni,j=0 с элементамиai,j = ci cj dmax{i,j} ,0 ≤ i,j ≤ N.Доказательство. Левая часть неравенства в лемме равнаNXn=0dn xnxn + 2nXq=1xn−q!x0∗ x1 = .. . xNd0d1· · · dNx0 d1 d1 · · · dN x1 .. .. ... .
. ...... . xNdN dN · · · dNВыполним замену переменной xn = cn yn . Тогда неравенство в лемме равносильно следующему:y ∗Ay ≤ λkyk2.Матрица A самосопряжённая по построению, поэтому минимальное значение λ совпадает сkAk.Если P ≥ N/3, что предполагается в постановке задаче, то Nk ≤ 2 при всех k и P .Таким образом, лемму 17 надо применить для случаев N = 0, N = 1 и N = 2. В первыхдвух случаях задача решается аналитически.
В последнем случае получается кубическоеуравнение, для которого можно подобрать аналитическое приближённое решение.117В итоге, можно сформулировать следующее утверждениеТеорема 10. Для каждого P ∈ [Pmin , Pmax ]|φ0 (P ) − φb0 (P )| ≤ λksw k2 ,гдеλ=max λk ,0≤k≤P −1λk = kAk,ai,j = ci cj dmax{i,j} ,иcn = wt0k +nP ,dn = γP (t0k + nP ).0 ≤ i,j ≤ Nk − 1,118Глава 4. Экспериментальные результатыПредложенный алгоритм оценивания ЧОТ назван «метод минимизации дисперсии шума» (МДШ), или «Noise Variance Minimization» (NVM). На основе полученных результатовсоздана реализация алгоритма поиска ЧОТ в пакете прикладных программ Matlab [72]. Экспериментальная проверка показала возможность применения результатов к решению прикладных задач.В текущей главе продемонстрированы результаты работы алгоритма, а так же сравнение с существующими алгоритмами оценивания ЧОТ.Сверхточное определение значения ЧОТ является специальной задачей, избыточнойво многих приложениях.
Так как в полигармонических моделях речевого сигнала частотывходящих в неё гармоник кратны ЧОТ, то и ошибка, возникающая при определении ЧОТкратно распространяется на старшие гармоники. На рис. 4.1 представлены: участок спектраисходного сигнала, участок спектра восстановленного сигнала с точным значением ЧОТ иучасток спектра восстановленного сигнала, где ошибка определения ЧОТ равна 1Гц.Хорошо видно, что ошибка в 1Гц значительно смещает положение кратных гармоник.В области частот начиная от 4кГц восстановленный и исходный сигналы уже совершенно непохожи друг на друга.4.1Демонстрация работы алгоритма определения ЧОТПродемонстрируем все шаги работы алгоритма быстрого оценивания ЧОТ на конкретном сигнале.
В качестве примера возьмём участок голосового сигнала (рис. 4.2).спектр исходного сигналаточная оценка ЧОТнеточная оценка ЧОТ33.544.55частота, кГц5.566.5Рисунок 4.1 — Влияние ошибки определения ЧОТ на спектр восстановленного сигнала1190.30.20.10-0.1-0.2-0.3-0.4020406080100120140NРисунок 4.2 — Участок голосового сигналаНа первом этапе алгоритма требуется зафиксировать множество M (см. раздел 3.3.3).В этом примере выберем два варианта:% M_1l_required = [0 1 2J_num= [3 3 35 % M_2l_required = [0 1 220];J_num=[1 0 01];3 4 5 6 7 8];3 3 3 2 2 2];3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 193 3 2 2 2 2 12101101001Стоит отметить, что сложность алгоритма напрямую зависит от выбора множества M (следствие 3). При выборе M1 и M2 сложность двух вариантов алгоритмов совпадает.Следующий шаг — приближённое вычисление φ0 (P ).















