rpmordermap

次数分析のための次数-RPM マップ

構文

map = rpmordermap(x,fs,rpm)

map = rpmordermap(x,fs,rpm,res)

map = rpmordermap(___,Name,Value)

[map,order,rpm,time,res] = rpmordermap(___)

rpmordermap(___)

説明

map = rpmordermap(x,fs,rpm) は、入力ベクトル x に対する次数分析の結果、次数-RPM マップ行列 map を返します。x は、設定された rpm (1 分あたりの回転数で表される回転速度) において測定されます。fs はサンプル レート (Hz) です。map の各列は、rpm の各値に存在する次数の平方根平均二乗 (RMS) 振幅推定値を含みます。rpmordermap は、一定のサイクルあたりのサンプル数のレートで x をリサンプリングし、短時間フーリエ変換を使用してリサンプリングした信号のスペクトル成分を解析します。

map = rpmordermap(x,fs,rpm,res) は、マップの次数分解能を次数を単位として指定します。

map = rpmordermap(___,Name,Value) は、前の構文の入力引数に加えて Name,Value のペアを使用してオプションを指定します。

[map,order,rpm,time,res] = rpmordermap(___) は、次数、回転速度および次数マップが計算された時点をもつベクトルを返します。また、使用される次数分解能も返します。

出力引数なしの rpmordermap(___) は、次数マップを、回転速度および時間の関数として対話型 Figure 上にプロットします。

例

4 つの次数をもつチャープの次数-RPM マップ

600 Hz で 5 秒間サンプリングされたシミュレーション信号を作成します。テスト中のシステムは、観察期間中に回転速度を 1 秒あたり 10 回転から 40 回転に上昇させます。

タコメーターの読み取り値を生成します。

fs = 600;
t1 = 5;
t = 0:1/fs:t1;

f0 = 10;
f1 = 40;
rpm = 60*linspace(f0,f1,length(t));

信号は、1、0.5、4 および 6 の次数をもつ、調和的に関連した 4 つのチャープで構成されています。次数 4 のチャープの振幅は他のチャープの振幅の 2 倍です。チャープを生成するために、台形則を使用して位相を回転速度の積分で表します。

o1 = 1;
o2 = 0.5;
o3 = 4;
o4 = 6;

ph = 2*pi*cumtrapz(rpm/60)/fs;

