Этот пример показывает, как измерить степень детерминированных периодических сигналов. Несмотря на то, что непрерывный вовремя, периодические детерминированные сигналы производят дискретные спектры мощности.
Мы обеспечиваем два примера того, как измерить среднюю степень сигнала. Примеры используют синусоиды и принимают импеданс загрузки 1 Ома.
В общих сигналах может быть классифицирован в три широких категории, сигналы степени, энергетические сигналы или ни одного. Детерминированные сигналы, которые составлены из синусоид, являются примером сигналов степени, которые имеют бесконечную энергию, но конечную среднюю степень. Случайные сигналы также имеют конечную среднюю силу и попадают в категорию сигналов степени. Переходный сигнал является примером энергетических сигналов, которые запускаются и заканчиваются нулевой амплитудой. Существуют все еще другие сигналы, которые не могут быть охарактеризованы или как степень или как энергетические сигналы.
В нашем первом примере мы оцениваем среднюю степень синусоидального сигнала с пиковой амплитудой 1 вольта и частотной составляющей на уровне 128 Гц.
Fs = 1024; t = 0:1/Fs:1-(1/Fs); A = 1; F1 = 128; x = A*sin(2*pi*t*F1);
Давайте посмотрим на фрагмент сигнала во временном интервале.
idx = 1:128; plot(t(idx),x(idx)) ylabel('Amplitude') xlabel('Time (seconds)') axis tight grid
Теоретическая средняя степень (среднее квадратичное) каждой комплексной синусоиды , который в нашем примере является 0.25 или-6.02 дБ. Так, составление степени в положительных и отрицательных частотах приводит к средней степени .
power_theoretical = (A^2/4)*2
power_theoretical = 0.5000
Мы можем также вычислить в дБ степень, содержавшуюся в положительных частотах:
pow2db(power_theoretical/2)
ans = -6.0206
Чтобы измерить среднюю степень сигнала, мы вызываем periodogram
и задаем опцию 'power'
.
periodogram(x,hamming(length(x)),[],Fs,'centered','power') ylim([-10 -5.5])
Когда мы видим от увеличившего масштаб фрагмент графика, каждая комплексная синусоида имеет среднюю силу примерно-6 дБ.
Другой способ вычислить среднюю степень сигнала путем "интеграции" области под кривой PSD.
periodogram(x,hamming(length(x)),[],Fs,'centered','psd')
Одна вещь заметить в этом графике состоит в том, что peaks графика спектра не имеет той же высоты как тогда, когда мы построили спектр мощности. Причина состоит в том, потому что при взятии власти Спектральной Плотности (PSD) измерения это - область под кривой (который является мерой средней степени), который имеет значение. Мы можем проверить, что путем вызывания функции bandpower
, которая использует прямоугольное приближение, чтобы объединяться под кривой, чтобы вычислить среднюю степень.
[Pxx_hamming,F] = periodogram(x,hamming(length(x)),[],Fs,'psd'); power_freqdomain = bandpower(Pxx_hamming,F,'psd')
power_freqdomain = 0.5000
С тех пор согласно теореме Парсевэла общая средняя степень в синусоиде должна быть равной, вычисляется ли это во временном интервале или частотном диапазоне, мы можем проверить предполагаемую общую среднюю степень нашего сигнала путем подведения сигнала во временном интервале.
power_timedomain = sum(abs(x).^2)/length(x)
power_timedomain = 0.5000
Для второго примера мы оцениваем общую среднюю степень сигнала, содержащего энергию в нескольких частотных составляющих: один в DC, один на уровне 100 Гц, и другом на уровне 200 Гц.
Fs = 1024; t = 0:1/Fs:1-(1/Fs); Ao = 1.5; A1 = 4; A2 = 3; F1 = 100; F2 = 200; x = Ao + A1*sin(2*pi*t*F1) + A2*sin(2*pi*t*F2);
Давайте посмотрим на фрагмент сигнала.
idx = 1:128; plot(t(idx),x(idx)) grid ylabel('Amplitude') xlabel('Time (seconds)')
Как в предыдущем примере, теоретическая средняя степень каждой комплексной синусоиды . Средняя степень DC сигнала равна своей пиковой мощности, поскольку это постоянно и поэтому дано . Составление степени в положительных и отрицательных частотах приводит к общей средней степени сигнала (сумма средней степени каждого гармонического компонента) .
power_theoretical = Ao^2 + (A1^2/4)*2 + (A2^2/4)*2
power_theoretical = 14.7500
Путем вычисления средней степени каждой уникальной частотной составляющей в дБ мы видим, что теоретические результаты совпадают со среднеквадратическим графиком спектра ниже.
pow2db([Ao^2 A1^2/4 A2^2/4])
ans = 1×3
3.5218 6.0206 3.5218
Чтобы измерить среднюю степень сигнала, мы еще раз используем функцию periodogram
, чтобы вычислить и построить спектр мощности сигнала.
periodogram(x,hamming(length(x)),[],Fs,'centered','power') ylim([0 7])
Как в первом примере, мы можем оценить общую среднюю степень сигнала путем "интеграции" под кривой PSD.
periodogram(x,hamming(length(x)),[],Fs,'centered','psd')
Еще раз высота peaks спектрального графика плотности в определенной частотной составляющей не может совпадать с теми графика спектра мощности по причинам, отмеченным в первом примере.
[Pxx, F] = periodogram(x, hamming(length(x)),[],Fs,'centered','psd'); power_freqdomain = bandpower(Pxx,F,'psd')
power_freqdomain = 14.7500
Снова, мы можем проверить предполагаемую среднюю степень сигнала путем вызова теоремы Parseval's и подведения сигнала во временном интервале.
power_timedomain = sum(abs(x).^2)/length(x)
power_timedomain = 14.7500
Вы, возможно, заметили, что, в то время как высота peaks степени и степени спектральные графики плотности отличаются, они отличаются постоянным отношением.
Pxx = periodogram(x,hamming(length(x)),[],Fs,'centered','psd'); Sxx = periodogram(x,hamming(length(x)),[],Fs,'centered','power'); plot(F,Sxx./Pxx) grid axis tight xlabel('Frequency (Hz)') title('Ratio of Power Spectrum to Power Spectral Density')
ratio = mean(Sxx./Pxx)
ratio = 1.3638
Это отношение связано с двухсторонней эквивалентной шумовой пропускной способностью (ENBW) окна. Можно вычислить это отношение непосредственно путем вызова ENBW с окном и его соответствующей частотой дискретизации.
bw = enbw(hamming(length(x)),Fs)
bw = 1.3638
В предыдущих разделах степень была измерена от один или несколько синусоид, имеющих частоту, которая совпала с интервалом. Оценки пиковой мощности обычно менее точны, когда частота сигнала вне интервала. Чтобы видеть этот эффект, создайте синусоиду с количеством нецелого числа циклов за один второй период.
Fs = 1024; t = 0:1/Fs:1-(1/Fs); A = 1; F = 20.4; x = A*sin(2*pi*F*t); nfft = length(x); power_theoretical = pow2db(A^2/4*2);
Создайте Окно Хэмминга и окно с плоской вершиной.
w1 = hamming(length(x)); w2 = flattopwin(length(x));
Вычислите периодограмму x
с помощью Окна Хэмминга. Увеличьте масштаб пика.
h1 = figure; stem(F,power_theoretical,'BaseValue',-50); [Pxx1,f1] = periodogram(x,w1,nfft,Fs,'power'); hold on plot(f1,pow2db(Pxx1),'+-') axis([0 40 -45 0]) legend('Theoretical','Periodogram') xlabel('Frequency (Hz)') ylabel('Power (dB)') title('Periodogram Power Spectrum Estimate') grid
Оценка пиковой мощности ниже теоретического пика, и частота пиковой оценки отличается от истинной частоты.
[Pmax,imax] = max(Pxx1); dPmax_w1 = pow2db(Pmax) - power_theoretical
dPmax_w1 = -1.1046
dFreq = f1(imax) - F
dFreq = -0.4000
Уменьшайте амплитудную ошибку с дополнением нуля
Чтобы видеть, почему это происходит, вычислите периодограмму с помощью большего числа интервалов БПФ.
[Pxx2,f2] = periodogram(x,w1,100*nfft,Fs,'power'); figure stem(F,power_theoretical,'BaseValue',-50) hold on plot(f1,pow2db(Pxx1),'+') plot(f2,pow2db(Pxx2)) hold off axis([0 40 -40 0]) legend('Theoretical Peak','nfft = 1024','nfft = 102400') xlabel('Frequency (Hz)') ylabel('Power (dB)') title('Periodogram Power Spectrum Estimate') grid
В исходной периодограмме спектральный пик расположен между двумя интервалами, и по этой причине предполагаемый пик ниже теоретического пика. Увеличение числа интервалов БПФ дает лучшее изображение спектра, несмотря на то, что это может быть в вычислительном отношении дорогим способом улучшить пиковые измерения.
Уменьшайте амплитудную ошибку с окном с плоской вершиной
Другой способ произвести лучшую оценку для пиковой амплитуды состоит в том, чтобы использовать различное окно. Вычислите периодограмму x
с помощью окна с плоской вершиной.
[Pxx,F1] = periodogram(x,w2,nfft,Fs,'power'); figure(h1) plot(F1,pow2db(Pxx)) legend('Theoretical','Hamming','Flat Top') hold off
Окно с плоской вершиной является широким и плоским. Это производит пиковую оценку ближе для теоретического значения, когда x
не содержит целое число циклов, и следовательно спектральный пик не падает точно на интервал.
dPmax_w2 = pow2db(max(Pxx)) - power_theoretical
dPmax_w2 = -6.2007e-04
Более широкий пик, что продукты окна с плоской вершиной могли быть недостатком при попытке разрешить близко расположенный peaks и частоту измеренного пика, снова отличается от частоты теоретического пика.
Уменьшайте амплитудную ошибку с переприсвоенной периодограммой
Теперь добавьте флаг 'reassigned'
в periodogram
. Переназначение периодограммы использует информацию о фазе, которая обычно отбрасывается, чтобы повторно присвоить сигнал его центру энергии. Процедура может привести к более резким спектральным оценкам. Постройте переприсвоенную периодограмму x
и увеличьте масштаб пика. Используйте Окно Хэмминга и окно с плоской вершиной.
[RPxx1,~,~,Fc1] = periodogram(x,w1,nfft,Fs,'power','reassigned'); [RPxx2,~,~,Fc2] = periodogram(x,w2,nfft,Fs,'power','reassigned'); stem(F,power_theoretical,'*','BaseValue',-40) hold on stem(Fc1,pow2db(RPxx1),'BaseValue',-50) stem(Fc2,pow2db(RPxx2),'BaseValue',-50) hold off legend('Theoretical','Hamming Reassignment','Flattop Reassignment') xlabel('Frequency (Hz)') ylabel('Power (dB)') title('Periodogram Power Spectrum Estimate') axis([19.5 21 -4 -2]) grid
Переприсвоенные оценки степени ближе к теоретическому значению для обоих окон с окном с плоской вершиной, производящим лучшее пиковое измерение.
[RPxx1max,imax1] = max(RPxx1); [RPxx2max,imax2] = max(RPxx2); dPmax_reassign_w1 = pow2db(RPxx1max) - power_theoretical
dPmax_reassign_w1 = -0.0845
dPmax_reassign_w2 = pow2db(RPxx2max) - power_theoretical
dPmax_reassign_w2 = -1.1131e-05
Оценки частоты также улучшены с помощью переприсвоенной периодограммы с окном с плоской вершиной, снова дающим лучшие результаты.
Fc1(imax1)-F
ans = -0.0512
Fc2(imax2)-F
ans = 5.6552e-04