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

rpmtrack

振動信号の RPM プロファイルの追跡および抽出

構文

rpm = rpmtrack(x,fs,order,p)

rpm = rpmtrack(xt,order,p)

rpm = rpmtrack(___,Name=Value)

[rpm,tout] = rpmtrack(___)

rpmtrack(___)

説明

rpm = rpmtrack(x,fs,order,p) は、fs のレートでサンプリングした振動信号 x からの回転速度 rpm の時間依存推定を返します。

2 列行列 p には、与えられた order に対応する時間-周波数リッジ上にある一連の点が含まれます。p の各行は、1 対の座標を指定します。order および p の両方を指定せずに rpmtrack を呼び出すと、関数によって時間-周波数マップを表示する対話型プロットが開かれ、点を選択できるようになります。

タコメーター パルス信号がある場合、tachorpm を使用して rpm を直接抽出します。

rpm = rpmtrack(xt,order,p) は、MATLAB^® timetable xt に保存されている信号から回転速度の時間依存推定を返します。

rpm = rpmtrack(___,Name=Value) は、名前と値の引数を使用して、上記の任意の構文に追加オプションを指定します。オプションには、時間-周波数マップや RPM プロファイルの開始時間の推定に使用するメソッドがあります。

[rpm,tout] = rpmtrack(___) は、RPM プロファイルが計算される時間ベクトルも返します。

出力引数なしの rpmtrack(___) は、パワー時間-周波数マップおよび推定 RPM プロファイルを対話型 Figure 上にプロットします。

例

振動信号の RPM プロファイル

3 つの高調波成分を含む振動信号を生成します。信号は 1 kHz で 16 秒間サンプリングされます。この信号の瞬時周波数は、エンジンの試運転と惰行試験に似ています。台形則を使用して周波数を積分することにより、瞬間位相を計算します。

fs = 1000;
t = 0:1/fs:16;

ifq = 20 + t.^6.*exp(-t);
phi = 2*pi*cumtrapz(t,ifq);

信号の高調波成分は、1、2 および 3 の次数に対応します。次数 2 の正弦波の振幅は他の 2 倍です。

ol = [1 2 3];
amp = [5 10 5];

