Моделирование выборочных данных суммой экспоненциальных функций Лекция 12

1. Digital Signal Processing

Лекция 12
DSP

2. Моделирование выборочных данных суммой экспоненциальных функций (метод Прони)

• Введение
• Метод наименьших квадратов Прони
• Модифицированный метод наименьших квадратов Прони
• Спектральная интерпретация метода Прони
• Примеры спектральных оценок на основе метода Прони
DSP

3. Цифровой спектральный анализ

Введение
Метод Прони — это метод моделирования последовательности отсчетов
данных с помощью линейной комбинации экспоненциальных функций был
предложен французским ученым Гаспаром Рише (бароном де Прони) в 1795
году. Он пришел к выводу, что законы, описывающие расширение газов,
могут быть представлены с помощью суммы экспоненциальных функций и
предложил метод для интерполяции данных своих измерений, основанный
на согласовании параметров экспоненциальной модели с измеренными.
Исходная процедура точно согласует экспоненциальную кривую содержащую
p затухающих экспонент Ajexp(ajt), каждая из которых характеризуется двумя
параметрами Aj и aj, с 2p результатами измерений данных. Современный
вариант метода Прони обобщен на модели, состоящие из затухающих
синусоид (комплексных экспонент), кроме этого, в нем используется
процедура оценивания параметров модели по методу наименьших квадратов
для приближенной подгонки модели в тех случаях, когда число точек данных
N>2p – превышает минимально необходимое их число для определения
параметров p экспонент. Эта процедура получила название обобщенного
метода Прони.
DSP

4. Цифровой спектральный анализ

Метод Прони, строго говоря, не является методом спектрального
оценивания. Тем не менее он тесно связан с алгоритмами линейного
предсказания по методу наименьших квадратов, используемыми при
спектральном оценивании на основе моделей авторегрессии. В отличие от
стохастических параметрических АРСС – моделей, в методе Прони для
аппроксимации данных используется детерминированная экспоненциальная
модель, вычисление спектральной плотности энергии (СПЭ) которой и
составляет суть спектральной интерпретации метода Прони.
Заметим, что периодограммную оценку спектральной плотности мощности
(СПМ) можно считать эквивалентной среднеквадратичной аппроксимации
данных с помощью ряда Фурье, т.е. гармонического набора комплексных
синусоид. Так для N отсчетов данных x[0],…,x[N–1], разделенных интервалом
T, аппроксимирующая последовательность имеет вид
N 1
где
xˆ[n] am exp( j 2 f m nT ),
n 0, N 1,
m 0
1
am
N
N 1
x[n] exp( j 2 f
n 0
m
nT ),
m 0, N 1,
если коэффициенты am определяются из условия минимизации суммарной
среднеквадратичной ошибки аппроксимации
N 1
n 0
DSP
x[n] xˆ[n] ,
2
а частоты fm гармонически связаны между собой:
fm m
NT
,
m 0, N 1.

5. Цифровой спектральный анализ

Таким образом, в гармонической модели частоты и число синусоид задаются
заранее, поэтому необходимо оценивать только мощность этих
синусоидальных составляющих на основе соотношения
2
1 N 1
am x[n] exp( j 2 mn / N ) , m 0, N 1,
N n 0
соответствующего вычислению СПМ дискретной периодограммы.
В свою очередь, негармоническая модель, используемая в методе Прони,
требует оценки не только мощности, но и числа синусоид и их частот.
С другой стороны, гармоническая модель наблюдаемых данных
предполагает их периодическое продолжение вне интервала наблюдения,
что далеко не всегда соответствует реальному поведению процесса и
связано с отрицательным проявлением эффектов окна. В негармонической
модели Прони искажающее действие окна исключено, поэтому точность
оценки СПМ по сравнению со стандартным подходом на основе
преобразования Фурье может значительно улучшиться.
2
DSP

6. Цифровой спектральный анализ

Метод наименьших квадратов Прони
Предположим, что имеется N комплексных отсчетов данных x[n], .
Обобщенный метод Прони позволяет оценить x[n] с помощью набора из p
экспоненциальных функций с произвольными амплитудами Ak, частотами fk,
фазами k и коэффициентами
затухания ak:
p
p
(6.66)
xˆ[n] Ak exp[( k j 2 f k )( n 1)T j k ] hk zkn 1 ,
k 1
где
hk Ak exp( j k ),
k 1
zk exp[( ak j 2 f k )T ].
Заметим, что hk – это комплексная амплитуда, представляющая собой
независящий от времени параметр, а zk – это комплексная экспонента,
которая описывает параметр, зависящий от времени.
Отыскание параметров Ak, fk, k, ak и r, минимизирующих сумму квадратов
N
2
ошибок
[ n] ,
где
n 1
p
[n] x[n] xˆ[n] x[n] hk zkn 1 ,
k 1
DSP
представляет трудную нелинейную задачу аппроксимации по методу
наименьших квадратов. Для ее решения могут быть использованы
различные итеративные алгоритмы, требующие больших вычислительных
затрат и не всегда сходящиеся к глобальному минимуму. Альтернативное
субоптимальное решение, в котором используются решения двух систем
линейных уравнений, основано на методе наименьших квадратов Прони.

