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

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

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

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

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

В качестве первого примера оцените среднюю степень синусоидального сигнала с пиковой амплитудой 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. The axes 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. The axes 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. The axes 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

Теоретическая степень нескольких синусоидов

Для второго примера оцените общую среднюю степень сигнала, содержащего энергию, в нескольких частотных составляющих: один при постоянном токе с амплитудой 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. The axes contains an object of type line.

Как и в предыдущем примере, теоретическая средняя степень каждой комплексной синусоиды является A2/4. Средняя степень сигнала постоянного тока равна его пиковой степени (поскольку она постоянна) и поэтому задается как 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. The axes 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. The axes 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

Снова проверьте предполагаемую среднюю степень сигнала, вызвав теорему Парсеваля и суммируя сигнал во временном интервале.

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. The axes 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. The axes 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. The axes 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. The axes 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. The axes 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

См. также

| | |