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

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

Мы обеспечиваем два примера того, как измерить среднюю степень сигнала. Примеры используют синусоиды и принимают импеданс загрузки 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

Теоретическая средняя степень (среднее квадратичное) каждой комплексной синусоиды 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])

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

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

Другой способ вычислить среднюю степень сигнала путем "интеграции" области под кривой 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)')

Как в предыдущем примере, теоретическая средняя степень каждой комплексной синусоиды 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])

Оценка Степени Нескольких Синусоид через PSD

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

Отношение между спектром мощности, степень спектральная плотность и 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')

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