信号の微分係数の取得

最新のリリースでは、このページがまだ翻訳されていません。このページの最新版は英語でご覧になれます。

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

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

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

drift: 床の変位 (cm)
t: 時間 (秒)
Fs: サンプルレート (1 kHz 相当)

load(fullfile(matlabroot,'examples','signal','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)')网格次要情节(2,1,2)情节甘氨胆酸(,广告)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

信号の微分係数の取得

参考

関連するトピック

Signal Processing Toolbox ドキュメンテーション

サポート

MATLAB、Simulink、その他の製品をお試しください