シミュレーション データを使用した複数クラス故障検出
この例では、Simulink モデルを使用して健全な状態のデータと故障状態のデータを生成する方法を説明します。このデータを使用して、さまざまな故障の組み合わせを検出する複数クラス分類器を開発します。この例は 3 重往復ポンプ モデルを使用し、漏れ、閉塞、およびベアリングの故障を含みます。
モデルの設定
この例では、zip ファイルに保存されている多くのサポート ファイルを使用します。ファイルを解凍してサポート ファイルにアクセスし、モデル パラメーターを読み込み、往復ポンプ ライブラリを作成します。
if ~exist('+mech_hydro_forcesPS','dir') unzip('pdmRecipPump_supportingfiles.zip') end % Load Parameters pdmRecipPump_Parameters %Pump CAT_Pump_1051_DataFile_imported %CAD % Create Simscape library if needed if exist('mech_hydro_forcesPS_Lib','file')~=4 ssc_build mech_hydro_forcesPS end
往復ポンプ モデル
往復ポンプは、電気モーター、ポンプ ハウジング、ポンプ クランク、およびポンプ プランジャーで構成されています。
mdl = 'pdmRecipPump';
open_system(mdl)
open_system([mdl,'/Pump'])
ポンプ モデルは、シリンダーの漏れ、吸込口の閉塞、およびベアリング摩擦の増加という 3 タイプの故障をモデル化するように設定されます。これらの故障はワークスペース変数としてパラメーター化され、ポンプ ブロック ダイアログによって設定されます。
故障データと健全データのシミュレーション
3 つの故障タイプそれぞれについて、故障なしから重大な故障までの故障重大度を表す値の配列を作成します。
% Define fault parameter variations numParValues = 10; leak_area_set_factor = linspace(0.00,0.036,numParValues); leak_area_set = leak_area_set_factor*TRP_Par.Check_Valve.In.Max_Area; leak_area_set = max(leak_area_set,1e-9); % Leakage area cannot be 0 blockinfactor_set = linspace(0.8,0.53,numParValues); bearingfactor_set = linspace(0,6e-4,numParValues);
ポンプ モデルはノイズを含めるように設定されているため、同じ故障パラメーター値でモデルを実行すると、結果としてさまざまなシミュレーション出力が得られます。これは、同じ故障状態と重大度に対して複数のシミュレーション結果が可能であることを意味するので、分類器の開発に役立ちます。そのような結果が得られるシミュレーションを設定するには、故障パラメーター値のベクトルを作成します。ここで値は、故障なし、1 つの故障、2 つの故障の組み合わせ、および 3 つの故障の組み合わせを表します。故障なし、1 つの故障、などの各グループについて、上記で定義されている故障パラメーター値から、故障値の 125 の組み合わせを作成します。これにより故障パラメーター値の合計 1000 の組み合わせが得られます。
nPerGroup = 125; % Number of elements in each fault group % No fault simulations leakArea = repmat(leak_area_set(1),nPerGroup,1); blockingFactor = repmat(blockinfactor_set(1),nPerGroup,1); bearingFactor = repmat(bearingfactor_set(1),nPerGroup,1); % Single fault simulations idx = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idx)']; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)]; idx = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;blockinfactor_set(idx)']; bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)]; idx = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;bearingfactor_set(idx)']; % Double fault simulations idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idxA)']; blockingFactor = [blockingFactor;blockinfactor_set(idxB)']; bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)]; idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idxA)']; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;bearingfactor_set(idxB)']; idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;blockinfactor_set(idxA)']; bearingFactor = [bearingFactor;bearingfactor_set(idxB)']; % Triple fault simulations idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); idxC = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idxA)']; blockingFactor = [blockingFactor;blockinfactor_set(idxB)']; bearingFactor = [bearingFactor;bearingfactor_set(idxC)'];
故障パラメーターの組み合わせを使用して Simulink.SimulationInput オブジェクトを作成します。各シミュレーション入力につき、異なる結果が生成されるようにランダム シードを必ず異なる値に設定します。
for ct = numel(leakArea):-1:1 simInput(ct) = Simulink.SimulationInput(mdl); simInput(ct) = setVariable(simInput(ct),'leak_cyl_area_WKSP',leakArea(ct)); simInput(ct) = setVariable(simInput(ct),'block_in_factor_WKSP',blockingFactor(ct)); simInput(ct) = setVariable(simInput(ct),'bearing_fault_frict_WKSP',bearingFactor(ct)); simInput(ct) = setVariable(simInput(ct),'noise_seed_offset_WKSP',ct-1); end
関数 generateSimulationEnsemble を使用して、上で定義された Simulink.SimulationInput オブジェクトによって定義されるシミュレーションを実行し、ローカル サブフォルダーに結果を格納します。次に、格納された結果から simulationEnsembleDatastore を作成します。
これら 1000 のシミュレーションを並列で実行する場合、標準のデスクトップでは約 1 時間かかり、約 620 MB のデータが生成されることに注意してください。便宜上、最初の 10 個のシミュレーションのみを実行するオプションが用意されています。
% Run the simulation and create an ensemble to manage the simulation results runAll = true; if runAll [ok,e] = generateSimulationEnsemble(simInput,fullfile('.','Data'),'UseParallel',true); else [ok,e] = generateSimulationEnsemble(simInput(1:10),fullfile('.','Data')); %#ok<UNRCH> end
[09-Apr-2018 09:01:38] Checking for availability of parallel pool... [09-Apr-2018 09:01:38] Loading Simulink on parallel workers... Analyzing and transferring files to the workers ...done. [09-Apr-2018 09:01:38] Configuring simulation cache folder on parallel workers... [09-Apr-2018 09:01:38] Running SetupFcn on parallel workers... [09-Apr-2018 09:01:39] Loading model on parallel workers... [09-Apr-2018 09:01:39] Transferring base workspace variables used in the model to parallel workers... [09-Apr-2018 09:01:41] Running simulations... [09-Apr-2018 09:02:28] Completed 1 of 1000 simulation runs [09-Apr-2018 09:02:33] Completed 2 of 1000 simulation runs [09-Apr-2018 09:02:37] Completed 3 of 1000 simulation runs [09-Apr-2018 09:02:41] Completed 4 of 1000 simulation runs [09-Apr-2018 09:02:46] Completed 5 of 1000 simulation runs [09-Apr-2018 09:02:49] Completed 6 of 1000 simulation runs [09-Apr-2018 09:02:54] Completed 7 of 1000 simulation runs [09-Apr-2018 09:02:58] Completed 8 of 1000 simulation runs [09-Apr-2018 09:03:01] Completed 9 of 1000 simulation runs...
ens = simulationEnsembleDatastore(fullfile('.','Data'));
シミュレーション結果からの特徴の処理と抽出
このモデルはポンプの吐出圧力、出流量、モーター速度、およびモーター電流をログに記録するよう設定されます。
ens.DataVariables
ans = 8×1 string array
"SimulationInput"
"SimulationMetadata"
"iMotor_meas"
"pIn_meas"
"pOut_meas"
"qIn_meas"
"qOut_meas"
"wMotor_meas"
アンサンブルのメンバーごとに、ポンプ出流量を前処理し、ポンプ出流量に基づいて状態インジケーターを計算します。状態インジケーターは後で故障分類のために使用されます。前処理のため、出流量の最初の 0.8 秒を削除します。これにはシミュレーションおよびポンプ始動からの過渡状態が含まれるためです。前処理の一部として、出流量のパワー スペクトルを計算し、SimulationInput を使って故障変数の値を返します。
対象とする変数だけが読み取りで返されるようにアンサンブルを設定し、この例の最後で定義されている関数 preprocess を呼び出します。
ens.SelectedVariables = ["qOut_meas", "SimulationInput"]; reset(ens) data = read(ens)
data=1×2 table
qOut_meas SimulationInput
__________________ ______________________________
[2001×1 timetable] [1×1 Simulink.SimulationInput]
[flow,flowP,flowF,faultValues] = preprocess(data);
流量と流量スペクトルをプロットします。プロットされるデータは故障なしの状態のものです。
% Figure with nominal subplot(211); plot(flow.Time,flow.Data); subplot(212); semilogx(flowF,pow2db(flowP)); xlabel('Hz')
流量スペクトルは、期待された周波数での共振ピークを示しています。具体的には、ポンプのモーター速度は 950 rpm、つまり 15.833 Hz であり、ポンプはシリンダーを 3 つもつため、流量は 3*15.833 Hz つまり 47.5 Hz で基本波であり、47.5 Hz の倍数において高調波であると予想されます。流量スペクトルは、期待される共振ピークを明確に示しています。ポンプの 1 つのシリンダーに故障が生じると、ポンプのモーター速度 15.833 Hz およびその高調波で共振が発生します。
流量スペクトルおよび低速信号は、可能な状態インジケーターをある程度示唆しています。具体的には、平均や分散などの一般的な信号統計量、およびスペクトル特性です。ピーク振幅の周波数、15.833 Hz 周辺のエネルギー、47.5 Hz 周辺のエネルギー、および 100 Hz を超えるエネルギーなどの、期待される高調波に関連するスペクトル状態インジケーターが計算されます。スペクトル尖度ピークの周波数も計算されます。
状態インジケーターのデータ変数および故障変数値の状態変数を使用して、アンサンブルを設定します。その後、関数 extractCI を呼び出して特徴を計算し、writeToLastMemberRead コマンドを使用して特徴および故障変数の値をアンサンブルに追加します。関数 extractCI はこの例の最後で定義されています。
ens.DataVariables = [ens.DataVariables; ... "fPeak"; "pLow"; "pMid"; "pHigh"; "pKurtosis"; ... "qMean"; "qVar"; "qSkewness"; "qKurtosis"; ... "qPeak2Peak"; "qCrest"; "qRMS"; "qMAD"; "qCSRange"]; ens.ConditionVariables = ["LeakFault","BlockingFault","BearingFault"]; feat = extractCI(flow,flowP,flowF); dataToWrite = [faultValues, feat]; writeToLastMemberRead(ens,dataToWrite{:})
上記のコードは、シミュレーション アンサンブルの最初のメンバーの状態インジケーターを前処理して計算します。アンサンブルの hasdata コマンドを使用して、アンサンブルのすべてのメンバーについてこれを繰り返します。さまざまな故障状態下におけるシミュレーション結果を把握するには、アンサンブルを 100 要素ごとにプロットします。
%Figure with nominal and faults figure, subplot(211); lFlow = plot(flow.Time,flow.Data,'LineWidth',2); subplot(212); lFlowP = semilogx(flowF,pow2db(flowP),'LineWidth',2); xlabel('Hz') ct = 1; lColors = get(lFlow.Parent,'ColorOrder'); lIdx = 2; % Loop over all members in the ensemble, preprocess % and compute the features for each member while hasdata(ens) % Read member data data = read(ens); % Preprocess and extract features from the member data [flow,flowP,flowF,faultValues] = preprocess(data); feat = extractCI(flow,flowP,flowF); % Add the extracted feature values to the member data dataToWrite = [faultValues, feat]; writeToLastMemberRead(ens,dataToWrite{:}) % Plot member signal and spectrum for every 100th member if mod(ct,100) == 0 line('Parent',lFlow.Parent,'XData',flow.Time,'YData',flow.Data, ... 'Color', lColors(lIdx,:)); line('Parent',lFlowP.Parent,'XData',flowF,'YData',pow2db(flowP), ... 'Color', lColors(lIdx,:)); if lIdx == size(lColors,1) lIdx = 1; else lIdx = lIdx+1; end end ct = ct + 1; end
さまざまな故障状態および重大度の下において、期待される周波数でスペクトルに高調波が含まれることに注意してください。
ポンプの故障の検出と分類
前の節では、シミュレーション アンサンブルの全メンバーの流量信号から状態インジケーターを前処理して計算しました。これはさまざまな故障の組み合わせと重大度によるシミュレーション結果に対応します。状態インジケーターを使用して、ポンプの流量信号からポンプの故障を検出し分類することができます。
状態インジケーターを読み取るようにシミュレーション アンサンブルを設定し、tall コマンドと gather コマンドを使用してすべての状態インジケーターと故障変数値をメモリに読み込むことができます｡
% Get data to design a classifier. reset(ens) ens.SelectedVariables = [... "fPeak","pLow","pMid","pHigh","pKurtosis",... "qMean","qVar","qSkewness","qKurtosis",... "qPeak2Peak","qCrest","qRMS","qMAD","qCSRange",... "LeakFault","BlockingFault","BearingFault"]; idxLastFeature = 14; % Load the condition indicator data into memory data = gather(tall(ens));
Starting parallel pool (parpool) using the 'local' profile ...
Preserving jobs with IDs: 1 2 because they contain crash dump files.
You can use 'delete(myCluster.Jobs)' to remove all jobs created with profile local. To create 'myCluster' use 'myCluster = parcluster('local')'.
connected to 6 workers.
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 45 sec
Evaluation completed in 45 sec
data(1:10,:)
ans=10×17 table
fPeak pLow pMid pHigh pKurtosis qMean qVar qSkewness qKurtosis qPeak2Peak qCrest qRMS qMAD qCSRange LeakFault BlockingFault BearingFault
______ _______ ______ ______ _________ ______ ______ _________ _________ __________ ______ ______ ______ ________ _________ _____________ ____________
43.909 0.86472 117.63 18.874 276.49 35.572 7.5242 -0.72832 2.7738 13.835 1.1494 35.677 2.2326 42690 1e-09 0.8 0
43.909 0.44477 125.92 18.899 12.417 35.576 7.869 -0.7094 2.6338 13.335 1.1449 35.686 2.3204 42697 1e-09 0.8 0
43.909 1.1782 137.99 17.526 11.589 35.573 7.4367 -0.72208 2.7136 12.641 1.1395 35.678 2.2407 42695 1e-09 0.8 0
44.151 156.74 173.84 21.073 199.5 33.768 12.466 -0.30256 2.4782 17.446 1.2138 33.952 2.8582 40518 2.4e-06 0.8 0
43.848 0.71756 110.92 18.579 197.02 35.563 7.5781 -0.72377 2.793 14.14 1.1504 35.669 2.2671 42682 1e-09 0.8 0
43.909 0.43673 119.56 20.003 11.589 35.57 7.5028 -0.74797 2.7913 13.833 1.1551 35.676 2.2442 42689 1e-09 0.8 0
43.788 0.31617 135.3 19.724 476.82 35.568 7.4406 -0.70964 2.6884 14.685 1.1473 35.673 2.2392 42687 1e-09 0.8 0
43.848 0.72747 121.63 19.733 11.589 35.523 7.791 -0.72736 2.7864 14.043 1.1469 35.633 2.2722 42633 1e-09 0.8 0
43.848 0.62777 128.85 19.244 11.589 35.541 7.5698 -0.6953 2.6942 13.451 1.1415 35.647 2.2603 42654 1e-09 0.8 0
43.848 0.4631 134.83 18.918 12.417 35.561 7.8607 -0.68417 2.6664 13.885 1.1504 35.671 2.3078 42681 1e-09 0.8 0
アンサンブル メンバーごとの故障変数値 (データ テーブルの 1 行) を故障フラグに変換し、これらの故障フラグを、各メンバーのさまざまな故障ステータスを捕捉する 1 つのフラグに組み合わせることが可能です。
% Convert the fault variable values to flags
data.LeakFlag = data.LeakFault > 1e-6;
data.BlockingFlag = data.BlockingFault < 0.8;
data.BearingFlag = data.BearingFault > 0;
data.CombinedFlag = data.LeakFlag+2*data.BlockingFlag+4*data.BearingFlag;状態インジケーターを入力として受け取り、組み合わせた故障フラグを返す分類器を作成します。2 次多項式カーネルを使用するサポート ベクター マシンを学習させます。cvpartition コマンドを使用して、アンサンブルのメンバーを学習用のセットと検証用のセットに分割します。
rng('default') % for reproducibility predictors = data(:,1:idxLastFeature); response = data.CombinedFlag; cvp = cvpartition(size(predictors,1),'KFold',5); % Create and train the classifier template = templateSVM(... 'KernelFunction', 'polynomial', ... 'PolynomialOrder', 2, ... 'KernelScale', 'auto', ... 'BoxConstraint', 1, ... 'Standardize', true); combinedClassifier = fitcecoc(... predictors(cvp.training(1),:), ... response(cvp.training(1),:), ... 'Learners', template, ... 'Coding', 'onevsone', ... 'ClassNames', [0; 1; 2; 3; 4; 5; 6; 7]);
検証データを使用して学習後の分類器の性能をチェックし、結果を混同プロットにプロットします。
% Check performance by computing and plotting the confusion matrix actualValue = response(cvp.test(1),:); predictedValue = predict(combinedClassifier, predictors(cvp.test(1),:)); confdata = confusionmat(actualValue,predictedValue); figure, labels = {'None', 'Leak','Blocking', 'Leak & Blocking', 'Bearing', ... 'Bearing & Leak', 'Bearing & Blocking', 'All'}; h = heatmap(confdata, ... 'YLabel', 'Actual leak fault', ... 'YDisplayLabels', labels, ... 'XLabel', 'Predicted fault', ... 'XDisplayLabels', labels, ... 'ColorbarVisible','off');
混同プロットは、故障の各組み合わせにつき、故障の組み合わせが正しく予測された回数 (プロットの対角エントリ) と、故障の組み合わせが誤って予測された回数 (非対角エントリ) を示しています。
混同プロットから、分類器は一部の故障状態を正しく分類しなかった (非対角項) ことがわかります。しかし、故障なしの状態は正しく予測されています。故障があっても故障なしの状態が予測された箇所 (1 列目) がいくつかありますが、それ以外には、厳密には正しい故障状態ではないにしても、故障が予測されています。全体的な検証の正確性は 84% で、故障があることを予測する正確性は 98% です。
% Compute the overall accuracy of the classifier
sum(diag(confdata))/sum(confdata(:))ans = 0.6150
% Compute the accuracy of the classifier at predicting % that there is a fault 1-sum(confdata(2:end,1))/sum(confdata(:))
ans = 0.9450
故障なしと予測され、故障が存在したケースについて調べます。まず、検証データで、実際の故障が閉塞故障でありながら、故障なしと予測されたケースを見つけます。
vData = data(cvp.test(1),:); b1 = (actualValue==2) & (predictedValue==0); fData = vData(b1,15:17)
fData=11×3 table
LeakFault BlockingFault BearingFault
_________ _____________ ____________
1e-09 0.77 0
1e-09 0.77 0
1e-09 0.71 0
1e-09 0.77 0
1e-09 0.77 0
1e-09 0.62 0
1e-09 0.77 0
1e-09 0.77 0
1e-09 0.71 0
8e-07 0.74 0
1e-09 0.74 0
検証データで、実際の故障がベアリング故障でありながら、故障なしとして予測されたケースを見つけます。
b2 = (actualValue==4) & (predictedValue==0); vData(b2,15:17)
ans = 0×3 empty table
故障なしと予測されたにもかかわらず故障が存在したケースを調べてみると、これは閉塞故障値の 0.77 がそのノミナル値 0.8 に近い場合、またはベアリング故障値の 6.6e-5 がそのノミナル値 0 に近い場合に生じることが明らかになります。閉塞故障値が小さいケースのスペクトルをプロットし、故障なしの状態と比較すると、スペクトルがよく似ているために検出が困難になることがわかります。分類器を再学習させて、閉塞値の 0.77 を故障なしの状態として含めると、故障検出器の性能は大幅に改善されます。あるいは、追加のポンプ測定値を使用してさらに多くの情報を提供することにより、小さな閉塞故障を検出できるように機能を改善することができます。
% Configure the ensemble to only read the flow and fault variable values ens.SelectedVariables = ["qOut_meas","LeakFault","BlockingFault","BearingFault"]; reset(ens) % Load the ensemble member data into memory data = gather(tall(ens));
Evaluating tall expression using the Parallel Pool 'local': - Pass 1 of 1: Completed in 38 sec Evaluation completed in 38 sec
% Look for the member that was incorrectly predicted, and % compute its power spectrum idx = ... data.BlockingFault == fData.BlockingFault(1) & ... data.LeakFault == fData.LeakFault(1) & ... data.BearingFault == fData.BearingFault(1); flow1 = data(idx,1); flow1 = flow1.qOut_meas{1}; [flow1P,flow1F] = pspectrum(flow1); % Look for a member that does not have any faults idx = ... data.BlockingFault == 0.8 & ... data.LeakFault == 1e-9 & ... data.BearingFault == 0; flow2 = data(idx,1); flow2 = flow2.qOut_meas{1}; [flow2P,flow2F] = pspectrum(flow2); % Plot the power spectra semilogx(... flow1F,pow2db(flow1P),... flow2F,pow2db(flow2P)); xlabel('Hz') legend('Small blocking fault','No fault')
まとめ
この例では Simulink モデルを使用して、往復ポンプの故障をモデル化し、さまざまな故障の組み合わせと重大度の下でモデルをシミュレートして、ポンプの出流量から状態インジケーターを抽出し、その状態インジケーターを使用してポンプの故障を検出するための分類器を学習させる方法を示しました。ここでは分類器を使用した故障検出の性能を調べ、小さな閉塞故障は故障なしの状態によく似ているため正確に検出できないことを確認しました。
サポート関数
function [flow,flowSpectrum,flowFrequencies,faultValues] = preprocess(data) % Helper function to preprocess the logged reciprocating pump data. % Remove the 1st 0.8 seconds of the flow signal tMin = seconds(0.8); flow = data.qOut_meas{1}; flow = flow(flow.Time >= tMin,:); flow.Time = flow.Time - flow.Time(1); % Ensure the flow is sampled at a uniform sample rate flow = retime(flow,'regular','linear','TimeStep',seconds(1e-3)); % Remove the mean from the flow and compute the flow spectrum fA = flow; fA.Data = fA.Data - mean(fA.Data); [flowSpectrum,flowFrequencies] = pspectrum(fA,'FrequencyLimits',[2 250]); % Find the values of the fault variables from the SimulationInput simin = data.SimulationInput{1}; vars = {simin.Variables.Name}; idx = strcmp(vars,'leak_cyl_area_WKSP'); LeakFault = simin.Variables(idx).Value; idx = strcmp(vars,'block_in_factor_WKSP'); BlockingFault = simin.Variables(idx).Value; idx = strcmp(vars,'bearing_fault_frict_WKSP'); BearingFault = simin.Variables(idx).Value; % Collect the fault values in a cell array faultValues = {... 'LeakFault', LeakFault, ... 'BlockingFault', BlockingFault, ... 'BearingFault', BearingFault}; end function ci = extractCI(flow,flowP,flowF) % Helper function to extract condition indicators from the flow signal % and spectrum. % Find the frequency of the peak magnitude in the power spectrum. pMax = max(flowP); fPeak = flowF(flowP==pMax); % Compute the power in the low frequency range 10-20 Hz. fRange = flowF >= 10 & flowF <= 20; pLow = sum(flowP(fRange)); % Compute the power in the mid frequency range 40-60 Hz. fRange = flowF >= 40 & flowF <= 60; pMid = sum(flowP(fRange)); % Compute the power in the high frequency range >100 Hz. fRange = flowF >= 100; pHigh = sum(flowP(fRange)); % Find the frequency of the spectral kurtosis peak [pKur,fKur] = pkurtosis(flow); pKur = fKur(pKur == max(pKur)); % Compute the flow cumulative sum range. csFlow = cumsum(flow.Data); csFlowRange = max(csFlow)-min(csFlow); % Collect the feature and feature values in a cell array. % Add flow statistic (mean, variance, etc.) and common signal % characteristics (rms, peak2peak, etc.) to the cell array. ci = {... 'qMean', mean(flow.Data), ... 'qVar', var(flow.Data), ... 'qSkewness', skewness(flow.Data), ... 'qKurtosis', kurtosis(flow.Data), ... 'qPeak2Peak', peak2peak(flow.Data), ... 'qCrest', peak2rms(flow.Data), ... 'qRMS', rms(flow.Data), ... 'qMAD', mad(flow.Data), ... 'qCSRange',csFlowRange, ... 'fPeak', fPeak, ... 'pLow', pLow, ... 'pMid', pMid, ... 'pHigh', pHigh, ... 'pKurtosis', pKur(1)}; end
参考
関連するトピック
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
Select web siteYou can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
