Непараметрические методы

Следующие разделы обсуждают периодограмму, измененную периодограмму, валлийский и методы мультизаострения непараметрической оценки, наряду со связанной функцией CPSD, оценкой передаточной функции и функцией когерентности.

Периодограмма

Скрипт Open Live Script

В общих чертах один способ оценить PSD процесса состоит в том, чтобы просто найти преобразование Фурье дискретного времени выборок процесса (обычно делавшимся на сетке с БПФ) и соответственно масштаб в квадрате из результата. Эта оценка называется периодограммой.

Оценка периодограммы PSD сигнала $x_{L} (n)$ из длины L

$P_{x x} (f) = \frac{1}{L F_{s}} {| \sum_{n = 0}^{L - 1} x_{L} (n) e^{- j 2 π f n / F_{s}} |}^{2},$

где Фс является частотой дискретизации.

На практике, фактический расчет $P_{x x} (f)$ может быть выполнен только в конечном числе точек частоты, и обычно использует БПФ. Большинство реализаций метода периодограммы вычисляет $N$ - укажите оценку PSD на частоты

$f_{k} = \frac{k F_{s}}{N}, k = 0, 1, \dots, N - 1 .$

В некоторых случаях расчет периодограммы через Алгоритм бпф более эффективен, если количество частот является степенью двойки. Поэтому весьма распространено заполнить входной сигнал нулями, чтобы расширить его длину к степени двойки.

Как пример периодограммы, рассмотрите следующий xn сигнала с 1001 элементом, который состоит из двух синусоид плюс шум:

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes (row vector)
f = [150;140];            % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
% The three last lines are equivalent to
% xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t));

Оценка периодограммы PSD может быть вычислена с помощью periodogram. В этом случае вектор данных умножается на Окно Хэмминга, чтобы произвести модифицированную периодограмму.

[Pxx,F] = periodogram(xn,hamming(length(xn)),length(xn),fs);
plot(F,10*log10(Pxx))
xlabel('Hz')
ylabel('dB')
title('Modified Periodogram Power Spectral Density Estimate')

Алгоритм

Периодограмма вычисляет и масштабирует выход БПФ, чтобы произвести энергию по сравнению с частотным графиком можно следующим образом.

Если входной сигнал с действительным знаком, величина получившегося БПФ симметрична относительно нулевой частоты (DC). Для БПФ ровной длины, только первое (1 + nfft/2) точки уникальны. Определите количество уникальных значений и сохраните только те уникальные точки.
Возьмите величины в квадрате уникальных значений БПФ. Масштабируйте величины в квадрате (за исключением DC) $2 / (F_{s} N)$ , где N является длиной сигнала до любого нулевого дополнения. Масштабируйте значение DC $1 / (F_{s} N)$ .
Создайте вектор частоты из количества уникальных точек, nfft и частоты дискретизации.
Постройте БПФ получившейся величины в квадрате против частоты.

Эффективность периодограммы

Скрипт Open Live Script

Следующие разделы обсуждают эффективность периодограммы относительно проблем утечки, разрешения, смещения и отклонения.

Спектральная утечка

Рассмотрите PSD конечной длины (длина $L$ ) сигнал $x_{L} (n)$ . Часто полезно интерпретировать $x_{L} (n)$ как результат умножения бесконечного сигнала, $x (n)$ , конечной длиной прямоугольное окно, $w_{R} (n)$ :

$x_{L} (n) = x (n) w_{R} (n) .$

Поскольку умножение во временном интервале соответствует свертке в частотном диапазоне, ожидаемое значение периодограммы в частотном диапазоне