vib = amp*cos(ol'.*phi);

次数 2 のリッジ上の点を使用して信号の RPM プロファイルを抽出し、可視化します。

time = 3;
order = 2;
p = [time order*ifq(t==time)];

rpmtrack(vib,fs,order,p)

空ぶかしエンジンの RPM プロファイル

自動車エンジンの空ぶかしによって発生する振動に似た信号を生成します。この信号は 1 kHz で 30 秒間サンプリングされ、それぞれ振幅 5、4、および 0.5 で次数が 1、2.4、および 3 である 3 つの高調波成分が含まれています。信号を単位分散のホワイト ガウス ノイズに組み込み、MATLAB® timetable に保存します。瞬時周波数を 60 倍にして RPM プロファイルを取得します。RPM プロファイルをプロットします。

fs = 1000;
t = (0:1/fs:30)';

fit = @(a,x) (t-x).^6.*exp(-(t-x)).*((t-x)>=0)*a';

fis = fit([0.4 1 0.6 1],[0 6 13 17]);
phi = 2*pi*cumtrapz(t,fis);

ol = [1 2.4 3];
amp = [5 4 0.5]';
vib = cos(phi.*ol)*amp + randn(size(t));

xt = timetable(seconds(t),vib);

plot(t,fis*60)

振動信号から RPM プロファイルを導きます。5 秒間隔で 4 ポイントを使用して 2.4 の次数に対応するリッジを指定します。出力 timetable の概要を表示します。

ndx = (5:5:20)*fs;
order = ol(2);

p = [t(ndx) order*fis(ndx)];

rpmest = rpmtrack(xt,order,p);

summary(rpmest)

RowTimes:

    tout: 30001x1 duration
        Values:
            Min           0 sec     
            Median        15 sec    
            Max           30 sec    
            TimeStep      0.001 sec 

Variables:

    rpm: 30001x1 double

        Values:

            Min       2.2204e-16
            Median        4327.2
            Max           8879.8

再構成した RPM プロファイルと再構成に使用したポイントをプロットします。

hold on
plot(seconds(rpmest.tout),rpmest.rpm,'.-')
plot(t(ndx),fis(ndx)*60,'ok')
hold off
legend('Original','Reconstructed','Ridge points','Location','northwest')

抽出した RPM プロファイルを使用して信号の次数-RPM マップを生成します。

rpmordermap(vib,fs,rpmest.rpm)

信号を構成する時間領域の波形を再構成してプロットします。過渡状態が減衰した後に発生する時間間隔を拡大します。

xrc = orderwaveform(vib,fs,rpmest.rpm,ol);

figure
plot(t,xrc)
legend([repmat('Order = ',[3 1]) num2str(ol')])
xlim([5 20])

ファン スイッチオフ RPM プロファイル

スイッチオフ後ファン ブレードが遅くなっていくときの RPM プロファイルを推定します。

20,000 rpm で回転する工業用天井型扇風機をオフにします。空気抵抗 (ベアリングの摩擦からの寄与は無視できるとします) により約 6 秒でファンのローターが停止します。高速カメラでファン ブレードの 1 つの x 座標を 1 kHz のレートで測定します。

fs = 1000;
t = 0:1/fs:6-1/fs;

rpm0 = 20000;

このファン ブレードを、ローターを中心として半径 50 cm を回転する質点として理想化します。ブレードには速度に比例する抵抗力がかかり、その結果、位相角は次の式で求められます。

$$\varphi =2\pi f_{0}T(1-e^{-t/T}),$$

ここで、${\mathit{f}}_0$ は初期周波数、$\mathit{T}=0.75$ 秒は減衰時間です。

a = 0.5;
f0 = rpm0/60;
T = 0.75;

phi = 2*pi*f0*T*(1-exp(-t/T));

ブレードの x 座標と y 座標を計算し、プロットします。分散 $0.1^2$ のホワイト ガウス ノイズを追加します。

x = a*cos(phi) + randn(size(phi))/10;
y = a*sin(phi) + randn(size(phi))/10;

plot(t,x,t,y)

関数 rpmtrack を使用して、RPM プロファイルを決定します。関数

rpmtrack(x,fs)

をコマンド ラインで入力し、対話型 Figure を開きます。

スライダーを使用して時間-周波数マップの周波数分解能を約 11 Hz に調節します。信号成分が次数 1 に相当すると仮定し、リッジ抽出の終了時間を 3.0 秒に設定します。時間-周波数マップで十字カーソルを使用し、[追加] ボタンを使用してリッジ上に 3 つのポイントを追加します または、選択した場所でカーソルをダブルクリックしてポイントを追加します。[推定] をクリックして RPM プロファイルを追跡および抽出します。

RPM プロファイルが指数関数的に減衰することを確認します。[エクスポート] タブで、[エクスポート] をクリックして [Generate MATLAB Script] を選択します。スクリプトがエディターで開きます。

% MATLAB Code from rpmtrack GUI

% Generated by MATLAB 9.12 and Signal Processing Toolbox 8.7

% Generated on 12-Oct-2021 09:36:49

% Set sample rate
fs = 1000.0000;
% Set order of ridge of interest
order = 1.0000;
% Set ridge points on ridge of interest
ridgePoints = [...
    0.4077 194.6023;...
    0.9781 89.4886;...
    2.3678 15.6250];
% Estimate RPM
[rpmOut,tOut] = rpmtrack(x,fs,order,ridgePoints,...
    'Method','stft',...
    'FrequencyResolution',11.1612,...
    'PowerPenalty',Inf,...
    'FrequencyPenalty',0.0000,...
    'StartTime',0.0000,...
    'EndTime',3.0000);

スクリプトを実行します。RPM プロファイルを片対数プロットに表示します。

semilogy(tOut,rpmOut)
ylim([500 20000])

入力引数

x — 入力信号

ベクトル

入力信号。ベクトルで指定します。

例: cos(pi/4*(0:159))+randn(1,160) では、2π Hz でサンプリングしたノイズを含んだ正弦波を指定します。

データ型: single | double

fs — サンプル レート

正の実数のスカラー

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

データ型: single | double

order — リッジ次数

正の実数のスカラー

リッジ次数。正の実数スカラーとして指定します。

データ型: single | double

p — リッジ ポイント

2 列の行列

リッジ ポイント。各行上に 1 つの時間-周波数座標を含む 2 列の行列として指定します。座標は、目的の次数リッジに属する時間-周波数マップ上のポイントを表現します。

データ型: single | double

xt — 入力 timetable

timetable

入力 timetable。xt には、duration 型の、増加する有限の、等間隔に配置された行時間が含まれなければなりません。timetable には、信号値をもつ数値データ ベクトルを 1 つのみ含めなければなりません。

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

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

データ型: single | double

名前と値の引数

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

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

例: "Method","fsst","PowerPenalty",10 は、リッジ上の隣接するポイント間で最大 10 デシベルのパワー差を許容して、フーリエ シンクロスクイーズド変換を使用して時間-周波数マップを推定するように指定します。

Method — 時間-周波数マップのタイプ

"stft" (既定値) | "fsst"

推定プロセスで使用する時間-周波数マップのタイプ。"stft" または "fsst" のいずれかで指定します。

• "stft" — 短時間フーリエ変換を使用してパワー スペクトログラム時間-周波数マップを計算します。短時間フーリエ変換の詳細については、pspectrum を参照してください。

• "fsst" — フーリエ シンクロスクイーズド変換を使用して時間-周波数マップを計算します。フーリエ シンクロスクイーズド変換の詳細については、fsst を参照してください。

FrequencyResolution — 周波数分解能帯域幅

数値スカラー

時間-周波数マップの計算に使用する周波数分解能帯域幅。Hz 単位で表される数値スカラーとして指定します。

データ型: single | double

PowerPenalty — 隣接するリッジ ポイント間のパワーの最大差

Inf (既定値) | dB 単位の数値スカラー

隣接するリッジ ポイント間のパワーの最大差。dB 単位で表される数値スカラーとして指定します。

このパラメーターを使用して、rpmtrack のリッジ抽出アルゴリズムが対応する次数の正しいリッジを検出することを確認します。PowerPenalty は、目的の次数のリッジが他のリッジと交差している、または周波数的に他のリッジに非常に近いが、パワー レベルが異なる場合に役立ちます。

データ型: single | double

FrequencyPenalty — おおまかなリッジ抽出におけるペナルティ

0 (既定値) | 非負のスカラー

おおまかなリッジ抽出におけるペナルティ。非負のスカラーとして指定します。

このパラメーターを使用して、rpmtrack のリッジ抽出アルゴリズムが、リッジ推定を不適切な時間-周波数位置に移動させる可能性がある大きいジャンプを回避することを確認します。FrequencyPenalty は、交差している次数リッジ、または周波数的に間隔が近い次数リッジを区別する場合に役立ちます。

データ型: single | double

StartTime — RPM プロファイル推定の開始時間

入力信号開始時間 (既定値) | 秒単位のスカラー | duration スカラー

RPM プロファイル推定の開始時間。数値または duration スカラーとして指定します。

データ型: single | double | duration

EndTime — RPM プロファイル推定の終了時間

入力信号終了時間 (既定値) | 秒単位のスカラー | duration スカラー

RPM プロファイル推定の終了時間。数値または duration スカラーとして指定します。

データ型: single | double | duration

出力引数

rpm — 回転速度推定

ベクトル | timetable

回転速度推定。1 分あたりの回転数で表されるベクトルとして返されます。

rpmtrack への入力が timetable である場合、rpm も単一変数ラベル付き rpm を含む timetable になります。timetable の行時間は、ラベル付き tout で duration 型です。

tout — 時間値

ベクトル

RPM プロファイルが推定される時間値。ベクトルとして返されます。

アルゴリズム

rpmtrack は、2 ステップ (おおまか-細かい) の推定手法を使用します。

1. x の時間-周波数マップを計算して、リッジ上の指定された一連のポイント p、そのリッジに対応する order、およびオプションのペナルティ パラメーター PowerPenalty と FrequencyPenalty に基づいて時間-周波数リッジを抽出します。抽出されたリッジは、RPM プロファイルのおおまかな推定を示します。

2. ボルド・カルマン フィルターを使用して抽出されたリッジに対応する次数波形を計算し、この波形から新しい時間-周波数マップを計算します。新しい時間-周波数マップから分離した次数リッジは、RPM プロファイルの細かな推定を示します。

参照

[1] Urbanek, Jacek, Tomasz Barszcz, and Jerome Antoni. "A Two-Step Procedure for Estimation of Instantaneous Rotational Speed with Large Fluctuations." Mechanical Systems and Signal Processing. Vol. 38, 2013, pp. 96–102.

拡張機能

C/C++ コード生成

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

使用上の注意および制限:

• 名前と値の引数の文字または文字列の入力は、コンパイル時に定数でなければなりません。

バージョン履歴

R2018a で導入

参考

関数

• ordertrack | rpmfreqmap | rpmordermap | tachorpm