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

envspectrum

機械診断の包絡線スペクトル

構文

es = envspectrum(x,fs)

es = envspectrum(xt)

es = envspectrum(___,Name=Value)

[es,f,env,t] = envspectrum(___)

envspectrum(___)

説明

es = envspectrum(x,fs) は、fs のレートでサンプリングされた信号 x の包絡線スペクトルを返します。x が行列の場合、この関数は各列の包絡線スペクトルを個別に計算し、結果を対応する es の列に返します。

es = envspectrum(xt) は、MATLAB^® timetable xt に保存された信号の包絡線スペクトルを返します。

es = envspectrum(___,Name=Value) は、名前と値の引数を使用して、上記の任意の構文に追加オプションを指定します。オプションには、包絡線信号の計算に使用するアルゴリズムやスペクトルを推定する周波数帯域があります。

[es,f,env,t] = envspectrum(___) は、es が計算された周波数のベクトル f を返します。env は包絡線信号です。t は env が計算された時間です。

出力引数を設定せずに envspectrum(___) を使用すると、現在の Figure に包絡線信号と包絡線スペクトルがプロットされます。

例

振動信号の包絡線スペクトル

2 つの振動信号、正常なベアリングの信号と破損したベアリングの信号をシミュレーションします。これらの包絡線スペクトルを計算して比較します。

ピッチの直径が 12 cm のベアリングは 8 つの回転要素を持ちます。各回転要素の直径は 2 cm です。内輪が 1 秒あたり 25 回駆動される間、外輪は静止状態を保ちます。加速度計はベアリングの振動を 10 kHz でサンプリングします。

fs = 10000;
f0 = 25;
n = 8;
d = 0.02;
p = 0.12;

正常なベアリングの振動信号には、駆動点周波数の次数が複数含まれます。0.1 秒分のデータをプロットします。