$E {{P_{}^{ˆ}}_{x x} (f)} = \frac{1}{F_{s}} \int_{- F_{s} / 2}^{F_{s} / 2} \frac{\sin^{2} (L π (f - f^{'}) / F_{s})}{L \sin^{2} (π (f - f^{'}) / F_{s})} P_{x x} (f^{'}) d f^{'},$

показ, что ожидаемое значение периодограммы является сверткой истинного PSD с квадратом ядра Дирихле.

Эффект свертки лучше всего изучен для синусоидальных данных. Предположим это $x (n)$ состоит из суммы $M$ комплексные синусоиды:

$x (n) = \sum_{k = 1}^{N} A_{k} e^{j ω_{k} n} .$

Его спектр

$X (ω) = \sum_{k = 1}^{N} A_{k} δ (ω - ω_{k}),$

который для последовательности конечной длины становится

$X (ω) = \int_{- π}^{π} (\sum_{k = 1}^{N} A_{k} δ (ε - ω_{k})) W_{R} (ω - ε) d ε .$

Предыдущее уравнение равно

$X (ω) = \sum_{k = 1}^{N} A_{k} W_{R} (ω - ω_{k}) .$

Таким образом в спектре сигнала конечной длины, дельты Дирака были заменены условиями формы $W_{R} (ω - ω_{k})$ , который соответствует частотной характеристике прямоугольного окна, сосредоточенного на частоте $ω_{k}$ .

Частотная характеристика прямоугольного окна имеет форму периодического sinc:

L = 32;
[h,w] = freqz(rectwin(L)/L,1);
y = diric(w,L);

plot(w/pi,20*log10(abs(h)))
hold on
plot(w/pi,20*log10(abs(y)),'--')
hold off
ylim([-40,0])
legend('Frequency Response','Periodic Sinc')
xlabel('\omega / \pi')

График отображает mainlobe и несколько боковых лепестков, самый большой из которых на приблизительно 13,5 дБ ниже пика mainlobe. Эти лепестки составляют эффект, известный как спектральную утечку. В то время как сигнал бесконечной длины имеет свою силу, сконцентрированную точно в дискретных точках частоты $f_{k}$ , оконное (или усеченный) сигнал имеет континуум степени, "пропущенной" вокруг дискретных точек частоты $f_{k}$ .

Поскольку частотная характеристика короткого прямоугольного окна является намного более плохим приближением к функции дельты Дирака, чем то из более длинного окна, спектральная утечка особенно очевидна, когда записи данных коротки. Рассмотрите следующую последовательность 100 выборок:

fs = 1000;                 % Sampling frequency
t = (0:fs/10)/fs;          % One-tenth second worth of samples
A = [1 2];                 % Sinusoid amplitudes
f = [150;140];             % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

Важно отметить, что эффект спектральной утечки является контингентом только на длине записи данных. Это не последствие того, что периодограмма вычисляется в конечном числе выборок частоты.

Разрешение

Разрешение относится к способности отличить спектральные функции и является ключевой концепцией согласно анализу спектральной эффективности средства оценки.

Для того, чтобы разрешить две синусоиды, которые находятся относительно близко друг к другу в частоте, необходимо для различия между этими двумя частотами быть больше ширины mainlobe пропущенных спектров для любой из этих синусоид. mainlobe ширина задана, чтобы быть шириной mainlobe в точке, где степень является половиной пика mainlobe степень (i.e., 3 дБ шириной). Эта ширина приблизительно равна $f_{s} / L$ .

Другими словами, для двух синусоид частот $f_{1}$ и $f_{2}$ , условие разрешающей способности требует этого

$f_{2} - f_{1} > \frac{F_{s}}{L} .$

В примере выше, где две синусоиды разделяются только на 10 Гц, запись данных должна быть больше 100 выборок, чтобы позволить разрешение двух отличных синусоид периодограммой.

Рассмотрите случай, где этому критерию не соответствуют, что касается последовательности 67 выборок ниже:

fs = 1000;                  % Sampling frequency
t = (0:fs/15)/fs;           % 67 samples
A = [1 2];                  % Sinusoid amplitudes
f = [150;140];              % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

Вышеупомянутая дискуссия о разрешении не рассматривала эффектов шума, поскольку отношение сигнал-шум (SNR) было относительно высоко к настоящему времени. Когда ОСШ является низким, истинные спектральные функции намного более трудно отличить, и шумовые артефакты появляются в спектральных оценках на основе периодограммы. Пример ниже иллюстрирует это:

fs = 1000;                  % Sampling frequency
t = (0:fs/10)/fs;           % One-tenth second worth of samples
A = [1 2];                  % Sinusoid amplitudes
f = [150;140];              % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 2*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

Смещение периодограммы

Периодограмма является смещенной оценкой PSD. Его ожидаемое значение, как ранее показывали, было

$E {{P_{}^{ˆ}}_{x x} (f)} = \frac{1}{F_{s}} \int_{- F_{s} / 2}^{F_{s} / 2} \frac{\sin^{2} (L π (f - f^{'}) / F_{s})}{L \sin^{2} (π (f - f^{'}) / F_{s})} P_{x x} (f^{'}) d f^{'} .$

Периодограмма является асимптотически несмещенной, который очевиден из более раннего наблюдения, что, когда длина записи данных стремится к бесконечности, частотная характеристика прямоугольного окна более тесно аппроксимирует функцию дельты Дирака. Однако в некоторых случаях периодограмма является плохим средством оценки PSD, даже когда запись данных длинна. Это происходит из-за отклонения периодограммы, как объяснено ниже.

Отклонение периодограммы

Отклонение периодограммы, как могут показывать,

$V a r ({P_{}^{ˆ}}_{x x} (f)) = {\begin{array}{c} P_{x x}^{2} (f), & 0 < f < F_{s} / 2, \\ 2 P_{x x}^{2} (f), & f = 0, F_{s} / 2, \end{array}$

который указывает, что отклонение не имеет тенденцию обнулять как длина данных $L$ стремится к бесконечности. В статистических терминах периодограмма не является сопоставимым средством оценки PSD. Тем не менее, периодограмма может быть полезным инструментом для спектральной оценки ситуации, где ОСШ высок, и особенно если запись данных длинна.

Модифицированная периодограмма

Скрипт Open Live Script

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

periodogram позволяет вам вычислять модифицированную периодограмму путем определения окна, которое будет использоваться на данных. Например, сравните прямоугольное окно по умолчанию и Окно Хэмминга. Задайте то же количество точек ДПФ в обоих случаях.

fs = 1000;                  % Sampling frequency
t = (0:fs/10)/fs;           % One-tenth second worth of samples
A = [1 2];                  % Sinusoid amplitudes
f = [150;140];              % Sinusoid frequencies
nfft = 1024;

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),nfft,fs)

periodogram(xn,hamming(length(xn)),nfft,fs)

Можно проверить, что несмотря на то, что боковые лепестки намного менее очевидны в Hamming-оконной периодограмме, два основных peaks более широки. На самом деле 3 дБ шириной из соответствия mainlobe Окну Хэмминга приблизительно дважды больше чем это прямоугольного окна. Следовательно, для фиксированной длины данных, разрешение PSD, достижимое с Окном Хэмминга, является приблизительно половиной настолько достижимого с прямоугольным окном. Конкурирующие интересы mainlobe ширины и высоты бокового лепестка могут быть разрешены в некоторой степени при помощи переменных окон, таких как окно Кайзера.

Непрямоугольная работа с окнами влияет на среднюю степень сигнала, потому что некоторые выборки времени ослабляются, когда умножено на окно. Компенсировать это, periodogram и pwelch нормируйте окно, чтобы иметь среднюю силу единицы. Это гарантирует, что измеренная средняя степень обычно независима от выбора окна. Если частотные составляющие не хорошо разрешены средствами оценки PSD, выбор окна действительно влияет на среднюю степень.

Модифицированная оценка периодограммы PSD

${P_{}^{ˆ}}_{x x} (f) = \frac{{| X (f) |}^{2}}{F_{s} L U},$

где U является постоянной нормализацией окна:

$U = \frac{1}{L} \sum_{n = 0}^{N - 1} {| w (n) |}^{2} .$

Для больших значений LU имеет тенденцию становиться независимым от длины окна. Сложение U когда постоянная нормализация гарантирует, что модифицированная периодограмма является асимптотически несмещенной.

Метод валлийцев

Скрипт Open Live Script

Улучшенное средство оценки PSD является тем, предложенным валлийским языком. Метод состоит из деления данных временных рядов в (возможно перекрывающийся) сегменты, вычисляя модифицированную периодограмму каждого сегмента, и затем составляя в среднем оценки PSD. Результатом является оценка PSD валлийцев. Функция тулбокса pwelch метод валлийцев реализаций.

Усреднение модифицированных периодограмм имеет тенденцию уменьшать отклонение оценки относительно одной оценки периодограммы целой записи данных. Несмотря на то, что перекрытие между сегментами вводит избыточную информацию, этот эффект уменьшен при помощи непрямоугольного окна, которое уменьшает важность или вес, данный в конец выборки сегментов (выборки то перекрытие).

Однако, как упомянуто выше, объединенное использование коротких записей данных и непрямоугольных окон приводит к уменьшаемому разрешению средства оценки. Таким образом, существует компромисс между сокращением отклонения и разрешением. Можно управлять параметрами в методе валлийцев, чтобы получить улучшенные оценки относительно периодограммы, особенно когда ОСШ является низким. Это проиллюстрировано в следующем примере.

Рассмотрите сигнал, состоящий из 301 выборки:

fs = 1000;             % Sampling frequency
t = (0:0.3*fs)/fs;     % 301 samples
A = [2 8];             % Sinusoid amplitudes (row vector)
f = [150;140];         % Sinusoid frequencies (column vector)

xn = A*sin(2*pi*f*t) + 5*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

Мы можем получить спектральную оценку валлийцев для 3 сегментов с 50%-м перекрытием с помощью прямоугольного окна.

pwelch(xn,rectwin(150),50,512,fs)

В периодограмме выше, шум и утечка делают одну из синусоид чрезвычайно неотличимой от искусственного peaks. В отличие от этого несмотря на то, что PSD, произведенный методом валлийцев, имеет более широкий peaks, можно все еще отличить эти две синусоиды, которые выделяются от "уровня шума".

Однако, если мы пытаемся уменьшать отклонение далее, потеря разрешения заставляет одну из синусоид быть потерянной в целом.

pwelch(xn,rectwin(100),75,512,fs)

Смещение и нормализация в методе валлийцев

Метод валлийцев дает к смещенной оценке PSD. Ожидаемое значение оценки PSD:

$E {P_{Welch} (f)} = \frac{1}{F_{s} L U} \int_{- F_{s} / 2}^{F_{s} / 2} {| W (f - f^{'}) |}^{2} P_{x x} (f^{'}) d f^{'},$

где L является длиной сегментов данных, U является той же нормализацией, постоянной существующий в определении модифицированной периодограммы, и W(f) является преобразованием Фурье функции окна. Как имеет место для всех периодограмм, средство оценки валлийцев является асимптотически несмещенным. Для записи данных фиксированной длины смещение оценки валлийцев больше, чем та из периодограммы, потому что длина сегментов меньше продолжительности целой выборки данных.

Отклонение средства оценки валлийцев затрудняет, чтобы вычислить, потому что это зависит и от используемого окна и от суммы перекрытия между сегментами. В основном отклонение обратно пропорционально количеству сегментов, модифицированные периодограммы которых усредняются.

Метод мультизаострения

Скрипт Open Live Script

Периодограмма может быть интерпретирована как фильтрация длины $L$ сигнал, $x_{L} (n)$ , через набор фильтров (ряд фильтров параллельно) $L$ КИХ-полосовые фильтры. Пропускная способность на 3 дБ каждого из этих полосовых фильтров, как могут показывать, приблизительно равна $f_{s} / L$ . Ответ величины каждого из этих полосовых фильтров напоминает ответ прямоугольного окна. Периодограмма может таким образом быть просмотрена как расчет степени каждого отфильтрованного сигнала (i.e., выход каждого полосового фильтра), который использует всего одну выборку каждого отфильтрованного сигнала и принимает что PSD $x_{L} (n)$ является постоянным по пропускной способности каждого полосового фильтра.

Как продолжительность увеличений сигнала, пропускная способность каждого полосового фильтра уменьшения, делая его большим количеством селективного фильтра, и улучшая приближение постоянного PSD по пропускной способности фильтра. Это обеспечивает другую интерпретацию того, почему оценка PSD периодограммы улучшается как продолжительность увеличений сигнала. Однако существует два фактора, очевидные из этой точки зрения, которые ставят под угрозу точность оценки периодограммы. Во-первых, прямоугольное окно дает к плохому полосовому фильтру. Во-вторых, расчет степени при выходе каждого полосового фильтра использует одну выборку выходного сигнала, производя очень грубое приближение.

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

Метод мультизаострения Томпсона (MTM) основывается на этих результатах обеспечить улучшенную оценку PSD. Вместо того, чтобы использовать полосовые фильтры, которые являются чрезвычайно прямоугольными окнами (как в методе периодограммы), метод MTM использует банк оптимальных полосовых фильтров, чтобы вычислить оценку. Эти оптимальные КИХ-фильтры выведены из набора последовательностей, известных как дискретные вытянутые сфероидальные последовательности (DPSSs, также известный как последовательности Slepian).

Кроме того, метод MTM предоставляет параметру пропускной способности времени, с которого можно сбалансировать отклонение и разрешение. Этот параметр дан продуктом пропускной способности времени, $N W$ и это непосредственно связано с количеством заострений, использовался для расчета спектра. Всегда существуют $2 N W - 1$ заострения раньше составляли мнение. Это означает это, как $N W$ увеличения, существует больше оценок спектра мощности и отклонение оценочных уменьшений. Однако пропускная способность каждого заострения также пропорциональна $N W$ , так как $N W$ увеличения, каждая оценка показывает больше спектральной утечки (i.e., более широкий peaks), и полная спектральная оценка более смещается. Для каждого набора данных обычно существует значение для $N W$ это позволяет оптимальный компромисс между смещением и отклонением.

Функцией Signal Processing Toolbox™, которая реализует метод MTM, является pmtm. Используйте pmtm вычислить PSD сигнала.

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
pmtm(xn,4,[],fs)

Путем понижения продукта пропускной способности времени можно увеличить разрешение за счет большего отклонения.

pmtm(xn,1.5,[],fs)

Этот метод является более в вычислительном отношении дорогим, чем метод валлийцев из-за стоимости вычисления дискретных вытянутых сфероидальных последовательностей. Для длинного ряда данных (10 000 точек или больше), полезно вычислить DPSSs однажды и сохранить их в MAT-файле. dpsssave, dpssload, dpssdir, и dpssclear обеспечиваются, чтобы сохранить базу данных сохраненного DPSSs в MAT-файле dpss.mat.

Перекрестная спектральная функция плотности

PSD является особым случаем перекрестной спектральной плотности (CPSD) функция, заданная между двумя сигналами x (n) и y (n) как

$P_{x y} (ω) = \frac{1}{2 π} \sum_{m = - \infty}^{\infty} R_{x y} (m) e^{- j ω m} .$

Как имеет место для корреляции и последовательностей ковариации, тулбокс оценивает PSD и CPSD, потому что длины сигнала конечны.

Оценить перекрестную спектральную плотность двух сигналов равной длины x и y с помощью метода валлийцев, cpsd функционируйте формирует периодограмму как продукт БПФ x и сопряженный из БПФ y. В отличие от PSD с действительным знаком, CPSD является комплексной функцией. cpsd обрабатывает разделение и работу с окнами x и y таким же образом как pwelch функция:

Sxy = cpsd(x,y,nwin,noverlap,nfft,fs)

Оценка передаточной функции

Скрипт Open Live Script

Одно применение метода валлийцев является непараметрической системой идентификации. Примите, что H является линейной, независимой от времени системой, и x (n) и y (n) являются входом к и выходом H, соответственно. Затем спектр мощности x (n) связан с CPSD x (n) и y (n)

$P_{y x} (ω) = H (ω) P_{x x} (ω) .$

Оценка передаточной функции между x (n) и y (n)

$H_{}^{ˆ} (ω) = \frac{{P_{}^{ˆ}}_{y x} (ω)}{{P_{}^{ˆ}}_{x x} (ω)} .$

Этот метод оценивает и величину и информацию о фазе. tfestimate функционируйте использует метод валлийцев, чтобы вычислить CPSD и спектр мощности, и затем формирует их частное для оценки передаточной функции. Используйте tfestimate тот же способ, которым вы используете cpsd функция.

Сгенерируйте сигнал, состоящий из двух синусоид, встроенных в белый Гауссов шум.

rng('default')

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

Отфильтруйте xn сигнала с КИХ-фильтром скользящего среднего значения. Вычислите фактический ответ величины и предполагаемый ответ.

h = ones(1,10)/10;	            % Moving-average filter
yn = filter(h,1,xn);

[HEST,f] = tfestimate(xn,yn,256,128,256,fs);
H = freqz(h,1,f,fs);

Постройте график результатов.

subplot(2,1,1)
plot(f,abs(H))
title('Actual Transfer Function Magnitude')
yl = ylim;
grid
subplot(2,1,2)
plot(f,abs(HEST))
title('Transfer Function Magnitude Estimate')
xlabel('Frequency (Hz)')
ylim(yl)
grid

Функция когерентности

Скрипт Open Live Script

Квадрат когерентности между двумя сигналами x (n) и y (n)

$C_{x y} (ω) = \frac{{| P_{x y} (ω) |}^{2}}{P_{x x} (ω) P_{y y} (ω)} .$

Это частное является вещественным числом между 0 и 1, который измеряет корреляцию между x (n) и y (n) на частоте $ω$ .

mscohere функционируйте берет последовательности xn и yn, вычисляет их спектры мощности и CPSD, и возвращает частное величины, в квадрате из CPSD и продукта спектров мощности. Его опции и операция похожи на cpsd и tfestimate функции.

Сгенерируйте сигнал, состоящий из двух синусоид, встроенных в белый Гауссов шум. Сигнал производится на уровне 1 кГц в течение 1 секунды.

rng('default')

fs = 1000;
t = (0:fs)/fs;
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

Отфильтруйте xn сигнала с КИХ-фильтром скользящего среднего значения. Вычислите и постройте функцию когерентности xn и фильтр выход yn как функция частоты.

h = ones(1,10)/10;
yn = filter(h,1,xn);

mscohere(xn,yn,256,128,256,fs)

Если входная длина последовательности, длина окна и количество наложения точек данных в окне таковы что mscohere работает только с одной записью, функция возвращает все единицы. Это вызвано тем, что функция когерентности для линейно зависимых данных является той.

Документация