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

Periodogram

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

Оценка периодограммы PSD сигнала $$x_L(n)$$ длины L равна

где Fs - частота дискретизации.

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

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

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

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')

Алгоритм

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

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

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

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

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

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

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

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

Его спектр

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

Предшествующее уравнение равно

Так в спектре сигнала конечной длины дельты Дирака были заменены терминами вида $$W_R(\omega -{\omega }_k)$$, что соответствует частотной характеристике прямоугольного окна с центром на частоте $${\omega }_k$$.

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

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')

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

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

Разрешение

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

В порядок разрешения двух синусоидов, которые относительно близки друг к другу по частоте, необходимо, чтобы различие между этими двумя частотами было больше, чем ширина основной области утечек спектров для любой из этих синусоид. Ширина майнлоба определяется как ширина мэнлоба в точке, где степень составляет половину пиковой степени мэнлоба (т.е. ширина 3 дБ). Эта ширина приблизительно равна $$f_s/L$$.

Другими словами, для двух синусоидов частот $$f_1$$ и $$f_2$$, условие разрешимости требует, чтобы

В приведенном выше примере, где две синусоиды разделены только 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)

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

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. Его ожидаемое значение ранее было показано как

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

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

Можно показать, что отклонение периодограммы

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

Измененная пародограмма

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

The 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)

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

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

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

где U - константа нормализации окна:

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

Метод Уэлча

Улучшенная оценка PSD предложена Welch. Метод состоит в том, чтобы разделить данные временных рядов на (возможно, перекрывающиеся) сегменты, вычислить модифицированную периодограмму каждого сегмента, а затем усреднить оценки PSD. Результатом является оценка PSD компании Welch. Функция тулбокса 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:

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

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

Метод мультитапера

Периодограмма может быть интерпретирована как фильтрация длины $$L$$ сигнал, $$x_L(n)$$, через банк фильтров (набор фильтров параллельно) $$L$$ Полосовые фильтры конечные импульсные характеристики. Ширина полосы 3 дБ каждого из этих полосно-пропускающих фильтров может быть показана приблизительно равной $$f_s/L$$. Амплитуда каждого из этих полосно-пропускающих фильтров напоминает величину прямоугольного окна. Таким образом, периодограмма может рассматриваться как расчет степени каждого фильтрованного сигнала (то есть выхода каждого полосно-пропускающего фильтра), которое использует только одну выборку каждого фильтрованного сигнала и принимает, что PSD $$x_L(n)$$ является постоянным по ширине полосы пропускания каждого полосно-пропускающего фильтра.

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

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

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

В сложение метод MTM предоставляет параметр полосы времени, с которым можно сбалансировать отклонение и разрешение. Этот параметр задается продуктом time-bandwidth, $$NW$$ и это непосредственно связано с количеством сужений, используемых для вычисления спектра. Всегда есть $$2NW-1$$ tapers используется для формирования оценки. Это означает, что, как $$NW$$ увеличивается, существует больше оценок спектра степени, и отклонение оценки уменьшается. Однако шумовая полоса каждой конусности также пропорциональна $$NW$$, так как $$NW$$ увеличивается, каждая оценка показывает больше спектральных утечек (то есть более широкого peaks) и общая спектральная оценка является более смещенной. Для каждого набора данных обычно существует значение для $$NW$$ что позволяет добиться оптимального компромисса между смещением и отклонением.

Функция 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 точек или более) полезно один раз вычислить DPSS и сохранить их в MAT-файле. dpsssave, dpssload, dpssdir, и dpssclear предусмотрены для хранения базы данных сохраненных DPSS в MAT-файле dpss.mat.

Функция кросс-спектральной плотности

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

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

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

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

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

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

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

Этот метод оценивает и величину, и информацию фазы. The 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

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

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

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

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