В следующих разделах рассматриваются периодограмма, модифицированная периодограмма, Welch и многопараметрические методы непараметрической оценки, а также связанная функция CPSD, оценка передаточной функции и функция когерентности.
В общих чертах, один из способов оценки PSD процесса состоит в том, чтобы просто найти дискретное преобразование Фурье выборок процесса (обычно выполняемое на сетке с БПФ) и соответствующим образом масштабировать величину в квадрате результата. Эта оценка называется периодограммой.
Оценка периодограммы PSD сигнала ) длины L равна
e-j2πfn/Fs | 2,
где Fs - частота дискретизации.
На практике фактическое вычисление ) может быть выполнено только в конечном количестве частотных точек и обычно использует БПФ. Большинство реализаций метода периодограммы вычисляют N-точечную оценку PSD на частотах
.., N-1.
В некоторых случаях вычисление периодограммы с помощью алгоритма БПФ является более эффективным, если число частот равно двум. Поэтому нередки случаи, когда входной сигнал заполняется нулями для увеличения его длины до степени, равной двум.
В качестве примера периодограммы рассмотрим следующий 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')
Алгоритм
Периодограмма вычисляет и масштабирует выходной сигнал БПФ для получения графика зависимости мощности от частоты следующим образом.
Если входной сигнал имеет действительное значение, то величина результирующего БПФ симметрична относительно нулевой частоты (DC). Для БПФ четной длины - только первый (1 + nfft/ 2) точки уникальны. Определите количество уникальных значений и сохраните только эти уникальные точки.
Возьмем значения в квадрате уникальных значений БПФ. Масштабировать значения в квадрате (за исключением постоянного тока) на ), где N - длина сигнала перед любым заполнением нуля. Масштабируйте значение постоянного тока на FsN).
Создайте частотный вектор из числа уникальных точек, nfft и частоты дискретизации.
Постройте график результирующей величины в квадрате БПФ относительно частоты.
В следующих разделах рассматривается выполнение периодограммы в отношении вопросов утечки, разрешения, смещения и отклонения.
Спектральная утечка
Рассмотрим PSD сигнала конечной длины (длина L) ). Часто полезно интерпретировать n) как результат умножения бесконечного сигнала (n) на прямоугольное конечной длины
wR (n).
Поскольку умножение во временной области соответствует свертке в частотной области, ожидаемое значение периодограммы в частотной области равно
Pxx (f ′) df ′,
показывая, что ожидаемое значение периодограммы является сверткой истинной PSD с квадратом ядра Dirichlet.
Эффект свертки лучше всего понимать для синусоидальных данных. Предположим, что ) состоит из суммы M комплексных синусоид:
=∑k=1NAkejωkn.
Его спектр составляет
λ-λ k),
который для последовательности конечной длины становится
λ-start) d
Предыдущее уравнение равно
λ-λ k).
Так в спектре сигнала конечной длины дельты Дирака были заменены слагаемыми вида ), что соответствует частотной характеристике прямоугольного окна, центрированного на частоте startk.
Частотная характеристика прямоугольного окна имеет форму периодического 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')
На графике показан мейнлобе и несколько сайдлобов, наибольший из которых примерно на 13,5 дБ ниже пика мейнлоба. Эти лепестки объясняют эффект, известный как спектральная утечка. В то время как сигнал бесконечной длины имеет свою мощность, сконцентрированную точно в дискретных частотных точках , оконный (или усеченный) сигнал имеет континуум мощности, «просочившийся» вокруг дискретных частотных точек .
Поскольку частотная характеристика короткого прямоугольного окна намного хуже аппроксимации дельта-функции Дирака, чем у более длинного окна, спектральная утечка особенно очевидна, когда записи данных короткие. Рассмотрим следующую последовательность из 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 дБ ширины). Эта ширина приблизительно равна .
Другими словами, для двух синусоид частот и условие разрешимости требует, чтобы
FSL.
В приведенном выше примере, когда две синусоиды разделены только 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) было относительно высоким до сих пор. Когда 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. Ранее было показано, что его ожидаемое значение
Pxx (f ′) df ′.
Периодограмма асимптотически не смещена, что видно из более раннего наблюдения, что по мере того, как длина записи данных стремится к бесконечности, частотная характеристика прямоугольного окна более близко приближается к дельта-функции Дирака. Однако в некоторых случаях периодограмма является плохой оценкой PSD, даже когда запись данных длинная. Это обусловлено дисперсией периодограммы, как объясняется ниже.
Отклонение периодограммы
Дисперсия периодограммы может быть показана как
), f = 0, Fs/2,
что указывает на то, что дисперсия не стремится к нулю, поскольку длина данных стремится к бесконечности. В статистическом отношении периодограмма не является последовательным оценщиком PSD. Тем не менее, периодограмма может быть полезным инструментом для спектральной оценки в ситуациях, когда SNR является высоким, и особенно, если запись данных является длинной.
Модифицированная периодограмма оканчивает сигнал временной области перед вычислением DFT для сглаживания краев сигнала. Это приводит к уменьшению высоты боковых обломков или спектральной утечки. Это явление порождает интерпретацию боковых лопастей как ложных частот, вводимых в сигнал резким усечением, которое происходит при использовании прямоугольного окна. Для непрямоугольных окон конечные точки усеченного сигнала ослабляются плавно, и, следовательно, вводимые ложные частоты гораздо менее серьезны. С другой стороны, непрямоугольные окна также расширяют основной блок, что приводит к уменьшению разрешения.
periodogram позволяет вычислить измененную периодограмму, указав окно, которое будет использоваться для данных. Например, сравните прямоугольное окно по умолчанию и окно Хэмминга. Укажите одинаковое количество точек DFT в обоих случаях.
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)
Вы можете проверить, что хотя боковые лопасти гораздо менее очевидны в периодограмме Хэмминга, два основных пика шире. На самом деле ширина на 3 дБ соответствия mainlobe окну Хэмминга приблизительно дважды больше чем это прямоугольного окна. Следовательно, для фиксированной длины данных разрешение PSD, достижимое с окном Хэмминга, составляет приблизительно половину разрешения, достижимого с прямоугольным окном. Конкурирующие интересы ширины основного блока и высоты бокового блока могут быть в некоторой степени разрешены с помощью переменных окон, таких как окно Кайзера.
Непрямоугольное оконное отображение влияет на среднюю мощность сигнала, поскольку некоторые отсчеты времени ослабляются при умножении на окно. Чтобы компенсировать это, periodogram и pwelch нормализовать окно, чтобы оно имело среднюю степень единства. Это гарантирует, что измеренная средняя мощность обычно не зависит от выбора окна. Если частотные компоненты плохо разрешены оценщиками PSD, выбор окна влияет на среднюю мощность.
Модифицированная оценка периодограммы ИПУ
2FSLU,
где U - константа нормализации окна:
2.
Для больших значений L, U стремится стать независимым от длины окна. Добавление U поскольку константа нормализации обеспечивает асимптотически несмещенную измененную периодограмму.
Усовершенствованная оценка PSD предложена Welch. Способ состоит в разделении данных временных рядов на (возможно, перекрывающиеся) сегменты, вычислении модифицированной периодограммы каждого сегмента, а затем усреднении оценок PSD. Результатом является оценка PSD Уэлча. Функция панели инструментов pwelch реализует метод Уэлча.
Усреднение модифицированных периодограмм имеет тенденцию уменьшать дисперсию оценки относительно одной оценки периодограммы всей записи данных. Хотя перекрытие между сегментами вводит избыточную информацию, этот эффект уменьшается за счет использования непрямоугольного окна, что уменьшает важность или вес, придаваемый конечным выборкам сегментов (выборкам, которые перекрываются).
Однако, как упоминалось выше, комбинированное использование коротких записей данных и непрямоугольных окон приводит к снижению разрешения оценщика. Таким образом, существует компромисс между уменьшением дисперсии и разрешением. Можно манипулировать параметрами в методе Уэлча для получения улучшенных оценок относительно периодограммы, особенно когда SNR низкий. Это проиллюстрировано в следующем примере.
Рассмотрим сигнал, состоящий из 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)
На приведенной выше периодограмме шум и утечка делают одну из синусоид по существу неотличимой от искусственных пиков. В противоположность этому, хотя выпускаемая методом Уэлча PSD имеет более широкие пики, всё же можно выделить две синусоиды, которые выделяются из «шумового пола».
Однако, если мы попытаемся еще больше уменьшить дисперсию, потеря разрешения приводит к полной потере одной из синусоид.
pwelch(xn,rectwin(100),75,512,fs)
Метод Уэлча дает смещенную оценку PSD. Ожидаемое значение оценки PSD:
(f ′) df ′,
где L - длина сегментов данных, U - та же самая константа нормализации, присутствующая в определении модифицированной периодограммы, и W (f) - преобразование Фурье оконной функции. Как и для всех периодограмм, оценщик Уэлча асимптотически не смещен. Для записи данных фиксированной длины смещение оценки Уэлча больше, чем у периодограммы, потому что длина сегментов меньше, чем длина всей выборки данных.
Дисперсию оценщика Уэлча трудно вычислить, поскольку она зависит как от используемого окна, так и от величины перекрытия между сегментами. В основном, дисперсия обратно пропорциональна количеству сегментов, модифицированные периодограммы которых усредняются.
periodogram может интерпретироваться как фильтрация длины сигнал, (n), через банк фильтра (ряд фильтров параллельно) полосовых фильтров ЕЛИ L. Ширина полосы пропускания 3 дБ каждого из этих полосовых фильтров может быть приблизительно равна fs/L. Амплитудный отклик каждого из этих полосовых фильтров напоминает отклик прямоугольного окна. Таким образом, периодограмма может рассматриваться как вычисление мощности каждого отфильтрованного сигнала (то есть выходного сигнала каждого полосового фильтра), которое использует только одну выборку каждого отфильтрованного сигнала и предполагает, что PSD n) является постоянной по ширине полосы пропускания каждого полосового фильтра.
По мере увеличения длины сигнала полоса пропускания каждого полосового фильтра уменьшается, делая его более избирательным фильтром, и улучшая аппроксимацию постоянной PSD над полосой пропускания фильтра. Это обеспечивает другую интерпретацию того, почему оценка PSD периодограммы улучшается по мере увеличения длины сигнала. Однако есть два фактора, очевидных с этой точки зрения, которые снижают точность оценки периодограммы. Во-первых, прямоугольное окно дает плохой полосовой фильтр. Во-вторых, вычисление мощности на выходе каждого полосового фильтра основывается на одной выборке выходного сигнала, создавая очень грубую аппроксимацию.
Методу Уэлча может быть дана аналогичная интерпретация в терминах банка фильтров. В реализации Уэлча несколько выборок используются для вычисления выходной мощности, что приводит к уменьшению дисперсии оценки. С другой стороны, ширина полосы пропускания каждого полосового фильтра больше, чем ширина полосы пропускания, соответствующая способу периодограммы, что приводит к потере разрешения. Таким образом, модель банка фильтров обеспечивает новую интерпретацию компромисса между дисперсией и разрешением.
Многоаппаратный метод Томпсона (MTM) основан на этих результатах, чтобы обеспечить улучшенную оценку PSD. Вместо использования полосовых фильтров, которые являются по существу прямоугольными окнами (как в методе периодограммы), метод MTM использует набор оптимальных полосовых фильтров для вычисления оценки. Эти оптимальные FIR-фильтры получены из набора последовательностей, известных как дискретные простейшие сфероидальные последовательности (DPSS, также известные как Slepian последовательности).
Кроме того, метод MTM обеспечивает параметр временной полосы пропускания, с помощью которого можно сбалансировать дисперсию и разрешение. Этот параметр задается продуктом временной полосы пропускания и непосредственно связан с количеством конусов, используемых для вычисления спектра. Для формирования оценки всегда используются 2NW-1 конусы. Это означает, что по мере увеличения NW появляется больше оценок спектра мощности, и дисперсия оценки уменьшается. Однако ширина полосы пропускания каждого сужения также пропорциональна , так что по мере увеличения NW каждая оценка демонстрирует большую спектральную утечку (т.е. более широкие пики), и общая спектральная оценка является более смещенной. Для каждого набора данных обычно существует значение , которое позволяет оптимальный компромисс между смещением и дисперсией.
Функция 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) как
− jstartm.
Как и в случае корреляционных и ковариационных последовательностей, инструментарий оценивает 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) посредством
Pxx (λ).
Оценка передаточной функции между x (n) и y (n) равна
Pˆxx (ω).
Этот метод оценивает как величину, так и информацию о фазе. 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) равна
) Pyy (λ).
Это частное является действительным числом между 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 работает только с одной записью, функция возвращает все записи. Это происходит потому, что функция когерентности для линейно зависимых данных равна единице.
cpsd | mscohere | periodogram | pmtm | pwelch | tfestimate