Вы хотите дифференцировать сигнал, не увеличивая шумовую степень. Функциональный diff
MATLAB® усиливает шум, и получившаяся погрешность ухудшается для производных высшего порядка. Чтобы решить эту проблему, используйте фильтр дифференциатора вместо этого.
Анализируйте смещение пола создания во время землетрясения. Найдите скорость и ускорение как функции времени.
Загрузите файл earthquake
. Файл содержит следующие переменные:
drift
: смещение Пола, измеренное в сантиметрах
t:
Время, измеренное в секундах
Fs
: Частота дискретизации, равная 1 кГц
load(fullfile(matlabroot,'examples','signal','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