Этот пример показывает, как измерить степень детерминированных периодических сигналов. Несмотря на непрерывность во времени, периодические детерминированные сигналы генерируют дискретные спектры степени. Пример также показывает, как улучшить измерения мощности с помощью метода переназначения.
В целом, сигналы могут быть классифицированы в три широкие категории, сигналы степени, энергетические сигналы или ни один из них. Детерминированные сигналы, которые состоят из синусоидов, являются примером сигналов степени, которые имеют бесконечную энергию, но конечную среднюю степень. Случайные сигналы также имеют конечную среднюю степень и попадают в категорию степеней. Переходный сигнал является примером энергетических сигналов, которые начинают и заканчивают с нулевой амплитудой. Существуют еще другие сигналы, которые не могут быть охарактеризованы как сигналы степени или энергетические сигналы.
В качестве первого примера оцените среднюю степень синусоидального сигнала с пиковой амплитудой 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
Для второго примера оцените общую среднюю степень сигнала, содержащего энергию, в нескольких частотных составляющих: один при постоянном токе с амплитудой 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);
Постройте график первых 128 выборок сигнала.
idx = 1:128; plot(t(idx),x(idx)) grid ylabel('Amplitude') xlabel('Time (seconds)')
Как и в предыдущем примере, теоретическая средняя степень каждой комплексной синусоиды является . Средняя степень сигнала постоянного тока равна его пиковой степени (поскольку она постоянна) и поэтому задается как . Учет степени в положительной и отрицательной частотах приводит к общему среднему значению степени (сумма средней степени каждого гармонического компонента) для сигнала.
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
Снова проверьте предполагаемую среднюю степень сигнала, вызвав теорему Парсеваля и суммируя сигнал во временном интервале.
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
bandpower
| enbw
| periodogram
| pow2db