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

Вы хотите дифференцировать сигнал, не увеличивая шумовую мощность. Функциональный 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')

Смотрите также

findpeaks | FVTool | designfilt | grpdelay | periodogram

Похожие темы

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