В этом примере показано, как измерить степень детерминированных периодических сигналов. Несмотря на то, что непрерывный вовремя, периодические детерминированные сигналы производят дискретные спектры мощности. Пример также показывает, как улучшить измерения мощности с помощью метода переназначения.
В общем случае сигналы могут быть классифицированы в три широких категории, сигналы степени, энергетические сигналы или ни одного. Детерминированные сигналы, которые составлены из синусоид, являются примером сигналов степени, которые имеют бесконечную энергию, но конечную среднюю степень. Случайные сигналы также имеют конечную среднюю силу и попадают в категорию сигналов степени. Переходный сигнал является примером энергетических сигналов, который начало и конец с нулевой амплитудой. Существуют все еще другие сигналы, которые не могут быть охарактеризованы или как сигналы степени или как энергетические сигналы.
Как первый пример, оцените среднюю степень синусоидального сигнала с пиковой амплитудой 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 с амплитудой 1.5, один на уровне 100 Гц с амплитудой 4, и один на уровне 200 Гц с амплитудой 3.
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