Взять производные сигнала

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

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

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

load('earthquake.mat')

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

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

Использовать designfilt для разработки дифференциатора FIR порядка 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

Связанные темы

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