Обнаружение искаженного сигнала в шуме

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

Рассмотрим, например, моделируемый выход нелинейного усилителя, который вводит искажение третьего порядка.

Входной сигнал является синусоидой с единичной амплитудой 180 Гц, дискретизированной на частоте 3,6 кГц. Сгенерируйте 10000 выборки.

N = 1e4;
n = 0:N-1;
fs = 3600;
f0 = 180;
t = n/fs;
y = sin(2*pi*f0*t);

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

rng default
noise = randn(size(y));

dispol = [0.5 0.75 1 0];
out = polyval(dispol,y+noise);

ns = 300:500; 

plot(t(ns),[out(ns);polyval(dispol,y(ns))])
xlabel('Time (s)')
ylabel('Signals')
axis tight
legend('With white noise','No white noise')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent With white noise, No white noise.

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

[pxx,f] = pwelch(out,[],[],[],fs);

pwelch(out,[],[],[],fs)

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

Поскольку усилитель вводит искажение третьего порядка, ожидается, что выход сигнал будет иметь:

  • Компонент постоянного тока (нулевая частота);

  • Основной компонент с той же частотой, что и вход, 180 Гц;

  • Две гармоники - частотные составляющие в два и три раза больше частоты входа, 360 и 540 Гц.

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

[pks,lox] = findpeaks(pxx,'NPeaks',4,'SortStr','descend');

hold on
plot(f(lox)/1000,10*log10(pks),'or')
hold off

legend('PSD','Frequency Components')

Figure contains an axes. The axes with title Welch Power Spectral Density Estimate contains 2 objects of type line. These objects represent PSD, Frequency Components.

components = sort([f(lox) f0*(0:3)'])'
components = 2×4

    0.8789  180.1758  360.3516  540.5273
         0  180.0000  360.0000  540.0000

pwelch работает путем деления сигнала на перекрывающиеся сегменты, вычисления периодограммы каждого сегмента и усреднения. По умолчанию функция использует восемь сегментов с 50% перекрытием. Для 10000 выборок это соответствует 2222 выборкам на сегмент.

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

pwelch(out,222,[],[],fs)

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

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

pwelch(out,4444,[],[],fs)

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

См. также

|

Похожие темы