CC BY
194
УДК 519.254+621.317.7:621.391
А.А. Львов, В.П. Глазков, В.П. Краснобельмов, Р.С. Коновалов, М.А. Соломин
ОЦЕНИВАНИЕ ПАРАМЕТРОВ КВАЗИГАРМОНИЧЕСКИХ СИГНАЛОВ МЕТОДОМ
МАКСИМАЛЬНОГО ПРАВДОПОДОБИЯ
Предложена методика оценивания параметров квазигармонического сигнала в случае, когда закон изменения его частоты неизвестен, основанная на методе максимального правдоподобия. Приведёны практические примеры, подтверждающие эффективность предлагаемой методики.
Квазигармонический сигнал, метод максимального правдоподобия, амплитуда, частота, фаза, сверхвысокочастотный сигнал, пьезорезонансный датчик давления
A.A. L'vov, V.P. Glazkov, V.P. Krasnobel'mov, R.S. Konovalov, M.A. Solomin
ESTIMATION OF QUASI-HARMONIC SIGNAL PARAMETERS BY THE MAXIMUM LIKELIHOOD METHOD
The technique for parameter estimation of the quasi-harmonic signal with unknown frequency based on the maximum likelihood method is proposed. The technique applications confirming its efficiency are given.
Quasi-harmonic signal, maximum likelihood method, amplitude, frequency, phase, microwave signal, piezo-resonant pressure transducer
Введение
Благодаря широкому практическому применению задача оценки параметров компонент сигнала из конечного числа зашумленных дискретных измерений до сих пор является активной областью исследований. Самый сложный случай оценивания, когда частота сигнала неизвестна. Несмотря на то, что существует множество схем оценки частоты, они могут быть классифицированы как непараметрические и параметрические методы.
К непараметрическим методам относятся преобразование Фурье и периодограммный анализ
[1], которые используют только информацию, содержащуюся в отсчетах анализируемого сигнала, без каких-либо дополнительных предположений о его. Существенным недостатком таких методов является невозможность получения одновременно высокого разрешения по времени и частоте. Это создает определенные проблемы при анализе нестационарных сигналов, сложных колебаний с широким спектром и низким соотношением сигнал/шум. Альтернативой являются параметрические методы, использование которых подразумевает наличие некоторой математической модели анализируемого сигнала.
Например, такие широко распространенные параметрические методы, как Кумерсана-Тафтса
[2] и Прони [3], позволяют достичь более высокой разрешающей способности, однако являются чувствительными к влиянию шумов. Другой метод, основанный на оценивании параметров сигнала, называемый, ESPRIT (Estimation of Signal Parameters via Rotational Invariance Technique) [4] существенно менее чувствителен к точности аппроксимации фоновых шумов и обеспечивают более высокое разрешение. Метод построения канонических форм матриц (Matrix Pencil) [5] похож на метод ESPRIT. Главное отличие состоит в том, что ESPRIT работает с сигнальным подпространством, определенным из корреляционной матрицы, в то время как Matrix Pencil обрабатывает данные напрямую, что даёт существенную экономию в числе операций [5]. В [6] было показано, что Matrix Pencil и метод Кумерсана-Тафтса являются частными случаями линейного предсказания, но первый метод является более эффективным при вычислении и менее чувствительным к шуму, чем второй. В тех случаях, когда предполагается, что модель является близким приближением к реальному объекту, параметрические методы обеспечивают более точные оценки и высокое разрешение, чем непараметрические [1,7]. Однако все перечисленные параметрические методы являются субоптимальными. Наилучшую оценку среди прочих параметрических методов для фиксированного набора данных и
147
базовой вероятностной модели даёт метод максимального правдоподобия (ММП), поскольку оценки, полученные с его помощью, являются состоятельными, асимптотически эффективными и асимптотически несмещёнными [1,8]. Но основной недостаток этого метода связан с трудностями расчёта оценок.
Тем не менее, в некоторых практически важных случаях квазигармонического сигнала, когда его частота меняется «достаточно медленно», удаётся найти оценки неизвестных параметров по ММП по дискретным отсчётам данного сигнала с помощью относительно несложных вычислительных процедур, рассмотренных в работе.
1. Постановка задачи
В некоторых приложениях [9-12] возникает задача, подобная следующей. Исследуется изменяющийся во времени сигнал вида:
z = A sin {w0 + v(t)]t + j}+ C +d ,
где: z - измеряемая величина; W0 = 2pf0 - известная круговая частота; v(t) - неизвестный параметр частотной модуляции; A и j- неизвестные амплитуда и фаза сигнала; C - постоянное смещение, обусловленное тем, что перед оцифровкой в аналого-цифровом преобразователе (АЦП) сигнал должен быть сделан неотрицательным; ё - аддитивный шум, который полагается распределенным по нормальному закону с нулевым математическим ожиданием и неизвестной дисперсией а2.
Пусть в дискретные равноотстоящие друг от друга моменты времени ti...,t„, были получены соответствующие измерения величины z, которые будем рассматривать как вектор данных
Z = (zj, z2,. ., zn)T . Кроме этого предполагается выполнение условия квазигармонического сигнала, то есть V() << wo. Поэтому считается, что в течение времени наблюдения T = nAt, где At - интервал времени между соседними отсчётами сигнала, величина v(t) остаётся постоянной. Тогда с учётом выражения, приведённого выше, имеем:
z j = z (tj )= A sin [(Wo + v)-tj + j]+ C + d j , (j = 1, n).
Задача заключается в нахождении оценок неизвестных величин A, j, C и v Смещение C тоже считается неизвестным, потому что в усилителях, используемых для масштабирования измеряемого сигнала перед АЦП возможен паразитный дрейф нуля. Решение подобной четырёх параметрической задачи, предлагаемое стандартом IEEE Std. 1057 [13], находится в результате итерационной вычислительной процедуры, основанной на решении системы уравнений (1) по методу наименьших квадратов (МНК) [8]. Алгоритм обеспечивает оценки при наличии гауссовского шума с точностью, близкой к нижней границе Крамера-Рао. В более сложных условиях, когда зашумленные отсчеты проведены только на части периода сигнала, алгоритм дает худшую точность оценивания [14]. В случае, когда мгновенная частота сигнала известна с небольшой погрешностью, возможно более быстрое и точное оценивание искомых параметров.
2. Процедура оценивания по ММП и её анализ
В работах [10, 12] показано, что при малых абсолютных значениях величины частотной модуляции v по сравнению с основной частотой W0, используя известные приближённые соотношения cos Vtj »1, sin vtj » ntj и сделав следующую замену переменных:
Ху = sin W)tj, x2j = cos W)tj, x3j = tj cos W)tj, x4j = -tj sin Wtj, X5j = 1, (j = 1, n),
q1 = A cosj, q2 = A sinj, q3 = Avcosj, q4 = Avsinj, q5 = С,
2)
нелинейная задача оценивания (1) сводится к линейной с квадратичным уравнением связи на вновь введённые оцениваемые переменные Q = (q¡, q2, q3, q4, q5)T:
zj = qxj + q2x2j + q3x3j + q4x4j + q5x5j + dj, (j = 1, n), 3)
qq- q2 q3 =0 4)
После этого решение системы (3) с ограничением (4) по ММП может быть найдено с помощью итерационной вычислительной процедуры [15]:
Q(k+1) = Q(k) Qk) GqQc) (xT X)-1GTQ(k) (5)
Q = Q - 2Q(%T (X t X^GQ® ( ) GQ ,
где к - номер итерации, « - обозначение обратной матрицы, X и С - матрицы плана эксперимента и квадратичной формы (4) соответственно:
^ =
x11 x12 x13 x14 x15 x21 x22 x23 x24 x25
"0 0 0 1 0"
0 0 -10 0
, G = 0 -1 0 0 0
1 0 0 0 0
0 0 0 0 0
_хп1 хп2 хп3 хп4 хп5 а нулевое приближение задаётся по МНК
Q(0)=(хт х )-1(хт z). (6)
В монографии [15] показано, что полученные оценки коэффициентов Q являются в линейном приближении несмещенными и эффективными как линейные оценки максимального правдоподобия,
причем, как правило, не требуется более двух итераций. Зная оценку вектора Q = (<71, 72, 1з, 44,15) из (5)-(6), с помощью (2) находятся оценки искомых параметров
A
я2 + Я2,
n
Л 2 -2
l<?3 + Я4
)/(я2 + Я2 )
с
Is,
(p-
где Q = arctg
% ^ Яз J
(7)
Q, при q3 > 0, q4 > 0 Q + ж, при q3 < 0 Q + 2ж, при q3 > 0, q4 < 0
Оценки параметров A, j и пиз (7) уже будут смещёнными, поскольку получаются в результате нелинейных преобразований над промежуточными переменными Q, но в [15] показано, что средний квадрат ошибки оценивания получается меньшим, чем без учёта ограничения (4). Таким образом, чисто теоретически предлагаемый метод обладает более высокой точностью, чем методы, описанные во введении, включая алгоритм стандарта IEEE Std. 1057, при более быстрой сходимости итерационного процесса (5)-(6) к решению.
Теперь рассмотрим примеры использования предлагаемого метода.
3. Оценка отклонения частоты сигналов с нелинейной частотной модуляцией
В некоторых приложениях, связанных с генерацией и анализом сверхвысокочастотного сигнала с внутриимпульсной нелинейной частотной модуляцией требуется точно измерять параметры угловой модуляции (частотной и фазовой), в том числе её основной параметр - девиацию частоты. Здесь мгновенная круговая частота сигнала может быть представлена в виде суммы нескольких составляющих [16]:
K
wt ) = w + Zw (t ),
k=1
где С - известная основная частота сигнала, не меняющаяся во времени, йк(0 - «медленные» моду-
лирующие составляющие, для которых выполняется соотношение
Zw (t )
k=1
< w.
При этом зависимость от времени модулирующих компонент определяются следующим равенством [16]:
w (t ) = 2pFkg (t ) cos(2pFkt ), (k = 1K ), 8)
где Fk - неизвестные коэффициенты, подлежащие оцениванию; g(t) - известная функция времени; K -число модуляционных составляющих.
При сделанных предположениях на начальном этапе задача оценивания коэффициентов Fk
K
сводится к описанной выше задаче оценивания закона модуляции v(t) = (t) квазигармониче-
k=1
ского сигнала по его дискретным отсчётам (1). При этом весь интервал наблюдения этого сигнала разбивается на M подынтевалов длины T, причём на каждом из них производится серия из n отсчётов сигнала (1), а величина v считается постоянной в течение времени Т и меняется скачкообразно при переходе на следующий подынтевал.
Тогда для первого подынтервала с учётом (1) можно записать
zj! = z (tj e Ti ) = Ai sin [(w0 + v(Tx ))• t} + j ]+ Сi + 8} , (j = 1, n),
где Ti - первый из подынтервалов равной длины T.
В уравнениях (9) по аналогии с задачей (1) неизвестными являются параметры A1, j, С1 и v(T1), которые могут быть оценены по ММП с помощью процедуры (5)-(7). Интерес для исследователя в этом случае представляет только параметр v(T1), оценка которого запоминается. Далее переходят к следующему подынтервалу T2, на котором частота равна O = 00 + v(?i). Аналогично из (5)-(7)
находится оценка V(T2), и т.д. В общем случае для подынтервала с номером m
z (tj e Tm )= Am sin [(Om_i + V (Tm )) • tj + jm ]+ Ст + , (j = Щ; m = \M)
jm
Повторив процедуру оценивания для всех M подынтервалов, можно получить последовательность оценок величин v(Tm), которые в совокупности дают закон изменения модулирующей частоты исследуемого сигнала.
На втором этапе решается задача нахождения оценок параметров Fk из (8) по дискретным отсчётам v(Tm), m = 1,M , не являющаяся целью настоящей работы.
В работе было проведено компьютерное моделирование процесса нахождения оценок v(Tm)
двумя способами: с помощью алгоритма IEEE Std. 1057 [14] и предлагаемым способом (5)-(7) по ММП. Точность оценивания характеризовалась дисперсией погрешности, которая рассчитывалась по
10000 экспериментам (( = ^ (ym — v /10000, где v - истинное значение отклонения), в зависимо-
m
сти от отношения сигнал/шум, задаваемого как A2/Г2 . На рис. 1 показаны зависимости дисперсии погрешностей оценивания от отношения сигнал/шум, выбираемого из ряда 20, 40, 50, 60, 80 (дБ). Видно, что предлагаемый алгоритм даёт значительный выигрыш в точности, что хорошо согласуется с теоретическими положениями, о которых упоминалось в конце раздела 2 работы. При алгоритм ММП сходился к решению после двух итераций, а IEEE Std. 1057 требовал не менее 5-6 итераций. 4. Измерения выходной частоты пьезорезонансных датчиков давления В настоящее время направление, связанное с разработкой и производством датчиков давления с пьезорезонансными чувствительными элементами интенсивно развивается в России и за рубежом [17]. Преимуществом резонансных датчиков перед датчиками с емкостными и пьезорезистивными чувствительными элементами является высокая точность и стабильность характеристик, малое изменение частоты колебаний при воздействии температуры в широком диапазоне от -60 °С до +150 °С.
Рисунок 1 - Зависимости дисперсии погрешности оценивания параметра частотной модуляции
от отношения сигнал/шум
Предлагается новый метод повышения точности измерения давления пьезорезонансным датчиком, основанный на использовании дискретного преобразования Фурье (ДПФ) [1,3] выходного
сигнала и нахождении оптимальной оценки его частоты по ММП.
На выходе пьезорезонансного датчика давления измеряются отсчёты сигнала неизвестной частоты. Задача заключается в оценивании именно этой частоты, которая связана с измеряемым давлением известным соотношением, зависящим от типа ЧЭ. Математическая модель выходного дискретного сигнала может быть представлена в виде:
z j = z (tj ) = A sin [2pftj + j ]+ C + 8j , (j = 1N), (10)
где все обозначения совпадают с уравнениями (1), а мгновенная частота f несущая информацию об измеряемом давлении, должна быть оценена путём цифровой обработки данных zj. N - число дискретных отсчётов сигнала.
Предлагается метод уточнения мгновенной частоты, основанный на ДПФ полученной последовательности дискретных отсчётов сигнала Z = (z 1, z2,..., zN )T в канале датчика.
Если вычислить ДПФ от данной последовательности, то в случае некратности периодов исследуемого сигнала Tm=1lf и дискретизации Ts=Dt получится дискретный амплитудный спектр сигнала z(t), вид которого схематично показан на рисунке 2. Истинная частота f не будет совпадать ни с одной из дискретных частот в спектре из-за известного явления растекания спектра [18]. Оценить значение истинной частоты на основании полученных отсчётов можно только приближённо. Воспользуемся методикой, описанной в разделе 2, для уточнения данной оценки.
Пусть в качестве начальной оценки мгновенной частоты было выбрано значение f0 - частота спектрального отсчёта максимальной амплитуды. Тогда истинное значение равно f=f0 + n, при этом погрешность оценивания n будет меньше одного дискрета по частоте. При достаточно большом количестве отсчётов n > 256 можно считать, что М << f *.
Рис. 2. ДПФ синусоиды при некратных периодах Тт и Тэ
С помощью алгоритма (5)-(7) можно оценить погрешность частоты V, при этом оценка иско-мои мгновенной частоты сигнала равна: / = / + V .
Теперь из модели (10) можно восстановить синусоиду выходного сигнала, что схематично изображено на рисунке 3,а, где сплошными вертикальными линиями показаны дискретные отсчёты На основании данных отсчётов и полученного периода зондирующего сигнала Тт = 1// можно так изменить период дискретизации Тх, чтобы он стал кратным Тт. Тогда за «новые отсчёты» сигнала г(0 можно взять их предсказанные значения на восстановленной синусоиде (показаны на рисунке 3 пунктирными линиями). По данным предсказанным отсчётам можно ещё раз взять ДПФ. В полученном энергетическом спектре явление растекания уже проявится меньше, поскольку период синусоиды кратен периоду дискретизации. Схематично эта ситуация показана на рисунке 3,б.
Рис. 3 а. Исходные (сплошные линии) Рис- 3 б
и предсказанные (пунктирные линии) отсчёты. Здесь уже центральный (максимальный) отсчёт будет ближе к истинному значению оцениваемой частоты. Приведённые соображения позволяют предложить новый метод уточнения частоты
выходного синусоидального тока пьезорезонансного датчика. Необходимо так изменять период дискретизации, чтобы он становился кратным периоду тока. Эта ситуация возникает, когда центральная спектральная составляющая имеет максимальную высоту относительно своих соседей. Поэтому, изменяя период дискретизации и оценивая параметры модели (10) по описанной методике, необходимо добиться максимальной относительной высоты центральной спектральной составляющей, что соответствует отсутствию растекания спектра, т.е. точной оценке частоты /
В работе было проведено математическое моделирование описанного метода в среде МаЛаЪ 7.0. «Истинная» частота на выходе датчика была взята равной /= 27 кГц, его амплитуда А = 0,7 мА, смещение В = 0, отношение сигнал/шум равно 5 дБ по мощности. Частота дискретизации сигнала в АЦП выбрана /5 = 675 кГц. Результат построения амплитудного спектра с помощью программы быстрого преобразования Фурье из пакета МаЛаЪ показан на рисунке 4.
Из графика видно, что наблюдается растекание спектра. За оценку частоты ^0 выбрано максимальное значение спектрального отсчёта. После этого по описанной методике была оценена частота сигнала. На рис. 5 показан результат вычисления амплитудного спектра после пяти итераций. Оценка
частоты была равна / = 27,01 кГц. Из полученного спектра видно, что явление растекания практически отсутствует, а точность оценки частоты высокая.
Рис. 4. Амплитудный спектр синусоидального сигнала частоты 1 = 27 кГц. Частота дискретизации не кратна измеряемой частоте
Рис. 5. Амплитудный спектр сигнала с рисунка 4 после оценки частоты Г по методике (5)-(7)
Для более высоких отношений сигнал/шум требуется меньшее число итераций, для достижения той же точности.
Полученные результаты говорят о перспективности предлагаемого метода повышения точности интеллектуальных пьезорезонансных датчиков давления.
Заключение
В работе предложен алгоритм оценивания параметров квазигармонического сигнала, частота которого изменяется по нелинейному закону, основанный на методе максимального правдоподобия. Предложена итерационная процедура решения получающихся нелинейных уравенений. Эффектив-152
ность предложенной методики демонстрируют преведённые примеры её использования, связанные с измерением отклонения осноной частоты высокочастотного сигнала и повышением точности измерения пьезорезонансных датчиков давления.
ЛИТЕРАТУРА
1. Stoica, P. Introduction to Spectral Analysis / P. Stoica, R. Moses. Upper Saddle River, NJ: Prentice-Hall, 1997. 427 p.
2. Tufts, D.W. Estimation of Frequencies of Multiple Sinusoids: Making Linear Prediction Perform like Maximum Likelihood / D.W. Tufts, R. Kumaresan // Proc. of IEEE, Vol. 70, № 9, pp. 975-989, Sep. 1982.
3. Марпл-мл., С. Л. Цифровой спектральный анализ и его приложения / С. Л. Марпл-мл. М.: Мир, 1990. 583 с.
4. Roy, R. ESPRIT Estimation of Signal Parameters via Rotational Invariance Techniques / R.Roy, T. Kailath // IEEE Trans. Acoust., Speech, Signal Process., vol. 37, № 7, pp.984-995, Jul. 1989.
5. Hua, Y. Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise / Y. Hua, T. Sarkar // IEEE Trans. Acoust., Speech, Signal Process., vol. 38, no. 4, pp. 814 -824, May 1990.
6. Sarkar, T. Using the Matrix Pencil Method to Estimate the Parameters of a Sum of Complex Exponentials / T. Sarkar, O. Pereira // IEEE Anten. and Propag. Mag., No.1, P. 48-55, Feb. 1995.
7. So, H.C. Linear Prediction Approach for Efficient Frequency Estimation of Multiple Real Sinusoids: Algorithms and Analyses // H.C. So, Kit Wing Chan, Y.T. Chan, K.C. Ho,. IEEE Trans. on Signal Processing, Vol. 53, No. 7, July 2005, P. 2290-2305.
8. Линник, Ю.В. Метод наименьших квадратов и основы теории обработки наблюдений / Ю.В. Линник. М.: ГИФМЛ, 1958, 336 с.
9. Казаков, К.В. Алгоритм двухканального оценивания параметров квазигармонических сигналов / К.В. Казаков, А.А. Львов, В.А. Пыльский // Вестник СГТУ, № 4(43), 2009. - С. 38-41.
10. Львов, А.А. Линейная петлевая схема обработки сигналов с датчиков / А.А. Львов, В.А. Пыльский // Вестник СГТУ, № 2(3), 2004. - С. 102-112.
11. Gureev, V.V. Improvement of the Current Loop Circuit for AC and DC Applications Based on Digital Signal Processing / V.V. Gureev, A.A. L'vov, V.A. Pyl'skiy // Proceed. of the IEEE Instrument. and Measur. Technol. Conf., Sorrento, Italy 24-27 April 2006, P. 1257-1261.
12. L'vov, A.A. Optimal Digital Signal Processing for Current Loop Circuit / A.A. L'vov, V.A. Pyl'skiy // Digest of the IEEE Conf. on Precis. Electromag. Measur., 2006, Torino, Italy, P. 1054-1060.
13. IEEE Standard for Digitizing Waveform Recorders, IEEE Std. 1057, 1994.
14. Händel, P. Properties of the IEEE-STD-1057 Four-Parameter Sine Wave Fit Algorithm / P. Händel // IEEE Trans. on. Instrument. Measur., Vol. 49, No. 6, Dec. 2000, P. 1189-1193.
15. Репин В.Г. Статистический синтез в условиях априорной неопределенности и адаптация информационных систем // В.Г. Репин, Г.П. Тартаковский. М.: Сов.радио, 1977, 242 с.
16. Краснобельмов, В.П. Синтез сигналов с нелинейной внутриимпульсной частотной модуляцией / В.П. Краснобельмов, Д. А. Баринов // Проблемы управления, обработки и передачи информации (АТМ-2011): сб. трудов II Междунар. науч. конф. - Саратов: Издательский Дом «Райт-Экспо», 2012. Т.2. С. 38-40.
17. Поляков, А. Перспективные кварцевые пьезорезонансные датчики давления / А. Поляков, В. Поляков, М. Одинцов // Компоненты и технологии, 2011, №1. C. 45-50.
18. Витязев, В.В. Цифровая частотная селекция сигналов / В.В. Витязев. М.: Радио и связь, 1993. 240 с.
Львов Алексей Арленович -
доктор технических наук, профессор кафедры «Информационные системы и технологии» Саратовского государственного технического университета имени Гагарина Ю.А. Глазков Виктор Петрович -доктор технических наук, профессор кафедры «Автоматизация, управление и мехатроника» Саратовского государственного
Alexey A. L'vov -
Dr. Sc., Professor Department Information Systems and Technologies,
Yuri Gagarin State Technical University of Saratov Viktor P. Glazkov -Dr. Sc., Professor
Department of Automation, Control and Mechatronics,
технического университета имени Гагарина Ю.А.
Краснобельмов Вадим Павлович -
аспирант
Саратовского государственного
технического университета имени Гагарина Ю.А.
Коновалов Роман Станиславович -
аспирант
Саратовского государственного
технического университета имени Гагарина Ю.А.
Соломин Максим Андреевич -
студент
Саратовского государственного
технического университета имени Гагарина Ю.А.
Статья пос
Yuri Gagarin State Technical University of Saratov
Vadim P. Krasnobel'mov -
Postgraduate
Yuri Gagarin State Technical University of Saratov
Roman S. Konovalov -
Postgraduate
Yuri Gagarin State Technical University of Saratov
Maxim A. Solomin -
Undergraduate,
Yuri Gagarin State Technical University of Saratov
пила в редакцию 15.10.14, принята к опубликованию 25.12.14
