DSP Лекция 12 Digital Signal Processing

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

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

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

DSP Цифровой спектральный анализ Таким образом, в гармонической модели частоты и число синусоид задаются заранее, поэтому необходимо оценивать только мощность этих синусоидальных составляющих на основе соотношения соответствующего вычислению СПМ дискретной периодограммы. В свою очередь, негармоническая модель, используемая в методе Прони, требует оценки не только мощности, но и числа синусоид и их частот. С другой стороны, гармоническая модель наблюдаемых данных предполагает их периодическое продолжение вне интервала наблюдения, что далеко не всегда соответствует реальному поведению процесса и связано с отрицательным проявлением эффектов окна. В негармонической модели Прони искажающее действие окна исключено, поэтому точность оценки СПМ по сравнению со стандартным подходом на основе преобразования Фурье может значительно улучшиться.

DSP Цифровой спектральный анализ Метод наименьших квадратов Прони Предположим, что имеется N комплексных отсчетов данных x[n],. Обобщенный метод Прони позволяет оценить x[n] с помощью набора из p экспоненциальных функций с произвольными амплитудами A k, частотами f k, фазами k и коэффициентами затухания a k : (6.66) где Заметим, что h k – это комплексная амплитуда, представляющая собой независящий от времени параметр, а z k – это комплексная экспонента, которая описывает параметр, зависящий от времени. Отыскание параметров A k, f k, k, a k и r, минимизирующих сумму квадратов ошибок где представляет трудную нелинейную задачу аппроксимации по методу наименьших квадратов. Для ее решения могут быть использованы различные итеративные алгоритмы, требующие больших вычислительных затрат и не всегда сходящиеся к глобальному минимуму. Альтернативное субоптимальное решение, в котором используются решения двух систем линейных уравнений, основано на методе наименьших квадратов Прони.

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

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

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

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

DSP Цифровой спектральный анализ Минимизируя сумму квадратов ошибок по каждому параметру h k,получаем следующее комплексное нормальное уравнение для их определения (6.78) где (p p) матрица и (N 1) вектор отсчетов данных определяются выражениями (6.79) Система уравнений (6.79) решается относительно неизвестных параметров h k, например, по методу Холецкого. По найденным параметрам z i и h i находятся коэффициенты затухания i, частоты f i, амплитуды A i и начальные фазы i с помощью соотношений: (6.80)

DSP Цифровой спектральный анализ Модифицированный метод наименьших квадратов Прони Для процессов, состоящих из p вещественных незатухающих ( =0) синусоид и шума, разработан модифицированный вариант метода Прони. В этом случае модель (6.66) можно записать в виде (6.81) где,,. Заметим, что z k являются величинами единичного модуля с произвольными частотами, которые появляются комплексно сопряженными парами до тех пор, пока либо. Соответствующий (6.67) и (6.68) характеристический полином для этого случая имеет вид (6.82) где a[0]=1, а a[k] – вещественные коэффициенты. Поскольку корни полинома (6.82) имеют единичный модуль и появляются в виде комплексно сопряженных пар, то уравнение (6.82) должно быть инвариантным относительно подстановки z -1 вместо z : (6.83)

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

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

DSP Цифровой спектральный анализ Спектральная интерпретация метода Прони Процедура Прони обычно завершается вычислением оценок параметров, определяющих амплитуды, частоты, фазы и коэффициенты затухания. Однако возможно вычисление и спектра Прони, соответствующего экспоненциальной аппроксимации. При этом можно получать разные варианты спектров в зависимости от принятых допущений относительно вида колебаний вне интервала наблюдения. Одно из допущений состоит в том, что сумма экспонент дискретного времени в соотношении (6.66) определяется на интервале – |z k |

DSP Цифровой спектральный анализ Другим возможным допущением является двусторонняя функция вида (6.88) где и Такое определение (6.88) обеспечивает симметрию затухающей части экспоненты относительно начала координат, а Z –преобразование от примет вид для Спектральная плотность энергии в этом случае определится как (6.89) для

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

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

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

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

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

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

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

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

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