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

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

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

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

Теоретическая средняя мощность (среднеквадратичная) каждой комплексной синусоиды равна ${\mathit{}}^{\mathrm{}}\mathrm{A2/4}$, что в данном примере равно 0,25 или -6,02 дБ. Так, учёт мощности в положительной и отрицательной частотах приводит к средней мощности $2{\mathit{\times }}^{\mathrm{}}\mathrm{}$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')

На этом графике пики графика спектра не имеют такой же высоты, как на графике спектра мощности. Высоты различны, потому что именно площадь под кривой - которая является мерой средней мощности - имеет значение при проведении измерений спектральной плотности мощности (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)')

Как и в предыдущем примере, теоретическая средняя мощность каждой комплексной синусоиды равна ${\mathit{}}^{\mathrm{}}\mathrm{A2/4}$. Средняя мощность постоянного тока сигнала равна его пиковой мощности (так как она постоянна) и поэтому задаётся ${\mathit{}}_{\mathrm{}}^{\mathrm{A02}}$. Учет мощности в положительной и отрицательной частотах приводит к общему среднему значению мощности (сумме средней мощности каждой гармонической составляющей) ${\mathit{}}_{\mathrm{}}^{\mathrm{A02}}\mathrm{+}{\mathit{2}}_{\mathrm{}}^{\mathrm{\times }}\mathrm{}\mathrm{}{\mathit{A12/4}}_{\mathrm{}}^{\mathrm{+}}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')

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

[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

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

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

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

Уменьшение амплитудной ошибки с помощью переназначенной периодограммы

Теперь добавьте '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