Вы хотите различать сигнал без увеличения мощности шума. Функция diff в MATLAB® усиливает шум, и результирующая неточность ухудшается для более высоких производных. Чтобы устранить эту проблему, используйте вместо нее дифференцирующий фильтр.
Проанализируйте смещение перекрытия здания во время землетрясения. Найдите скорость и ускорение как функции времени.
Загрузите файл earthquake. Файл содержит следующие переменные:
drift: Смещение пола, измеряемое в сантиметрах
t: Время, измеряемое в секундах
Fs: Частота дискретизации, равная 1 кГц
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))
tt = t(1:end-delay);
vd = vdrift;
vd(1: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')