7. Цифровой спектральный анализ

Ключевым моментом метода Прони является тот факт, что функция xˆ[ n]
является решением некоторого однородного линейного разностного
уравнения с постоянным коэффициентами, вид которого можно
определить следующим образом. Определим сначала полином Ф(z),
корнями которого являются экспоненты zk:
p
(6.67)
( z ) ( z zk ),
k 1
для которого справедливо эквивалентное
представление в виде степенного
p
ряда
(6.68)
( z ) a[m]z p m ,
m 0
с комплексными коэффициентами a(m), для которых a(0)=1. Осуществляя в
выражении (6.66) сдвиг индекса от n до n–m и домножая обе его части на
коэффициент a(m), получим
p
a[m]xˆ[n m] a[m] hk zkn m 1 ,
где 1 n m N.
k 1
Записывая аналогичные произведения a[0] xˆ[ n ],…, a[m 1]xˆ[n m 1] и
суммируя p+1 произведение, получаем
p
p
p
a[m]xˆ[n m] h a[m]z
m 0
где p 1 n N.
DSP
i 1
i
m 0
n m 1
i
,
(6.69)

8. Цифровой спектральный анализ

Осуществляя в (6.66) подстановку
получим уравнениеp
zin m 1 zin p 1 zip m ,
p
a[m]xˆ[n m] h z
m 0
i 1
n p 1
i i
p
a[m]z
m 0
p m
i
0,
(6.70)
в котором равенство нулю следует из факта, что вторая сумма в (6.69) есть
полином (zi), вычисленный в точке, соответствующей одному из его корней.
Таким образом, для аппроксимирующей последовательности x ( n )
справедливо разностное уравнение
p
xˆ[n] a[m]xˆ[n m],
(6.71)
m 1
определенное для p 1 n N.
Полином (z), ассоциированный с этим разностным уравнением, называют
характеристическим, а его корни zk определяют экспоненциальные
параметры в (6.66).
Если учесть, что разность между реальными измеренными данными x[n] и
их аппроксимацией xˆ[ n ] есть величина ошибки [n], так, что
(6.72)
x[n] xˆ[n] [n],
n 0, N 1,
то подстановка
(6.72) в (6.71) дает
соотношение
p
p
p
x[n] a[m]xˆ[n m] [n] a[m]x[n m] a[m] [n m], (6.73)
m 1
DSP
m 1
m 0
где использовано равенство xˆ[n m] x[n m] [n m] и положено a[0]=1.

9. Цифровой спектральный анализ

Соотношение (6.73) можно трактовать как моделирующее процесс x[n] с
помощью модели авторегрессии и скользящего среднего (АРСС–модели) с
идентичными АР– и СС– параметрами, возбуждаемой шумовым процессом
e[n].
В обобщенном методе Прони вводится новая ошибка
p
e[n] a[m] [n m], n p 1, , N ,
(6.74)
m 0
и моделирующее x[n] уравнение принимает вид
p
x[n] a[m]x[n m] e[n],
(6.75)
m 1
который идентичен уравнению для ошибки линейного предсказания вперед
e[n] с коэффициентами фильтра линейного предсказания a[m].
Выбирая параметры a[m] из условия минимизации суммы квадратов
N
ошибок линейного предсказания
e[n]
n p 1
2
, мы тем самым сводим
нелинейную задачу минимизации суммы квадратов ошибок аппроксимации
N
[n] , к линейной системе нормальных ковариационных
2
n 1
DSP
уравнений линейного предсказания.
Таким образом, первый этап обобщенной процедуры Прони сводится к
процедуре оценивания АР–параметров a[m] на основе ковариационного
метода линейного предсказания с оценкой числа экспонент p по правилам
выбора порядка АР–модели.

10. Цифровой спектральный анализ

