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

Следующие разделы обсуждают периодограмму, измененную периодограмму, валлийский и методы мультизаострения непараметрической оценки, наряду со связанной функцией 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 степень (т.е. 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{matrix} P_{x x}^{2} (f), & 0 < f < F_{s} / 2, \\ 2 P_{x x}^{2} (f), & f = 0, F_{s} / 2, \end{matrix}$

который указывает, что отклонение не имеет тенденцию обнулять как длина данных $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 ширины и высоты бокового лепестка могут быть разрешены в некоторой степени при помощи переменных окон, таких как окно Kaiser.

Непрямоугольная работа с окнами влияет на среднюю степень сигнала, потому что некоторые выборки времени ослабляются, когда умножено на окно. Чтобы компенсировать это, 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} .$

Для больших значений L U имеет тенденцию становиться независимым от длины окна. Сложение 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_{Валлийский язык} (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$ . Ответ значения каждого из этих полосовых фильтров напоминает ответ прямоугольного окна. Периодограмма может таким образом быть просмотрена как вычисление степени каждого отфильтрованного сигнала (т.е. вывод каждого полосового фильтра), который использует всего одну выборку каждого отфильтрованного сигнала и принимает что PSD $x_{L} (n)$ является постоянным по пропускной способности каждого полосового фильтра.

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

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

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

Кроме того, метод MTM предоставляет параметру пропускной способности времени, с которого можно сбалансировать отклонение и разрешение. Этот параметр дан продуктом пропускной способности времени, $N W$ и это непосредственно связано с количеством заострений, используемых, чтобы вычислить спектр. Всегда существуют $2 N W - 1$ заострения раньше составляли мнение. Это означает это, как $N W$ увеличения, существует больше оценок спектра мощности и отклонение оценочных уменьшений. Однако пропускная способность каждого заострения также пропорциональна $N W$ , так как $N W$ увеличения, каждая оценка показывает больше спектральной утечки (т.е. более широкий 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 и продукта спектров мощности. Его опции и операция подобны функциям tfestimate и cpsd.

Сгенерируйте сигнал, состоящий из двух синусоид, встроенных в белый Гауссов шум. Сигнал выбирается на уровне 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 работает только с одной записью, функция возвращает все единицы. Это вызвано тем, что функция когерентности для линейно зависимых данных является той.

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