このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

信号の微分係数の取得

ノイズ パワーを増大させずに信号を微分する必要があるとします。MATLAB® 関数 diff を使用するとノイズが増幅されるため、高次導関数において精度がさらに低下します。この問題を解決するには、微分器フィルターを代わりに使用します。

地震発生時の建物の床の変位を解析します。変位の速度と加速度を時間の関数として求めます。

earthquake ファイルを読み込みます。ファイルには、以下の変数が含まれています。

• drift:床の変位 (cm)

• t:時間 (秒)

• Fs:サンプル レート (1 kHz 相当)

load('earthquake.mat')

関数 pwelch を使用して、信号のパワー スペクトルの推定値を表示します。信号エネルギーの大部分が 100 Hz 未満の周波数下にあることに注意してください。

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

designfilt を使用して、50 次の FIR 微分器を設計します。信号エネルギーの大部分が含まれるように、100 Hz の通過帯域周波数と 120 Hz の阻止帯域周波数を指定します。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) = [];

出力には過渡特性も含まれており、この過渡特性の長さはフィルター次数に等しいか群遅延の 2 倍です。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

ドリフト速度を微分して加速度を求めます。ラグは 2 倍の長さになります。遅延を補正するには 2 倍の数のサンプルを破棄し、さらに同数を過渡特性を消去するため破棄します。速度と加速度をプロットします。

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

関連するトピック

• 実践に即したデジタル フィルター処理の紹介