CC BY
122
ЭЛЕКТРОНИКА "\ Ж 1 ^Ж 1 ■ ▼
УДК 621.385
О РЕАЛИЗАЦИИ МЕТОДА ПРОНИ ДЛЯ СЛАБО ЗАШУМЛЕННЫХ
ВЫБОРОК
ГРИЦУНОВ А.В., КОПОТЬ М.А., ЛЕБЕДЕВ О.Г, ФРОЛОВА Т.И.
Обсуждается возможность использования гармонической интерполяции (исходного метода Прони) для гармонического анализа выходного сигнала в спектральных моделях электронных приборов. Она более устойчива по сравнению с методом наименьших квадратов Прони (гармонической аппроксимацией) при нахождении частот гармонических компонент для слабо зашумленных выборок. Исследуется частотное разрешение предлагаемой методики для различного числа близко расположенных гармонических составляющих.
1. Введение
Одной из составных частей алгоритма спектрального моделирования вакуумных и полупроводниковых электронных приборов [1] является гармонический анализ выходного сигнала, полученного в виде дискретной временной выборки мгновенных значений напряжения или тока на выходных клеммах прибора.
Классическое дискретное преобразование Фурье, использовавшееся ранее для этой цели [2], обладает предельно высокой устойчивостью (робастностью) оценок и дает информацию не только об амплитуде, но и о фазе каждой из базисных гармонических компонент, что существенно при анализе усилителей (так как в этом случае необходимо вычислять “электронное смещение фазы” сигнала). Однако ему присущи недостатки в виде ограниченного частотного разрешения на коротких временных последовательностях, получаемых при моделировании переходных процессов, а также явления “утечки” [3]. Поэтому развитие спектрального моделирования стало возможным только после освоения параметрических методов гармонического и спектрального анализа [4].
Очевидно, что статистические модели спектрального оценивания (например, авторегрессионная) плохо пригодны для решения указанной задачи, поскольку гармоническое разложение выходного сигнала прибора требует сохранения суммарной мощности процесса при переходе к дискретным гармоническим компонентам. Кроме того, такие модели обычно дают оценку спектральной плотности мощности сигнала без учета его фазы. Поэтому анализ выходных сигналов следует проводить с
использованием детерминированных моделей. Наиболее распространены среди них гармоническая интерполяция и аппроксимация функций, известные еще как исходный метод Прони и метод наименьших квадратов Прони соответственно [4].
Из общих соображений следует, что при декомпозиции зашумленных функций, какими обычно являются результаты моделирования приборов методами Монте-Карло, более целесообразна аппроксимация временных выборок, так как она непосредственно обеспечивает статистически усредненные оценки. Однако практика применения метода наименьших квадратов Прони показала, что при достаточно “чистом” (т.е. незашумленном) выходном сигнале, который обычно имеет место в “холодном” приборе, а также в первые моменты переходного процесса данный метод имеет тенденцию к ненормальному завершению. Оно происходит, как правило, на первом этапе алгоритма (решение системы нормальных ковариационных уравнений) и вызвано тем, что предполагаемый порядок модели существенно превышает фактическое количество гармонических компонент в анализируемой функции.
В этой связи имеет смысл попытаться вместо гармонической аппроксимации использовать для анализа слабо зашумленных выборок выходных сигналов гармоническую интерполяцию (исходный метод Прони), которая, как предполагается, должна обладать меньшей критичностью к выбору порядка модели. В данной статье описывается реализация такого подхода, который может оказаться полезным и в других приложениях (например, в радиолокации).
2. Постановка задачи
Определена временная выборка исходной функции fl = f(lAt), где At — временной интервал дискретизации; 1=0, 1, ..., L—1; L > 3 — длина выборки, кратная 3. Функция f может быть комплексной или вещественной.
Для обоих случаев необходимо заменить заданную выборку последовательностью M незатухающих комплексных экспонент, где M=2L/3 — порядок модели, т.е.
M M .
fl = Z Am Єхр(ІЮтШ) = X Amzm • (1)
m=1 m=1
Здесь rom и Am — круговая частота и комплексная амплитуда m-й неизвестной комплексной экспоненты соответственно; zm = exp(iromAt) — z-преобразование исходной временной последовательности.
M уравнений (1) с неизвестными величинами zm
и Am могут быть записаны в обобщенной матричной форме:
||z||A = f , (2)
РИ, 2002, № 4
27
где ||z|| — (M x M) матрица коэффициентов zm*, записываемая в виде:
( 0 0 о ^
z1 z2 •• zM
1 1 1
z|| = z1 z2 - zM
M-1 M-1 M-1
1 z1 z2 • zM )
A и f — векторы Am и f (здесь l=0...M—1).
Система M/2 линейных разностных уравнений для M/2 коэффициентов bj...Ьм/2 имеет вид:
M/2 *
Z (bmfl+m + b mfl_m) = “fl , (7)
m=1
где l=M/2...L—M/2—1. Чтобы преобразовать эту систему к пригодному для решения виду, коэффициенты b и fдолжны быть явно выражены через их вещественные и мнимые части, т.е. b = br + ib1; f = fr + if1. В результате система (7) может быть переписана в виде:
Для вещественной функции f уравнения (1) и (2) не изменяются, но множества величин zm и Am могут быть разделены на два равных подмножества из M/2 членов, каждое из которых является комплексным сопряжением другого. Чтобы обеспечить унификацию алгоритма, желательно рассматривать спектр вещественной функции как двусторонний (т.е. с отрицательными и положительными частотами) . Алгоритм гармонической интерполяции далее будет описан только для комплексной выборки, так как модификация его для вещественной последовательности очевидна.
3. Реализация метода
Первой стадией алгоритма гармонической интерполяции является нахождение коэффициентов характеристического полинома. Этот полином с M нулями zj... zm может быть записан так:
'M/2
Z (fri+m + fri-m)brm + m=1
M/2 . . .
+ Z (—fh+m + f1l-m)b1m = -f ri ; m=1
M/2 . .
Z (f.l+m + f.l-m)brm + (8)
m=1 M/2
+ £ (frl+m - frl-m)b.m = f .
* m=1
Для решения системы (8) может быть рекомендован алгоритм гауссова исключения с частичным выбором ведущего элемента [5]. В итерационном уточнении, как правило, нет необходимости, поскольку точность расчета мгновенных значений выходного сигнала в моделях электронных приборов обычно невысока.
MM
П (z - zm) =ZamzM-m, (3)
m=1 m=0
где am — коэффициенты многочлена (комплексные или вещественные, в зависимости от типа функции f) и ао = 1. Используя процедуру, примененную в [4] для модифицированного метода наименьших квадратов Прони, можно записать систему из M/2 линейных разностных уравнений для M коэффициентов а1 ...ам :
M/2
Z(aM/2-mfl+m + aM/2+mfl-m) = “aM/2fl > (4)
m=1
здесь l=M/2...L—M/2—1. Дополнительным условием является комплексная сопряженность коэффициентов:
am = aM • a м-m , (5)
позволяющая в два раза уменьшить количество неизвестных переменных.
Если ввести нормализованные коэффициенты bm = aM/2-m /aM/2, характеристический полином (3) может быть переписан так:
M/2
Z bmz
M/2+m
m=1
+ z
M/2
M/2*
+ Z b mz
m=1
M/2-m
(6)
В выражении (6) учтено, что в соответствии с условием (5) aM/2+m /aM/2 = (aM/2-m /aM/2) .
Вторым этапом алгоритма является факторизация характеристического полинома (6) с нахождением
всех его нулей z 1... z м . Для полиномов с комплексными и вещественными коэффициентами наиболее пригодны методы Дженкинса-Трауба [6] и Дженкинса [7] соответственно. Типичная длина выборки выходного сигнала для нестационарных и спектральных моделей электронных приборов равна 30...255 (т.е. ближайшее меньшее, чем 25 ...28 значение, кратное 3). Поэтому порядок модели M обычно равен 20...170. К сожалению, оригинальные алгоритмы Дженкинса-Трауба и Дженкинса работают удовлетворительно только со степенями полиномов, не превышающими 49 и 99 соответственно. Для больших порядков точность нахождения последних нулей недопустимо снижается.
Причина некорректного функционирования оригинальных алгоритмов факторизации очевидна с учетом понижения ими степени полинома после нахождения следующего нуля. Новые коэффициенты полинома накапливают ошибки округления. Вследствие высокой чувствительности корней многочлена к значениям этих коэффициентов [5] последние найденные нули могут быть достаточно далеки от соответствующих корней исходного (нередуцированного) полинома.
Радикальное устранение этого недостатка в контексте известных алгоритмов, вероятно, невозможно. Тем не менее, его влияние может быть ослаблено
28
РИ, 2002, № 4
путем уменьшения среднеквадратичной погрешности вычисления корней. Соответствующий метод описан в [8]. Он заключается в использовании двух критериев сходимости (“сильного” и “слабого”) вместо первоначального одного. Это позволяет увеличить среднеквадратичную точность факторизации полинома ценой увеличения времени счета (примерно в 2-3 раза). Модифицированные алгоритмы Дженкинса-Трауба и Дженкинса работают удовлетворительно с полиномами степени до 256 включительно.
Третьей стадией алгоритма гармонической интерполяции является решение системы линейных уравнений (2). Найденные нули zj... zM должны быть предварительно нормализованы к единичному модулю и затем отсортированы по частоте с удалением кратных. Соответственно, порядок модели M может быть уменьшен. Метод решения системы (2) такой же, как на первом этапе, т.е. гауссово исключение с частичным выбором ведущего элемента и без итерационного уточнения.
Поскольку на каждом этапе описанного алгоритма могут происходить ошибки из-за плохой обусловленности исходной выборки, важно обеспечить надежное его функционирование для устранения возможности аварийного завершения процесса моделирования. Одним из способов может быть итеративное снижение порядка модели M на 2 с соответствующим уменьшением длины выборки L на 3 до тех пор, пока будет получен успешный результат. Добавление аддитивного белого шума со среднеквадратическим значением около—100-120 дБ от амплитуды анализируемой функции также может быть полезным для обеспечения устойчиво -сти работы алгоритма со слабо зашумленными выборками. В целом же, как и предполагалось, решение системы (8) существенно реже подвержено ненормальному завершению, чем решение системы нормальных ковариационных уравнений на соответствующем этапе метода наименьших квадратов Прони.
Демонстрационная программа на Фортране — 90 для гармонической интерполяции комплексной выборки приведена на рис.1. Функция Number_Of_Harmonics возвращает количество гармонических компонент для анализируемой последовательности после удаления кратных нулей характеристического полинома. Функции Harmonic_Frequency и Harmonic_Amplitude возвращают относительную линейную частоту (в единицах частоты дискретизации) и комплексную амплитуду любой компоненты соответственно. Это упрощенная версия профессиональной программы без итеративного снижения порядка модели для плохо обусловленных выборок и указанной выше модификации алгоритма факторизации полинома. Поэтому максимальное значение L для данной программы равно 72. В программе использованы подпрограммы библиотеки IMSL DLSLRG, DZPOCC и DLSLCG.
4. Проверка метода
Приведенная демонстрационная программа была применена к 64-точечной тест-последовательности из книги [4], а также к аналогичной последовательности, у которой имелись только четыре одиночные гармонические компоненты с относительными линейными частотами —0.15, 0.10, 0.20, 0.21 и модулями амплитуд 0.1, 0.1, 1.0 и 1.0 соответственно, а шумовые составляющие отсутствовали. Во втором случае линейные частоты и модули комплексных амплитуд исходных гармоник были восстановлены так:
иц /2л = -0.15000000..., Ю18 /2л = 0.10000000..., ®23 /2л = 0.20000000..., ю24/2л = 0.20999999...,
A11 = 0.099999999...; A18 = 0.10000000...; A23| = 1.0000000...; A24| = 0.99999999...,
т.е. для чисто гармонической выборки точность интерполяции вполне удовлетворительна. Однако применительно к исходной зашумленной последо -вательности результаты иные:
ю15/2л = -0.15001436..., ю22/2л = 0.099987216..., ю25/2л = 0.20004258..., ю26/2л = 0.20996366...,
A15I
A22
A25
A26
0.093785705...; 0.089861647...; 1.2376770...; 0.045267775....
Как видно, погрешность нахождения частот гармонических компонент увеличилась незначительно. В то же время точность оценки комплексных амплитуд теперь совершенно неудовлетворительна. Следовательно, как и предполагалось, описанная методика мало пригодна для выборок с большим отношением сигнал/шум. Одной из причин этого является неполное использование анализируемой последовательности на заключительном этапе алгоритма, так как не более 2/3 исходной выборки задействовано во время решения системы (2). Приемлемым выходом может быть замена гармонической интерполяции функции fi с использованием системы (2) ее гармонической аппроксимацией на всех L точках исходной последовательности, т.е. комбинация исходного подхода Прони для этапов формирования и факторизации характеристического полинома с методом наименьших квадратов Прони на этапе оценки амплитуд гармоник.
В описанных расчетах использован Compaq Visual Fortran Professional Edition v. 6.5 для Intel 386 с исходными настройками для конфигурации Win32 Debug. Заметим, что в [4] содержится ошибка: мнимая часть элемента тест-последовательности с номером 58 должна быть положительна.
Для иллюстрации возможностей предлагаемой методики на рис. 2 показана зависимость частотного разрешения описанного алгоритма от числа гармо -нических компонент во временной последовательности. В спектре исходной незашумленной комплексной выборки с L=255 присутствует группа гармоник со средней относительной линейной частотой 1/16. Она состоит из h=2...8 эквидистант-
РИ, 2002, № 4
29
module Constants
complex (kind=8), parameter :: Imagl = (0.0D+00,1.0D+00) real (kind=8), parameter :: Pi = 3.141592653589793D+00 real (kind=8), parameter :: Pi2 = 6.283185307179586D+00 end module Constants
module Global_Data
complex (kind=8) : A (1:48) ! комплексные амплитуды
complex (kind=8) : B (1:24) ! коэффициенты полинома
integer (kind=4) : L ! длина выборки (L)
integer (kind=4) : M ! порядок модели (M)
complex (kind=8) : Z (1:48) ! нули полинома
end module Global_Data
integer (kind=4) function Number_Of_Harmonics ( SSize, F ) use Constants
use Global_Data use IMSL
integer (kind=4) complex (kind=8) real (kind=8) :: real (kind=8) :: integer (kind=4)
intent(in) intent(in) : Freq1 Freq2 Limit Ll Mm Mm1
Offset Switch TempC
allocatable allocatable allocatable :: allocatable :: allocatable ::
SSize
F (0:SSize-1)
integer (kind=4) integer (kind=4) integer (kind=4) integer (kind=4) integer (kind=4) complex (kind=8) complex (kind=8), allocatable :: TempCV (:) complex (kind=8), allocatable :: TempCM (:,:) real (kind=8), allocatable :: TempRV1 (:) real (kind=8), allocatable :: TempRV2 (:) real (kind=8), allocatable :: TempRM (:,:) real (kind=8), parameter :: Near = 1.0D-06 !—Инициализация процедуры: if ( SSize < 3 ) then
print *, 'Ошибка: слишком короткая выборка..
Number_Of_Harmonics = 0 ! ошибка
return end if
размер выборки анализируемая выборка (f) частота гармоники 1 частота гармоники 2 для сортировки Шелла индекс выборки (l) индекс порядка (m) индекс порядка (альтернатив.) для сортировки Шелла для сортировки Шелла временный комплексный скаляр временный комплексный вектор временная комплексная матрица временный веществен. вектор 1 временный веществен. вектор 2 временная веществен. матрица критерий кратного нуля
M = min ( (SSize/3)*2, 48 ) ! порядок модели
L = 3*M/2 ! полезная длина выборки
!—Решение системы линейных уравнений для коэффициентов характеристич. полинома: allocate ( TempRV1(0:M-1), TempRV2(1:M), TempRM(0:M-1,1:M) ) do Ll = M/2, L-M/2-1
TempRV1(Ll-M/2) = - real(F(Ll))
TempRV1(Ll) = - imag(F(Ll)) do Mm = 1, M/2
TempRM(Ll-M/2,Mm) = real(F(Ll+Mm)) + real(F(Ll-Mm))
TempRM(Ll-M/2,M/2+Mm) = -imag(F(Ll+Mm)) + imag(F(Ll-Mm))
TempRM(Ll,Mm) = imag(F(Ll+Mm)) + imag(F(Ll-Mm))
TempRM(Ll,M/2+Mm) = real(F(Ll+Mm)) - real(F(Ll-Mm)) end do end do
call DLSLRG ( N=M, A=TempRM, B=TempRV1, X=TempRV2, LDA=M, IPath=1 ) B(1:M/2) = cmplx ( TempRV2(1:M/2), TempRV2(M/2+1:M), 8 ) deallocate ( TempRV1, TempRV2, TempRM )
!—Факторизация характеристического полинома: allocate ( TempCV(0:M) )
TempCV(M/2+1:M:+1) = B(1:M/2)
TempCV(M/2) = 1.0D+00
TempCV(M/2-1:0:-1) = conjg(B(1:M/2))
call DZPOCC ( NDeg=M, Coeff=TempCV, Root=Z )
deallocate ( TempCV )
!—Нормализация нулей полинома:
Z = Z / abs(Z)
!—Сортировка Шелла нулей полинома:
Offset = M/2
do while ( Offset > 0 )
Limit = M - Offset Switch = 1 do while ( Switch > 0 )
30
РИ, 2002, № 4
Switch = 0 do Mm = 1, Limit
Freql = atan2 ( imag(Z(Mm)), real(Z(Mm)) ) / Pi2 Freq2 = atan2 ( imag(Z(Mm+Offset)), real(Z(Mm+Offset)) ) / Pi2 if ( Freq1 > Freq2 ) then TempC = Z (Mm)
Z(Mm) = Z(Mm+Offset)
Z(Mm+Offset) = TempC Switch = Mm end if end do
Limit = Switch - Offset end do
Offset = Offset / 2 end do
!—Удаление кратных нулей:
Mm = 2
do while ( Mm <= M )
Freq1 = atan2 ( imag(Z(Mm-1)), real(Z(Mm-1)) ) / Pi2 Freq2 = atan2 ( imag(Z(Mm)), real(Z(Mm)) ) / Pi2 if ( abs(Freq2-Freq1) > abs(Freq1)*Near ) then Mm = Mm + 1 else
do Mm1 = Mm, M-1
Z(Mm1) = Z(Mm1+1) end do
M = M - 1 ! новый порядок модели
end if end do
L = 3*M/2 ! новая полезная длина выборки
!—Решение системы линейных уравнений для комплексных амплитуд: allocate ( TempCV(0:M-1), TempCM(0:M-1,1:M) )
TempCV(0) = F(0)
T empCM(0, :) = 1.0D+00 do Ll = 1, M-1
TempCV(Ll) = F(Ll) do Mm = 1, M
TempCM(Ll,Mm) = TempCM(Ll-1,Mm)* Z(Mm) end do end do
call DLSLCG ( N=M, A=TempCM, B=TempCV, X=A, LDA=M, IPath=1 ) deallocate ( TempCV, TempCM )
!—Завершение процедуры:
Number_Of_Harmonics = M return
end function Number_Of_Harmonics
real (kind=8) function Harmonic_Frequency ( NHarm ) use Constants use Global_Data
integer (kind=4), intent(in) :: NHarm ! номер гармоники
if ( NHarm < 1 .or. NHarm > M ) then
print *, 'Ошибка: неправильный номер гармоники...'
Harmonic_Frequency = 0.0D+00 ! ошибка
return end if
Harmonic_Frequency = atan2 ( imag(Z(NHarm)), real(Z(NHarm)) ) / Pi2 return
end function Harmonic_Frequency
complex (kind=8) function Harmonic_Amplitude ( NHarm ) use Constants use Global_Data
integer (kind=4), intent(in) :: NHarm ! номер гармоники
if ( NHarm < 1 .or. NHarm > M ) then
print *, 'Ошибка: неправильный номер гармоники...'
Harmonic_Amplitude = (0.0D+00,0.0D+00) ! ошибка
return end if
Harmonic_Amplitude = A(NHarm) return
end function Harmonic_Amplitude
РИ, 2002, № 4
Рис. 1
31
ных комплексных экспонент единичной амплитуды с нулевыми начальными фазами. На вертикаль -ной оси отложен логарифм отношения частотного разрешения классического метода 2п / ДТ к разрешению описанного алгоритма Дю • ДТ — временная длина выборки. Здесь Дю определяется как минимальное расстояние по частоте между соседними экспонентами при условии, что приведенный метод восстанавливает все их с относительной ошибкой не более нескольких процентов. Несмотря на приближенную оценку этой зависимости, а также существенное снижение частотного разрешения метода при наличии даже незначительного аддитивного шума, преимущество алгоритма гармонической интерполяции перед обычным разложением в ряд Фурье очевидно. 5
5. Заключение
Описана методика гармонического анализа слабо зашумленного выходного сигнала в спектральных моделях вакуумных и полупроводниковых электронных приборов. При достаточно большом отношении сигнал/шум она обеспечивает на 1-3 порядка более высокую разрешающую способность по сравнению с классическими методами и в то же время существенно более устойчива, чем метод наименьших квадратов Прони. Одним из вариантов улучшения ее показателей для зашумленных выборок является комбинация интерполяции фун -кции при вычислении частот гармонических компонент с ее аппроксимацией для оценки их амплитуд и фаз.
Литература: 1. Gritsunov A. V. On spectral modeling of microwave devices // Telecomm. and Radio Engineering. 2001. Vol. 55, No. 8. P. 98—102. 2. Грицунов A.B., Галаган
A.B. Об использовании уравнения возбуждения второго порядка при моделировании автогенераторов со скрещенными полями // Радиотехника и электроника. 1989. Т. 34, № 8. С. 1719-1722. 3. ОтнесР, Эноксон Л. Прикладной анализ временных рядов. М.: Мир, 1982. 428 с. 4. Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения. М.: Мир, 1990. 584с. 5. Форсайт Дж, Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980. 280 с. 6. Jenkins M. A., Traub J. F. Algorithm 419: Zeros of a complex polynomial // Comm. of ACM. 1972. Vol. 15. P. 97-99. 7. Jenkins M. A. Algorithm 493: Zeros of a real polynomial // ACM Trans. on Math. Software. 1975. Vol. 1. P. 178-189. 8. GritsunovA.V. The spectral and harmonic analysis in microwave electronics // Proc. Int. Conf. “MIKON-2002”. Gdansk, Poland. May 2002.
Поступила в редколлегию 23.04.2002
Рецензент: д-р техн. наук, проф. Лучанинов А.И.
Грицунов Александр Валентинович, канд. техн. наук, доцент кафедры микроэлектроники, электронных приборов и устройств ХНУРЭ. Научные интересы: компьютерное моделирование СВЧ-приборов со скрещенными полями. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (0572) 40-93-62. E-mail: [email protected].
Копоть Михаил Андреевич, младший научный сотрудник кафедры физических основ электронной техники ХНУРЭ. Научные интересы: компьютерный анализ СВЧ-приборов. Адрес: Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (0572) 40-94-84.
Лебедев Олег Григорьевич, канд. техн. наук, доцент кафедры радиолокационных систем Харьковского института военно-воздушных сил. Научные интересы — компьютерный анализ и синтез современных радиолокационных систем. Адрес: Украина, 61165, Харьков, ул. Клочковская, 228, тел. (0572) 30-82-14.
Фролова Татьяна Ивановна, младший научный сотрудник кафедры физических основ электронной техники ХНУРЭ. Научные интересы: численное моделирование СВЧ-приборов М-типа. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (0572) 40-94-84.
32
РИ, 2002, № 4