Второй этап процедуры состоит в нахождении корней zi полинома (6.68)
сформированного из коэффициентов линейного предсказания a[m]
(факторизация полинома).
Наконец, когда значения параметров z1 , z p были определены с помощью
линейного предсказания по методу наименьших квадратов и факторизацией
полинома, то аппроксимация xˆ[ n ] , описываемая уравнением (6.66),
становится линейной относительно оставшихся неизвестных параметров
h1 , , h p . Матричная форма соотношения (6.66) имеет вид
(6.76)
Xˆ Z H ,
где (N p) матрица Z, (p 1) вектор H и (N 1) вектор X определяются
выражениями:
1
z
1
Z
N 1
z1
DSP
1
z1
z 2N 1
1
h1
h
z1
2
, H ,
z pN 1
h p
xˆ[1]
x[2]
ˆ
Xˆ
.
x
[
N
]
ˆ
(6.77)

11. Цифровой спектральный анализ

2
N
Минимизируя сумму квадратов ошибок
x[n] xˆ[n] по каждому параметру
n 1
hk,получаем следующее комплексное нормальное уравнение для их
определения
H
H
(6.78)
где (p p) матрица
выражениями
Z Z H Z
Z Z
X,
и (N 1) вектор отсчетов данных определяются
x[1]
11 1 p
x[2]
,
(6.79)
Z H Z ,
X
p1 pp
x[ N ]
H
z z N 1
N 1
n
j k
, z j z k 1,
jk z j z k kj , jk z j z k 1
n 0
N,
z j z k 1.
Система уравнений (6.79) решается относительно неизвестных параметров
hk, например, по методу Холецкого.
По найденным параметрам zi и hi находятся коэффициенты затухания i,
частоты fi, амплитуды Ai и начальные фазы i с помощью соотношений:
i ln zi , f i arctg Im zi / Re zi / 2 T ,
(6.80)
DSP
Ai hi ,
i arctg Im hi / Re hi .

12. Цифровой спектральный анализ

Модифицированный метод наименьших квадратов Прони
Для процессов, состоящих из p вещественных незатухающих ( =0) синусоид
и шума, разработан модифицированный вариант метода Прони. В этом
случае модель (6.66)p можно записать в виде
p
xˆ (n) 2 A cos 2 f (n 1)T
h z n 1 h ( z ) n 1 , (6.81)
k 1
k
k
k
где 1 n N, hk Ak exp j k , zk exp j2 f k
T .
k 1
k k
k
k
Заметим, что zk являются величинами единичного модуля с произвольными
частотами, которые появляются комплексно сопряженными парами до тех
пор, пока f k 0 либо f k 1 / 2T .
Соответствующий (6.67) и (6.68) характеристический полином для этого
случая имеет вид
p
2p
(6.82)
( z )
( z z )( z z ) a[k ]z 2 p k 0,
i
i
i 1
k 0
где a[0]=1, а a[k] – вещественные коэффициенты. Поскольку корни
полинома (6.82) имеют единичный модуль и появляются в виде комплексно
сопряженных пар, то уравнение (6.82) должно быть инвариантным
относительно подстановки z -1 вместо z:
1
z ( z ) z
2p
2p
2p
a[k ]z
k 0
DSP
k 2 p
2p
a[k ]z k 0.
k 0
(6.83)

13. Цифровой спектральный анализ

Сравнивая (6.82) и (6.83) можно видеть, что a[k ] a[2 p k ], при 0 k p ,
где a[0]=a[2p]=1. Следовательно, требование на комплексно сопряженные
пары корней единичного модуля реализуются посредством наложения на
коэффициенты полинома свойства симметрии относительно центрального
элемента. Таким образом, однородное линейное разностное уравнение,
для которого (6.81) рассматривается в качестве его решения, имеет вид
p
a[ p]xˆ[n p] a[ p k ]xˆ[n p k ] a[ p k ]xˆ[n p k ] 0 (6.84)
k 1
для 2 p 1 n N .
С учетом равенства a[ p k ] a[ p k ],k 1, p и введением коэффициентов
g p [k ] a[ p k ] / a[ p]
(6.84) преобразуется в уравнение
p
xˆ[n p] g p [k ] xˆ[n p k ] xˆ[n p k ] 0.
k 1
DSP
(6.85)

14. Цифровой спектральный анализ

В модифицированном методе Прони на первом этапе ошибка линейного
предсказания, определяемая уравнением (6.75), заменяется ошибкой
линейного сглаживания, использующей как предшествующие, так и
последующие значения отсчетов данных:
p
e[n] x[n] g p [k ] x[n k ] x[n k ] ,
k 1
определенной на интервале p 1 n N p (используются только
имеющиеся отсчеты данных), и минимизируется сумма квадратов ошибок
сглаживания
N p
p
e[n] .
2
n p 1
Получаемые нормальные уравнения для определения коэффициентов gp[k]
соответствуют уравнениям модифицированного ковариационного метода
оценки АР–параметров.
DSP

15. Цифровой спектральный анализ