t = 0:1/fs:1-1/fs;
z = [1 0.5 0.2 0.1 0.05]*sin(2*pi*f0*[1 2 3 4 5]'.*t)/5;

plot(t,z)
xlim([0.4 0.5])

ベアリングの外輪に欠陥があると、ベアリングで 5 ミリ秒の一連の影響を引き起こします。いずれ、これらの影響はベアリングの摩耗につながります。影響は、ベアリングの外輪転動体通過周波数 (BPFO: Ball Pass Frequency Outer race) で発生します。

$$BPFO=\frac{1}{2}n{f}_0[1-\frac{d}{p}\cos \theta ],$$

ここで、$$f_0$$ は駆動レート、n は回転要素の数、d は回転要素の直径、p はベアリングのピッチの直径、θ はベアリングの接触角です。接触角は 0 と仮定して BPFO を計算します。

ca = 0;
bpfo = n*f0/2*(1-d/p*cos(ca))

bpfo = 83.3333

フラット トップ ウィンドウによってウィンドウを適用した 3 kHz の正弦波として、各影響をモデル化します。櫛形関数との畳み込みにより、影響を周期的にします。0.1 秒分のデータをプロットします。

fImpact = 3000;
tImpact = 0:1/fs:5e-3-1/fs;
xImpact = sin(2*pi*fImpact*tImpact).*flattopwin(length(tImpact))'/10;

xComb = zeros(size(t));
xComb(1:fs/bpfo:end) = 1;

x = conv(xComb,xImpact,'same')/3;

plot(t,x+z)
xlim([0.4 0.5])

ホワイト ガウス ノイズを信号に付加します。1/30² のノイズ分散を指定します。0.1 秒分のデータをプロットします。

yGood = z + randn(size(z))/30;
yBad = x+z + randn(size(z))/30;
plot(t,yGood,t,yBad)
xlim([0.4 0.5])
legend('Healthy','Damaged')

包絡線信号とスペクトルを計算してプロットします。

envspectrum([yGood' yBad'],fs)
xlim([0 10*bpfo]/1000)

ピークの位置と BPFO の高調波の周波数を比較します。包絡線スペクトルの BPFO 高調波は、ベアリング摩耗の兆候です。

harmImpact = (1:10)*bpfo;
[X,Y] = meshgrid(harmImpact,ylim);

hold on
plot(X/1000,Y,':k')
legend('Healthy','Damaged','BPFO harmonics')
hold off

信号のウェルチ スペクトルを計算します。5 Hz の周波数分解能を指定します。

figure
pspectrum([yGood' yBad'],fs,'FrequencyResolution',5)
legend('Healthy','Damaged')

スペクトルの下端では、駆動周波数とその次数により、その他の特徴が不明確になります。異常がないベアリングのスペクトルと破損したベアリングのスペクトルは、区別できません。

xlim([0 10*bpfo]/1000)

故障したベアリングのスペクトルは、影響を及ぼす周波数によって変調された BPFO 高調波を示します。

xlim((bpfo*[-10 10]+fImpact)/1000)

timetable の包絡線スペクトル

10 ミリ秒ごとに 1 回転するベアリングの振動信号に似た 2 チャネル信号を生成します。信号は 10 kHz で 0.2 秒間サンプリングされます。これはベアリングの 20 回転に相当します。

fs = 10000;
tmax = 20;
mlt = 0.01;
t = 0:1/fs:mlt-1/fs;

各 10 ミリ秒間隔の間に:

• 最初のチャネルは減衰定数が 700、正弦波周波数が 600 Hz の減衰正弦波です。

• 2 番目のチャネルは減衰定数が 800、正弦波周波数が 500 Hz の別の減衰正弦波です。2 番目のチャネルは最初のチャネルから 5 ミリ秒遅れています。

信号をプロットします。

y1 = sin(2*pi*600*t).*exp(-700*t);
y2 = sin(2*pi*500*t).*exp(-800*t);
y2 = [y2(51:100) y2(1:50)];

T = (0:1/fs:mlt*tmax-1/fs)';
Y = repmat([y1;y2],1,tmax)';

plot(T,Y)

時間間隔 T を使用して duration 配列を作成します。この duration 配列と 2 チャネル信号で timetable を構築します。

dt = seconds(T);
ttb = timetable(dt,Y);

出力引数を設定せずに envspectrum を使用して、2 つのチャネルの包絡線信号と包絡線スペクトルを表示します。最後の 100 Hz 間隔を除くナイキスト区間全体のスペクトルを計算します。

envspectrum(ttb,'Band',[100 4900])

信号の包絡線スペクトルには、反復レート 1/0.01 = 0.1 kHz の整数倍の位置にピークがあります。これは予測とまったく同じです。envspectrum は、高周波数の正弦波成分を除去し、低周波数の反復動作に注目します。包絡線スペクトルが回転機の解析に役立つツールであるのは、このためです。

包絡線信号とその計算回数を計算します。出力変数の型を確認します。

[~,~,ttbenv,ttbt] = envspectrum(ttb,'Band',[100 4900]);
whos ttb*

  Name           Size            Bytes  Class        Attributes

  ttb         2000x1             48977  timetable              
  ttbenv      2000x1             48985  timetable              
  ttbt        2000x1             16002  duration               

時間ベクトルの型は、入力 timetable の時間値のような、duration です。出力 timetable のサイズは、入力 timetable と同じです。

入力 timetable の各チャネルを、別々の変数として保存します。包絡線信号と時間ベクトルを計算します。出力の型を確認します。

btb = timetable(dt,Y(:,1),Y(:,2));

[~,~,btbenv,btbt] = envspectrum(btb,'Band',[100 4900]);
whos btb*

  Name           Size            Bytes  Class        Attributes

  btb         2000x2             49199  timetable              
  btbenv      2000x2             49219  timetable              
  btbt        2000x1             16002  duration               

出力 timetable のサイズは、入力 timetable と同じです。

変調されたパルスの包絡線スペクトル

1 kHz で 5 秒間サンプリングされた信号を生成します。信号は、T = 0.25 秒ごとに繰り返される 0.01 秒の矩形パルスで構成されています。搬送周波数 150 Hz の正弦波で信号の振幅変調を行います。

fs = 1e3;
tmax = 5;

t = 0:1/fs:tmax;
y = pulstran(t,0:0.25:tmax,rectpuls=0.01);

fc = 150;
z = modulate(y,fc,fs);

元の信号と変調した信号をプロットします。最初の数サイクルのみを表示します。

plot(t,y,t,z,"-")
grid on
axis([0 1 -1.1 1.1])

信号の包絡線と包絡線スペクトルを計算します。複素数復調を使用して信号の包絡線を決定します。搬送周波数を中心として 20 Hz 間隔で包絡線スペクトルを計算します。

[q,f,e,te] = envspectrum(z,fs,Method="demod",Band=[fc-10 fc+10]);

包絡線信号と包絡線スペクトルをプロットします。0 から 50 Hz の区間を拡大します。

subplot(2,1,1)
plot(te,e)
xlabel("Time")
title("Envelope")

subplot(2,1,2)
plot(f,q)
xlim([0 50])
xlabel("Frequency")
title("Envelope Spectrum")

包絡線信号の周期 T = 0.25 秒は、元の信号と同じです。包絡線スペクトルのパルスは 1 / T = 4 Hz です。

計算を繰り返しますが、今度は関数 hilbert を使用して包絡線を計算します。10 次有限インパルス応答 (FIR) フィルターを使用して信号をバンドパス フィルター処理します。組み込み関数 envspectrum を使用して包絡線信号と包絡線スペクトルをプロットします。

envspectrum(z,fs,Method="hilbert",FilterOrder=10)

分散 1/3 のホワイト ガウス ノイズに信号を組み込みます。結果をプロットします。

zn = z + randn(size(z))/3;

plot(t,zn,"-")
grid on
axis([0 1 -1.1 1.1])

包絡線信号と包絡線スペクトルを計算して表示します。搬送周波数を中心として 10 Hz 間隔で、複素数復調を使用して包絡線スペクトルを計算します。0 から 50 Hz の区間を拡大します。

envspectrum(zn,fs,Band=[fc-5 fc+5])
xlim([0 50])

入力引数

x — 入力信号

ベクトル | 行列

ベクトルまたは行列として指定される入力信号。x がベクトルの場合、単一チャネルとして取り扱われます。x が行列の場合、envspectrum は各列の包絡線スペクトルを個別に計算し、結果を対応する es の列に返します。

例: cos(pi/4*(0:159))+randn(1,160) は単一チャネルの行ベクトル信号です。

例: cos(pi./[4;2]*(0:159))'+randn(160,2) は 2 チャネル信号です。

データ型: single | double
複素数のサポート: あり

fs — サンプル レート

正の実数のスカラー

サンプル レート。正の実数スカラーとして指定します。

データ型: single | double

xt — 入力 timetable

timetable

入力 timetable。xt に増加する有限の行時間を含めなければなりません。xt がマルチチャネル信号を表す場合は、行列を含む単一変数か、ベクトルで構成される複数の変数のどちらかをもたなければなりません。

timetable が欠損している場合や時間点が重複している場合、欠損または重複する時間および非等間隔の時間をもつ timetable の整理のヒントを使用して修正できます。

例: timetable(seconds(0:4)',randn(5,2)) は 1 Hz で 4 秒間サンプリングされた 2 チャネルの確率変数を指定します。

データ型: single | double

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、コンマを使用して名前と値をそれぞれ区切り、Name を引用符で囲みます。

例: 'Method','hilbert','FilterOrder',30,'Band',[0 fs/4] は、30 次バンドパス フィルターを使用し解析信号の包絡線を計算して、0 からナイキスト周波数の半分までの包絡線スペクトルを計算します。

Method — 包絡線信号の計算アルゴリズム

"demod" (既定値) | "hilbert"

包絡線信号の計算アルゴリズム。"hilbert" または "demod" として指定します。詳細については、アルゴリズムを参照してください。

Band — 包絡線スペクトルを計算する周波数帯域

[fs/4 fs*3/8] (既定値) | 2 要素ベクトル

包絡線スペクトルを計算する周波数帯域。0 とナイキスト周波数の間で正確に増加する値から成る 2 要素ベクトルとして指定します。

データ型: single | double
複素数のサポート: あり

FilterOrder — FIR フィルター次数

50 (既定値) | 正の整数スカラー

FIR フィルター次数。正の整数スカラーとして指定します。

• Method が "hilbert" の場合、この引数は FIR バンドパス フィルターの次数を指定します。

• Method が "demod" の場合、この引数は FIR ローパス フィルターの次数を指定します。

データ型: single | double

出力引数

es — 包絡線スペクトル

ベクトル | 行列

包絡線スペクトル。ベクトルまたは行列として返されます。

f — 周波数

ベクトル

包絡線スペクトルが計算される周波数。ベクトルとして返されます。

env — 包絡線信号

ベクトル | 行列 | timetable

包絡線信号。ベクトル、行列、または timetable として返されます。

envspectrum への入力が timetable である場合、env も timetable になります。env の時間値の形式は、入力 timetable の時間値と同じです。

• 入力が行列を含む単一変数をもつ timetable である場合、env は行列を含む単一変数をもちます。

• 入力がベクトルで構成される複数の変数をもつ timetable である場合、env はベクトルで構成される複数の変数をもちます。

t — 時間値

ベクトル

包絡線信号が計算される時間値。ベクトルとして返されます。

envspectrum への入力が timetable である場合、t の形式は、入力 timetable の時間値と同じです。

アルゴリズム

envspectrum は、最初に入力信号 x から DC バイアスを取り除き、包絡線信号を計算します。

• Method を "hilbert" に設定した場合、関数は次の処理を行います。

  1. 信号をバンドパス フィルター処理します。FIR フィルターの次数は FilterOrder と、ba(1) および ba(2) のカットオフ周波数で指定されます。ここで ba は Band を使用して指定される周波数帯域です。

  2. 関数 hilbert を使用して解析信号を計算します。

  3. 包絡線信号を解析信号の絶対値として計算します。

• Method を "demod" に設定した場合、関数は次の処理を行います。

  1. 信号の複素数復調を実行します。信号は、exp(j2πf₀t) で乗算されます。ここで、f₀ = (ba(1) + ba(2))/2 です。

  2. 復調した信号をローパス フィルター処理し、解析信号を計算します。FIR フィルターの次数は、FilterOrder と、(ba(2) – ba(1))/2 のカットオフ周波数で指定されます。

  3. 包絡線信号を解析信号の絶対値の 2 倍として計算します。

関数は、包絡線信号の計算後、包絡線から DC バイアスを取り除き、FFT を使用して包絡線スペクトルを計算します。

参照

[1] Randall, Robert Bond. Vibration-Based Condition Monitoring. Chichester, UK: John Wiley & Sons, 2011.

拡張機能

C/C++ コード生成
MATLAB® Coder™ を使用して C および C++ コードを生成します。

バージョン履歴

R2017b で導入

参考

envelope | hilbert | orderspectrum

トピック

• 回転機の振動解析
• 転動体ベアリングの故障診断 (Predictive Maintenance Toolbox)