x = [1 1 2 1]*cos([o1 o2 o3 o4]'*ph);

信号の次数-RPM マップを可視化します。

rpmordermap(x,fs,rpm)

ヘリコプターの振動データの次数-RPM マップ

ヘリコプターのコクピットに設置された加速度計から取得したシミュレーション データを解析します。

ヘリコプターのデータを読み込みます。振動測定値 vib は、500 Hz のレートで 10 秒間サンプリングされています。データの調査から、データには線形トレンドがあることがわかります。トレンドを除去して、トレンドによる次数推定の質の低下を防ぎます。

load('helidata.mat')

vib = detrend(vib);

非線形 RPM プロファイルをプロットします。ローターは、最大回転速度である 1 分あたり 27,600 回転に達するまで回転数を上昇させ、その後降下させます。

plot(t,rpm)
xlabel('Time (s)')
ylabel('RPM')

次数-RPM マップを計算します。0.015 の次数分解能を指定します。

[map,order,rpmOut,time] = rpmordermap(vib,fs,rpm,0.015);

マップを可視化します。

imagesc(time,order,map)
ax = gca;
ax.YDir = 'normal';
xlabel('Time (s)')
ylabel('Order')

より細かい次数分解能を使用して計算を繰り返します。rpmordermap の組み込み機能を使用してマップをプロットします。次数が低いほど、より明確に分解されます。

rpmordermap(vib,fs,rpm,0.005)

次数-RPM マップのウォーターフォール プロット

2 つの線形チャープと 1 つの 2 次チャープ (すべて 600 Hz で 5 秒間サンプリング) で構成される信号を生成します。テスト期間中、信号を発生させるシステムは回転速度を 1 秒あたり 10 回転から 40 回転に上昇させます。

タコメーターの読み取り値を生成します。

fs = 600;
t1 = 5;
t = 0:1/fs:t1;

f0 = 10;
f1 = 40;
rpm = 60*linspace(f0,f1,length(t));

線形チャープの次数は 1 および 2.5 です。次数 1 の成分の振幅は他の成分の振幅の 2 倍です。2 次チャープは次数 6 から始まり、測定の最後にこの次数に戻ります。振幅は 0.8 です。この情報を使用して信号を作成します。

o1 = 1;
o2 = 2.5;
o6 = 6;

x = 2*chirp(t,o1*f0,t1,o1*f1)+chirp(t,o2*f0,t1,o2*f1) + ...
    0.8*chirp(t,o6*f0,t1,o6*f1,'quadratic');

信号の次数-RPM マップを計算します。各測定セルのピーク振幅を使用します。0.25 次の分解能を指定します。サイドローブの減衰量が 80 dB のチェビシェフ ウィンドウをデータに適用します。

[map,or,rp] = rpmordermap(x,fs,rpm,0.25, ...
    'Amplitude','peak','Window',{'chebwin',80});

次数-RPM マップをウォーターフォール プロットとして描画します。

[OR,RP] = meshgrid(or,rp);
waterfall(OR,RP,map')

view(-15,45)
xlabel('Order')
ylabel('RPM')
zlabel('Amplitude')

対話型の次数-RPM マップ

出力引数なしで rpmordermap を呼び出して、対話型の次数-RPM マップをプロットします。

ファイル helidata.mat を読み込みます。このファイルには、ヘリコプターのコクピットに設置された加速度計から取得したシミュレーション振動データが含まれています。データは 500 Hz のレートで 10 秒間サンプリングされています。データの線形トレンドを除去します。rpmordermap を呼び出して、次数-RPM マップの対話型プロットを生成します。0.005 次の次数分解能を指定します。

load helidata.mat
rpmordermap(detrend(vib),fs,rpm,0.005)

Figure 下部の RPM 対時間のプロットに関するより詳細な説明は、アルゴリズムを参照してください。

Figure 内で十字カーソルを移動して、6 秒後の次数 0.053 の位置の RPM および RMS 振幅を特定します。

ツール バーの [X 軸ズーム] ボタン をクリックし、2 ～ 4 秒の時間領域を拡大します。RPM 対時間プロットのグレーの四角形部分は、関心領域を示します。この領域をスライドさせて時間全体にわたってパンできます。

[ウォーターフォール プロット] ボタン をクリックし、次数-RPM マップをウォーターフォール プロットとして表示します。可視性を向上させるため、[左に回転] ボタン を 3 回使用してプロットを時計回りに回転させます。5 ～ 7 秒の区間にパナーを移動します。

入力引数

x — 入力信号

ベクトル

入力信号。行ベクトルまたは列ベクトルとして指定します。

例: cos(pi/4*(0:159))+randn(1,160) は、ホワイト ガウス ノイズに含まれる正弦波を指定します。

fs — サンプル レート

正のスカラー

サンプル レート。正のスカラーとして Hz 単位で指定します。

rpm — 回転速度

正の値のベクトル

回転速度。1 分あたりの回転数で表される正の値のベクトルとして指定します。rpm は、x と同じ長さをもっていなければなりません。

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

• タコメーター パルス信号がない場合、rpmtrack を使用して振動信号から rpm を抽出します。

例: 100:10:3000 は、システムが、最初は 1 分あたり 100 回転で回転し、その後 10 回転ずつ回転数を上げて 1 分あたり 3000 回転になるように指定します。

res — 次数分解能

(15 × fs)/(16 × max(rpm)) (既定値) | 正のスカラー

次数-RPM マップの次数分解能。正のスカラーで指定します。res が指定されていない場合、rpmordermap により、サイクルあたりのサンプル数が一定の信号のサンプル レートを 256 で除算した値が設定されます。リサンプリングした入力信号長が十分ではない場合、関数はリサンプリングした信号長全体を使用して単一の次数推定を計算します。

実際の次数分解能は、指定された値とは若干異なる場合があります。詳細については、アルゴリズムを参照してください。

データ型: single | double

名前と値の引数

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

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

例: 'Scale','dB','Window','hann' は、次数マップの推定がデシベル単位でスケーリングされ、ハン ウィンドウを使用して決定されるように指定します。

Amplitude — 次数-RPM マップの振幅

'rms' (既定値) | 'peak' | 'power'

次数-RPM マップの振幅。'Amplitude' および 'rms'、'peak'、'power' のいずれかで構成されるコンマ区切りのペアとして指定します。

• 'rms' — 各推定次数の平方根平均二乗振幅を返します。

• 'peak' — 各推定次数のピーク振幅を返します。

• 'power' — 各推定次数のパワー レベルを返します。

OverlapPercent — 隣り合ったセグメント間のオーバーラップ率

50 (既定値) | 0 から 100 までのスカラー値

隣り合ったセグメント間のオーバーラップ率。'OverlapPercent' と、0 から 100 までのスカラー値で構成されるコンマ区切りのペアとして指定します。値 0 は隣り合ったセグメントがオーバーラップしないことを意味します。値 100 は隣り合ったセグメントが 1 サンプル分シフトされることを意味します。オーバーラップ率が大きくなると、滑らかなマップが生成されますが、計算時間が長くなります。詳細については、アルゴリズムを参照してください。

データ型: double | single

Scale — 次数-RPM マップのスケーリング

'linear' (既定値) | 'dB'

次数-RPM マップのスケーリング。'Scale' と 'linear' または 'dB' のいずれかで構成されるコンマ区切りのペアとして指定します。

• 'linear' — 線形にスケーリングされたマップを返します。

• 'dB' — 値をデシベル単位で表した対数マップを返します。

Window — 解析ウィンドウ

'flattopwin' (既定値) | 'chebwin' | 'hamming' | 'hann' | 'kaiser' | 'rectwin'

解析ウィンドウ。'Window' と次の値のいずれかで構成されるコンマ区切りペアとして指定します。

• 'flattopwin' はフラット トップ ウィンドウを指定します。詳細については、flattopwin を参照してください。

• 'chebwin' はチェビシェフ ウィンドウを指定します。cell 配列を使用して、サイドローブの減衰をデシベル単位で指定します。サイドローブの減衰量は 45 dB より大きくなければなりません。指定されていない場合、既定値の 100 dB が使用されます。詳細については、chebwin を参照してください。

• 'hamming' はハミング ウィンドウを指定します。詳細については、hamming を参照してください。

• 'hann' はハン ウィンドウを指定します。詳細については、hann を参照してください。

• 'kaiser' はカイザー ウィンドウを指定します。cell 配列を使用して、形状パラメーター β を指定します。形状パラメーターは正のスカラーでなければなりません。指定されていない場合、既定値の 0.5 が使用されます。詳細については、kaiser を参照してください。

• 'rectwin' は箱型ウィンドウを指定します。詳細については、rectwin を参照してください。

例: 'Window','chebwin' は、サイドローブの減衰量が 100 dB であるチェビシェフ ウィンドウを指定します。

例: 'Window',{'chebwin',60} は、サイドローブの減衰量が 60 dB であるチェビシェフ ウィンドウを指定します。

例: 'Window','kaiser' は、0.5 の形状パラメーターをもつカイザー ウィンドウを指定します。

例: 'Window',{'kaiser',1} は、1 の形状パラメーターをもつカイザー ウィンドウを指定します。

データ型: char | string | cell

出力引数

map — 次数-RPM マップ

行列

次数-RPM マップ。行列として返されます。

order — 次数

ベクトル

次数。ベクトルとして返されます。

rpm — 回転速度

ベクトル

回転速度。ベクトルとして返されます。

time — 時点

ベクトル

時点。ベクトルとして返されます。

res — 次数分解能

スカラー

次数分解能。スカラー値として返されます。

アルゴリズム

次数分析は、回転自体から生じる回転系における振動の調査です。多くの場合、これらの振動の周波数は回転速度に比例します。比例定数は "次数" です。

回転速度は通常、ほとんどの実験的状況下で個別に測定され、時間と共に変化します。回転から生じる振動を正確に解析するには、測定した信号をリサンプリングして内挿を行い、サイクルあたりのサンプル数を一定にする必要があります。この過程によって、回転速度の定数倍である周波数をもつ信号成分が一定のトーンに変換されます。変換により、周波数が時間によって急激に変化する場合に生じるスペクトル成分の不鮮明さが低減します。

関数 rpmordermap は以下のステップを実行します。

1. cumtrapz を使用し、位相角を回転速度の時間積分として推定します。

   $$\varphi (t)={\displaystyle {\int }_0^t\frac{\mathrm{RPM}(\tau )}{60}}d\tau .$$

2. resample を使用し、信号をアップサンプリングしてローパス フィルターで処理します。このステップにより、高周波数成分のエイリアシングを行わずに、サンプリングされていない時点の信号を関数で内挿できるようになります。rpmordermap は、信号を 15 の係数でアップサンプリングします。

3. interp1 を使用し、アップサンプリングされた信号を位相領域の等間隔グリッド上に線形に内挿します。測定で使用可能な最高次数は、サンプル レートと、系が到達する最高回転速度によって決定されます。

   $$O_{\max }=\frac{f_s/2}{\max \left(\frac{\mathrm{RPM}}{60}\right)}.$$

   この最高次数を正確に取得するには、少なくとも O_max の 2 倍で信号をサンプリングする必要があります。より正確な結果を得るため、rpmordermap は追加の係数 4 でオーバーサンプリングします。結果の位相-領域サンプル レート f_p は次になります。

   $$f_p=4\times 2O_{\max }=4\times 2\frac{f_s/2}{\max \left(\frac{\mathrm{RPM}}{60}\right)}.$$

   既定では、rpmordermap は次数-RPM 行列をターゲット次数分解能で計算するように構成されています。

   $$r=\frac{f_p}{256}=\frac{4\times 60}{256}\frac{2\times f_s/2}{\max (\text{RPM})}=\frac{15}{16}\frac{f_s}{\max (\text{RPM})},$$

   ただし、入力引数 res を使用して別の値を指定することができます。

4. spectrogram を使用し、内挿された信号の短時間フーリエ変換 (STFT) を計算します。既定では、関数は、信号を L サンプルのセグメントに分割し、それぞれにフラット トップ ウィンドウを適用します。隣り合ったセグメント間でオーバーラップするサンプルの数は次になります。

   $$N_{\text{overlap}}=\min \left(\lceil \frac{p_{\text{overlap}}}{100}\times L\rceil ,L-1\right)$$

   ここで、p_overlap は 'OverlapPercent' の名前と値のペアを使用して指定され、既定値は 50% です。DFT 長は L に設定されます。分解能は、以下によってサンプル レートとセグメント長に関連付けられます。

   $$r=\frac{f_p}{L}\times \text{ENBW},$$

   ここで、ウィンドウの等価ノイズ帯域幅 ENBW は enbw を使用して計算されます。分解能を調節して、近接した次数を識別します。r の値を小さくすると、より大きなセグメント長が必要になります。ENBW 自体は L によって異なるため、L は指定された r と f_p に対して再帰的に計算されなければなりません。通常、結果の L の値は整数ではないため、rpmordermap は ceil を使用して値を丸めます。したがって実際の次数分解能は、指定されたターゲット値とは若干異なる場合があります。特定の分解能を得る必要がある場合、信号がそのために十分なサンプルをもっていることを確認してください。

対話型 rpmordermap ウィンドウの下部にある RPM 対時間プロットの赤い点は、ウィンドウが適用された各セグメントの右端に対応します。プロットの青ラインは、時間の関数として描画された RPM 信号です。

参照

[1] Brandt, Anders. Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Chichester, UK: John Wiley & Sons, 2011.

拡張機能

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

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2015b で導入

R2023b: gpuArray オブジェクトの使用

関数 rpmordermap は gpuArray オブジェクトをサポートしています。この機能を使用するには、Parallel Computing Toolbox™ が必要です。

参考

orderspectrum | ordertrack | orderwaveform | rpmfreqmap | rpmtrack | spectrogram | tachorpm