Спектральная интерпретация метода Прони
Процедура Прони обычно завершается вычислением оценок параметров,
определяющих амплитуды, частоты, фазы и коэффициенты затухания.
Однако возможно вычисление и “спектра” Прони, соответствующего
экспоненциальной аппроксимации xˆ[ n ] . При этом можно получать разные
варианты спектров в зависимости от принятых допущений относительно вида
колебаний вне интервала наблюдения.
Одно из допущений состоит в том, что сумма экспонент дискретного времени
в соотношении (6.66) определяется на интервале – <n< как односторонняя
функция вида
p
hk z kn , n 0;
(6.86)
xˆ[n 1]
k 1
0,
n 0.
Z–преобразование от (6.86) имеет вид
hk
,
Xˆ 1 ( z )
1
k 1 1 z k z
для |z|>|zk|<1 (затухающие экспоненты), спектральная плотность энергии
Прони для рассматриваемой модели будет определяться выражением
p
Sˆ1 ( f ) TXˆ 1 exp( j 2 fT ) ,
2
DSP
которое определено на интервале частот 1 / 2T f 1 / 2T .
(6.87)

16. Цифровой спектральный анализ

Другим возможным допущением является двусторонняя функция вида
p
n
n 0,
hk z k ,
xˆ 2 [ n] p k 1
(6.88)
n
hk z k* , n 0,
k 1
* 1
z
exp(
T
j
2
f
T
)
и
z
k exp( k T j 2 f k T ).
k
k
где k
Такое определение (6.88) обеспечивает симметрию затухающей части
экспоненты относительно начала координат, а Z–преобразование от X 2 ( n )
примет вид
p
1
1
X 2 ( z ) hk
1
1
1 zk z
k 1
1 z k* z
( z k 1 / z k* ) z 1
hk
,
*
1
*
2
1 ( zk 1 / zk )z ( zk / zk )z
k 1
p
1
z
z zk 1.
k
для
Спектральная плотность энергии в этом случае определится как
2
(6.89)
Sˆ2 ( f ) TX 2 exp( j 2 fT ) ,
для 1 / 2T f 1 / 2T .
DSP

17. Цифровой спектральный анализ

В общем случае спектр S 2 ( f ) имеет более острые пики, чем спектр Sˆ1 ( f ).
Высота пиков СПЭ определяется величиной (2Ak/ k)2, а ширина пиков (по
уровню 6 дБ) величиной k/ , поэтому разрешение по частоте меняется в
зависимости от затухания. Для незатухающей синусоиды ( =0),
определенной на бесконечном временном интервале S 2 ( f ) имеет
бесконечное значение на частоте синусоиды и ведет себя подобно
дельта–функции.
DSP

18. Цифровой спектральный анализ

Примеры спектральных оценок на основе метода Прони
На рис. 31- 34 представлены спектральные оценки комплексной
64-точечной тест-последовательности Марпла, полученные на основе
обобщенного метода Прони с использованием односторонней и
двусторонней моделей при значениях порядка моделей 16 и 30
соответственно. На рис.35-37 представлены спектральные оценки
действительной 64-точечной тест-последовательности, полученные на
основе модифицированного метода Прони с использованием односторонней
модели и линейчатого спектра при значениях порядка моделей 16 и 30
соответственно. Действительная тест-последовательность, содержащая три
синусоиды с частотами 0.1, 0.2, и 0.21 при отношении сигнал/шум +10дб,
+20дб, +30дб соответственно, где отношение сигнал/шум определяется как
отношение мощности каждой синусоиды к полной мощности аддитивного
окрашенного шума, полученного фильтрацией белого гауссова шума,
заимствована из обзора Кея и Марпла. Полоса шумового процесса
центрирована относительно частоты 0.35.
DSP

19. Цифровой спектральный анализ

Рисунок 31. Обобщенный метод Прони. Порядок модели 16, односторонний спектр.
DSP

20. Цифровой спектральный анализ

Рисунок 32. Обобщенный метод Прони. Порядок модели 16, двусторонний спектр.
DSP

21. Цифровой спектральный анализ

Рисунок 33. Обобщенный метод Прони. Порядок модели 30, односторонний спектр.
DSP

22. Цифровой спектральный анализ

Рисунок 34. Обобщенный метод Прони. Порядок модели 30, двусторонний спектр.
DSP

23. Цифровой спектральный анализ

Рисунок 35. Модифицированный метод Прони. Порядок модели 16, односторонний спектр.
DSP

24. Цифровой спектральный анализ

Рисунок 36. Модифицированный метод Прони. Порядок модели 16, линейчатый спектр.
DSP

25. Цифровой спектральный анализ

Рисунок 37. Модифицированный метод Прони. Порядок модели 30, линейчатый спектр.
DSP