Измерьте степень детерминированных периодических сигналов

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

Классификация сигнала

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

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

Как первый пример, оцените среднюю степень синусоидального сигнала с пиковой амплитудой 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

Figure contains an axes object. The axes object contains an object of type line.

Теоретическая средняя степень (среднее квадратичное) каждой комплексной синусоиды A2/4, который в этом примере является 0.25 или-6.02 дБ. Так, составление степени в положительных и отрицательных частотах приводит к средней степени 2×A2/4.

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

Figure contains an axes object. The axes object with title Periodogram Power Spectrum Estimate contains an object of type line.

Когда вы видите от увеличившего масштаб фрагмент графика, каждая комплексная синусоида имеет среднюю силу примерно-6 дБ.

Оценка Степени Одной Синусоиды через PSD

Другой способ вычислить среднюю степень сигнала путем "интеграции" области под кривой PSD.

periodogram(x,hamming(length(x)),[],Fs,'centered','psd')

Figure contains an axes object. The axes object with title Power Spectral Density contains an object of type line.

В этом графике 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);

Постройте первые 128 выборок сигнала.

idx = 1:128;
plot(t(idx),x(idx))
grid
ylabel('Amplitude')
xlabel('Time (seconds)')

Figure contains an axes object. The axes object contains an object of type line.

Как в предыдущем примере, теоретическая средняя степень каждой комплексной синусоиды A2/4. Средняя степень DC сигнала равна своей пиковой мощности (поскольку это постоянно), и поэтому дают A02. Составление степени в положительных и отрицательных частотах приводит к общему среднему значению степени (сумма средней степени каждого гармонического компонента) A02+2×A12/4+2×A22/4 для сигнала.

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

Figure contains an axes object. The axes object with title Periodogram Power Spectrum Estimate contains an object of type line.

Оценка степени нескольких синусоид Используя PSD

Как в первом примере, оцените общую среднюю степень сигнала путем "интеграции" под кривой PSD.

periodogram(x,hamming(length(x)),[],Fs,'centered','psd')

Figure contains an axes object. The axes object with title Power Spectral Density contains an object of type line.

Еще раз высота 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

Отношение между спектром мощности, спектральной плотностью мощности и ENBW

Вы, возможно, заметили, что, в то время как высота 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')

Figure contains an axes object. The axes object with title Ratio of Power Spectrum to Power Spectral Density contains an object of type line.

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

Figure contains an axes object. The axes object with title Periodogram Power Spectrum Estimate contains 2 objects of type stem, line. These objects represent Theoretical, Periodogram.

Оценка пиковой мощности ниже теоретического пика, и частота пиковой оценки отличается от истинной частоты.

[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

Figure contains an axes object. The axes object with title Periodogram Power Spectrum Estimate contains 3 objects of type stem, line. These objects represent Theoretical Peak, nfft = 1024, nfft = 102400.

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

Уменьшайте амплитудную ошибку с окном с плоской вершиной

Другой способ произвести лучшую оценку для пиковой амплитуды состоит в том, чтобы использовать различное окно. Вычислите периодограмму x использование окна с плоской вершиной.

[Pxx,F1] = periodogram(x,w2,nfft,Fs,'power');

figure(h1)
plot(F1,pow2db(Pxx))
legend('Theoretical','Hamming','Flat Top')
hold off

Figure contains an axes object. The axes object with title Periodogram Power Spectrum Estimate contains 3 objects of type stem, line. These objects represent Theoretical, Hamming, Flat Top.

Окно с плоской вершиной является широким и плоским. Это производит пиковую оценку ближе для теоретического значения когда 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

Figure contains an axes object. The axes object with title Periodogram Power Spectrum Estimate contains 3 objects of type stem. These objects represent Theoretical, Hamming Reassignment, Flattop Reassignment.

Переприсвоенные оценки степени ближе к теоретическому значению для обоих окон с окном с плоской вершиной, производящим лучшее пиковое измерение.

[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

Смотрите также

| | |