Берите производные сигнала

Вы хотите дифференцировать сигнал, не увеличивая степень шума. Функциональные diff MATLAB ® усиливает шум, и результирующая неточность ухудшается для более высоких производных. Чтобы устранить эту проблему, используйте вместо этого дифференцирующий фильтр.

Анализируйте смещение пола создания во время землетрясения. Найдите скорость и ускорение как функции времени.

Загрузите файл earthquake. Файл содержит следующие переменные:

load('earthquake.mat')

Использование pwelch для отображения оценки спектра степени сигнала. Обратите внимание, что большая часть энергии сигнала содержится в частотах ниже 100 Гц.

pwelch(drift,[],[],[],Fs)

Использование designfilt для разработки конечной импульсной характеристики дифференциатора порядка 50. Для включения большей части энергии сигнала задайте частоту полосы пропускания 100 Гц и частоту полосы остановки 120 Гц. Осмотрите фильтр с помощью fvtool.

Nf = 50; 
Fpass = 100; 
Fstop = 120;

d = designfilt('differentiatorfir','FilterOrder',Nf, ...
    'PassbandFrequency',Fpass,'StopbandFrequency',Fstop, ...
    'SampleRate',Fs);

fvtool(d,'MagnitudeDisplay','zero-phase','Fs',Fs)

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

dt = t(2)-t(1);

vdrift = filter(d,drift)/dt;

Отфильтрованный сигнал задерживается. Использование grpdelay чтобы определить, что задержка составляет половину порядка фильтра. Компенсируйте его отбросом выборок.

delay = mean(grpdelay(d))

delay = 25

tt = t(1:end-delay);
vd = vdrift;
vd(1:delay) = [];

Выходы также содержат переходный процесс, длина которого равна порядку фильтра или вдвое больше задержки группы. delay выборки отбрасывали выше. Сбросьте delay больше для устранения переходного процесса.

tt(1:delay) = [];
vd(1:delay) = [];

Постройте график дрейфа и скорости дрейфа. Использование findpeaks чтобы проверить, что максимумы и минимумы дрейфа соответствуют нулевым пересечениям его производной.

[pkp,lcp] = findpeaks(drift);
zcp = zeros(size(lcp));

[pkm,lcm] = findpeaks(-drift);
zcm = zeros(size(lcm));

subplot(2,1,1)
plot(t,drift,t([lcp lcm]),[pkp -pkm],'or')
xlabel('Time (s)')
ylabel('Displacement (cm)')
grid

subplot(2,1,2)
plot(tt,vd,t([lcp lcm]),[zcp zcm],'or')
xlabel('Time (s)')
ylabel('Speed (cm/s)')
grid

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

adrift = filter(d,vdrift)/dt;

at = t(1:end-2*delay);
ad = adrift;
ad(1:2*delay) = [];

at(1:2*delay) = [];
ad(1:2*delay) = [];

subplot(2,1,1)
plot(tt,vd)
xlabel('Time (s)')
ylabel('Speed (cm/s)')
grid

subplot(2,1,2)
plot(at,ad)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid

Вычислите ускорение с помощью diff. Добавьте нули, чтобы компенсировать изменение размера массива. Сравните результат с полученным с фильтром. Заметьте количество высокочастотного шума.

vdiff = diff([drift;0])/dt;
adiff = diff([vdiff;0])/dt;

subplot(2,1,1)
plot(at,ad)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid
legend('Filter')
title('Acceleration with Differentiation Filter')

subplot(2,1,2)
plot(t,adiff)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid
legend('diff')

См. также

designfilt | findpeaks | FVTool | grpdelay | periodogram

Похожие темы

• Практическое введение в цифровую фильтрацию