Основы цифровой обработки сигналов
Параметрические и непараметрические методы спектрального анализа
Непараметричесике методы
В непараметрических методах для расчета спектра сигнала используется только информация, заключенная в отсчетах сигнала, без каких-либо дополнительных предположений. Рассмотрим два таких метода – периодограмму и метод Уэлча.
Периодограмма
Периодограмма – это оценка спектральной плотности мощности, полученная по отсчетам одной реализации случайного процесса.
Периодограмма рассчитывается по формуле:
где – число отсчетов, – частота дискретизации.
Периодограмма не является состоятельной оценкой спектральной плотности мощности, т.к. дисперсия такой оценки сравнима с квадратом ее математического ожидания. С ростом числа отсчетов значения периодограммы начинают все больше флуктуировать.
В качестве примера оценим в Engee спектральную плотность мощности сигнала, представляющего собой сумму двух синусоид с частотами 100 и 300 Гц. Для этого воспользуемся функцией periodogram, которая содержится в библиотеке DSP. Входными параметрами этой функции являются исходный сигнал x, логическая переменная
onesided, определяющая тип спектра (true – односторонний, false – двусторонний), длина БПФ nfft, частота дискретизации fs и тип используемого окна window. Выходом функции является объект y, содержащий вектор дискретных частот , вызываемый командой freq(y), и вектор спектральной плотности мощности , вызываемый командой power(y). Построим графики сигнала во времненой и в частотной области.
using Plots, DSP;
fs = 4000;
dt = 1/fs;
t = [0:dt:1;];
x = 2*cos.(2*pi*100*t) + cos.(2*pi*300*t);
plot(t[1:200], x[1:200], xlabel="t", ylabel="сигнал", legend=false)
y = DSP.periodogram(x, onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(freq(y)[1:400], power(y)[1:400], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
✏️Задание 1
Постройте в Engee аккорд ля мажор первой октавы, представляющий собой сумму синусоид с частотами 440, 550 и 660 Гц. Оцените спектральную плотность мощности этого сигнала с помощью функции periodogram. Постройте график спектральной плотности мощности.
Решение
using Plots, FFTW;
fc = [440 550 660]';
fs = 8000;
dt = 1/fs;
t = [0:dt:0.1;];
x = cos.(2*pi*fc[1]*t) + cos.(2*pi*fc[2]*t) + cos.(2*pi*fc[3]*t);
x = x/3;
y = DSP.periodogram(x, onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(freq(y)[1:100], power(y)[1:100], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
Метод Уэлча
При вычислении периодограммы по длинному фрагменту случайного сигнала она оказывается весьма изрезанной. Для уменьшения этой изрезанности необходимо применить какое-либо усреднение.
Метод, предложенный Уэлчем, использует весовую функцию и разбиение сигнала на перекрывающиеся сегменты. Вычисления по методу Уэлча организуются следующим образом:
- Вектор отсчетов сигнала делится на перекрывающиеся сегменты. Как правило, используется перекрытие 50%.
- Каждый сегмент умножается на используемую весовую функцию.
- Для взвешенных сегментов вычисляются модифицированные периодограммы.
- Периодограммы всех сегментов усредняются.
Теперь решим предыдущую задачу в Engee методом Уэлча с помощью функции welch_pgram, которая содержится в библиотеке DSP. Входными параметрами этой функции являются исходный сигнал x, число сегментов n, величина перекрытия сегментов noverlap (в отсчетах), логическая переменная
onesided, определяющая тип спектра (true – односторонний, false – двусторонний), длина БПФ nfft, частота дискретизации fs и тип используемого окна window. Выходом функции является объект y, содержащий вектор дискретных частот , вызываемый командой freq(y), и вектор спектральной плотности мощности , вызываемый командой power(y).
using Plots, DSP;
fs = 4000;
dt = 1/fs;
t = [0:dt:1;];
x = 2*cos.(2*pi*100*t) + cos.(2*pi*300*t);
n=div(length(x),8);
y = welch_pgram(x, n, div(n,2), onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(freq(y)[1:200], power(y)[1:200], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
✏️Задание 2
Аналгично Заданию 1, постройте в Engee аккорд ля мажор первой октавы, представляющий собой сумму синусоид с частотами 440, 550 и 660 Гц. Оцените спектральную плотность мощности этого сигнала с помощью функции welch_pgram.
Решение
using Plots, DSP;
fc = [440 550 660]';
fs = 8000;
dt = 1/fs;
t = [0:dt:0.1;];
x = cos.(2*pi*fc[1]*t) + cos.(2*pi*fc[2]*t) + cos.(2*pi*fc[3]*t);
x = x/3;
n=div(length(x),8);
y = welch_pgram(x, n, div(n,2), onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(freq(y)[1:100], power(y)[1:100], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
Параметрические методы
Использование параметрических методов подразумевает наличие некоторой математической модели анализируемого случайного процесса. Спектральный анализ сводится к решению оптимизационной задачи, т.е. поиску таких параметров модели, при которых она наиболее близка к реально наблюдаемому сигналу.
Авторегрессионная модель
В авторегрессионной модели сигнал формируется путем пропускания дискретного белого шума через рекурсивный фильтр -го порядка. Этот метод сводится к определению коэффициентов модели заданного порядка , оценке мощности белого шума и расчету спектральной плотности мощности по формуле:
Для определения коэффициентов модели производится минимизация ошибки линейного предсказания сигнала. Этот метод заключается в том, что сигнал пропускается через нерекурсивный фильтр. Взвешенную сумму предыдущих отсчетов входного сигнала называют линейным предсказанием следующего входного отсчета, а выходной сигнал – ошибкой предсказания.
Авторегрессионные методы дают хорошие результаты, когда спектр анализируемого сигнала имеет четко выраженные пики. К таким сигналам относится сумма нескольких синусоид с шумом.
Метод MUSIC
Метод MUSIC (MUltiple SIgnal Classifiction) предназначен для спектрального анализа сигналов, представляющих собой сумму нескольких синусоид с белым шумом. Целью спектрального анализа подобных сигналов обычно является не расчет спектра как такового, а расчет псевдоспектра, который представляет собой набор частот и уровней гармонических составляющих.
В основе метода лежит анализ собственных чисел и собственных векторов корреляционной матрицы сигнала.
Пусть сигнал представляет собой сумму комплексных экспонент, начальные фазы которых случайны, с белым шумом:
где – отсчеты дискретного белого шума, – период дискретизации, – начальные фазы (независимые случайные величины, равномерно распределенные на интервале от 0 до ).
Корреляционная функция такого сигнала будет иметь вид:
где – единичная импульсная функция, равная 1 при и 0 в остальных случаях, – дисперсия шума.
Далее из отсчетов корреляционной функции формируется корреляционная матрица.
Спектр комплексной экспоненты, рассчитанный аналитически, представляет собой дельта-функцию, расположенную на соответствующей частоте. Чтобы получить бесконечные выбросы на частотах , псевдоспектр в методе MUSIC рассчитывается следующим образом:
где – -й элемент -го собственного вектора корреляционной матрицы.