(57)【要約】  内燃機関の運転パラメータの検出値についての周波数成分分析を、機関の回転に同期して行う周波数成分分析装置が提供される。運転パラメータが所定時間間隔でサンプリングされ、サンプル値がデジタル値に変換される。検出値に含まれる複数の周波数成分に対応する、第１要素の強度及び第１要素に対して位相が９０度ずれた第２要素の強度が、所定数のサンプル値について算出される。第１要素強度及び第２要素強度を用いて複数周波数の周波数成分強度が機関の回転に同期して算出される。機関の回転速度が設定閾値以上であるときは、第１要素強度の積算値の一部及び第２要素強度の積算値の一部をそれぞれ前回算出値に置き換えて周波数成分強度の算出が行われる。
2010035634.tif 000002
【特許請求の範囲】
【請求項１】
内燃機関の運転パラメータの検出値についての周波数成分分析を、前記機関の回転に同期して行う周波数成分分析装置であって、
前記運転パラメータを所定時間間隔でサンプリングし、サンプル値をデジタル値に変換するサンプリング手段と、
前記検出値に含まれる複数の周波数成分に対応する、第１要素の強度及び前記第１要素に対して位相が９０度ずれた第２要素の強度を、所定数のサンプル値について算出する要素強度算出手段と、
前記第１要素強度及び前記第２要素強度を用いて前記複数周波数の周波数成分強度を前記機関の回転に同期して算出する周波数成分強度算出手段とを備え、
前記周波数成分強度算出手段は、前記機関の回転速度が設定閾値以上であるときは、前記第１要素強度の積算値の一部及び前記第２要素強度の積算値の一部をそれぞれ前回算出値に置き換えて前記周波数成分強度の算出を行う周波数成分分析装置。
【請求項２】
前記設定閾値は、前記所定数のサンプル値を得るのに要する時間に応じて設定される請求項１の周波数成分分析装置。
【請求項３】
前記設定閾値は、前記周波数成分強度の算出に必要とされる単位時間当たりの処理ステップ数に応じて設定される請求項１の周波数成分分析装置。
【請求項４】
前記周波数成分強度算出手段は、前記機関の回転に同期したトリガ信号の発生時期を中心とした所定期間内に得られたサンプル値を用いて前記周波数成分強度の算出を行う請求項１から３の何れか１項の周波数成分分析装置。
【請求項５】
内燃機関の運転パラメータの検出値についての周波数成分分析を、前記機関の回転に同期して行う周波数成分分析方法であって、
ａ）前記運転パラメータを所定時間間隔でサンプリングしてサンプル値をデジタル値に変換し、
ｂ）前記検出値に含まれる複数の周波数成分に対応する、第１要素の強度及び前記第１要素に対して位相が９０度ずれた第２要素の強度を、所定数のサンプル値について算出し、
ｃ）前記第１要素強度及び前記第２要素強度を用いて前記複数周波数の周波数成分強度を前記機関の回転に同期して算出する
ステップを備え、
前記機関の回転速度が設定閾値以上であるときは、前記第１要素強度の積算値の一部及び前記第２要素強度の積算値の一部をそれぞれ前回算出値に置き換えて前記周波数成分強度の算出が行われる周波数成分分析方法。
【請求項６】
前記設定閾値は、前記所定数のサンプル値を得るのに要する時間に応じて設定される請求項５の周波数成分方法。
【請求項７】
前記設定閾値は、前記周波数成分強度の算出に必要とされる単位時間当たりの処理ステップ数に応じて設定される請求項５の周波数成分分析方法。
【請求項８】
前記周波数成分強度算出手段は、前記機関の回転に同期したトリガ信号の発生時期を中心とした所定期間内に得られたサンプル値を用いて前記周波数成分強度の算出を行う請求項５から７の何れか１項の周波数成分分析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、内燃機関の運転パラメータの検出値についての周波数成分分析を、機関回転に同期して行う周波数成分分析装置に関する。
【背景技術】
【０００２】
特許文献１には、離散フーリエ変換（Discrete Fourier Transform）アルゴリズム（以下「ＤＦＴアルゴリズム」という）と、高速フーリエ変換（Fast Fourier Transform）アルゴリズム（以下「ＦＦＴアルゴリズム」という）とを用いて、内燃機関に装着されたノックセンサの出力信号の周波数成分分析を行う信号処理装置が示されている。ＦＦＴアルゴリズムは演算速度を高めるために特殊なアルゴリズムが採用されるため、必要な周波数成分（具体的にはノックセンサ出力信号の中心周波数成分）の強度が得られない場合がある。そこで、上記信号処理装置では、ノックセンサ出力信号の中心周波数成分（強度が最大の成分）についてはＤＦＴアルゴリズムを用いて強度が算出され、中心周波数成分以外の周波数成分についてはＦＦＴアルゴリズムを用いて強度が算出される。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平１１−３０３６７３号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ノッキングセンサ出力信号の周波数成分分析は、基本的に機関の回転に同期して行う必要があるため、機関回転速度が高い状態では、周波数成分分析の演算時間をより短縮することが必要となる。【０００５】
しかしながら、特許文献１に示された装置ではこの点が考慮されていないため、周波数成分分析によって得られる周波数成分強度を用いたノッキングの判別処理に確保できる時間が短くなり、正確な判定が困難となるおそれがある。【０００６】
本発明はこの点に着目してなされたものであり、検出パラメータの周波数成分分析を内燃機関の回転に同期して行う際に、周波数成分強度演算をより適切に実行し、特に機関回転速度が高い状態での演算速度を高めることができる周波数成分分析装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記目的を達成するため本発明は、内燃機関の運転パラメータの検出値（ＶＫＮＫ）についての周波数成分分析を、前記機関の回転に同期して行う周波数成分分析装置を提供する。この周波数成分分析装置は、サンプリング手段と、要素強度算出手段と、周波数成分強度算出手段とを備えている。前記サンプリング手段は、前記運転パラメータを所定時間（ＴＳＭＰ）間隔でサンプリングし、サンプル値をデジタル値に変換する。前記要素強度算出手段は、前記検出値（ＶＫＮＫ）に含まれる複数の周波数成分に対応する、第１要素の強度（ＤＭＦＴＳ）及び前記第１要素に対して位相が９０度ずれた第２要素の強度（ＤＭＦＴＣ）を、所定数（ＮＤ）のサンプル値について算出する。周波数成分強度算出手段は、前記第１要素強度（ＤＭＦＴＳ）及び前記第２要素強度（ＤＭＦＴＣ）を用いて前記複数周波数の周波数成分強度（ＳＴＦＴ）を前記機関の回転に同期して算出する。さらに前記周波数成分強度算出手段は、前記機関の回転速度（ＮＥ）が設定閾値（ＮＥＴＨ）以上であるときは、前記第１要素強度の積算値の一部及び前記第２要素強度の積算値の一部をそれぞれ前回算出値（ＴＲＭＳ１，ＴＲＭＣ１）に置き換えて前記周波数成分強度（ＳＴＦＴ）の算出を行う。【０００８】
この構成によれば、内燃機関の運転パラメータが所定時間間隔でサンプリングされ、サンプル値がデジタル値に変換される。得られたデジタル値に基づいて複数の周波数成分に対応する、第１要素の強度及び第１要素に対して位相が９０度ずれた第２要素の強度が、所定数のサンプル値について算出され、第１要素強度及び第２要素強度を用いて、複数周波数の周波数成分強度が機関の回転に同期して算出される。機関の回転速度が設定閾値以上であるときは、第１要素強度の積算値の一部及び第２要素強度の積算値の一部をそれぞれ前回算出値に置き換えて周波数成分強度の算出が行われる。したがって、周波数成分強度を算出するときに、第１及び第２要素強度の積算演算の一部を再度行う必要がなくなり、演算速度を高めることができる。その結果、算出された周波数成分強度を用いた後処理の時間を十分に確保することが可能となる。【０００９】
前記設定閾値（ＮＥＴＨ）は、前記所定数（ＮＤ）のサンプル値を得るのに要する時間（ＴＳ）に応じて設定されることが望ましい。【００１０】
この構成によれば、機関回転速度の設定閾値は、所定数のサンプル値を得るのに要する時間に応じて設定される。具体的には、機関回転速度が高くなって周波数成分強度の演算周期が、所定数のサンプル値を得るのに要する時間より短くなると、周波数成分強度の算出に適用する第１及び第２要素の強度の前回算出値と今回算出値の一部が重複するので、その重複している強度について前回算出値を使用することが可能となる。【００１１】
前記設定閾値（ＮＥ）は、前記周波数成分強度（ＳＴＦＴ）の算出に必要とされる単位時間当たりの処理ステップ数（ＮＳ）に応じて設定されることが望ましい。【００１２】
この構成によれば、機関回転速度の設定閾値は、周波数成分強度の算出に必要とされる単位時間当たりの処理ステップ数に応じて設定される。上述したように、周波数成分強度の演算周期が、所定数のサンプル値を得るのに要する時間より短くなると、周波数成分強度の算出に適用する第１及び第２要素の強度の前回算出値と今回算出値の一部が重複するので、その重複している強度について前回算出値を使用することが可能となる。しかし、実際の演算装置では、周波数成分強度の演算に要素強度の前回算出値を適用する処理を行うことにより処理ステップ数が増加し、かえって演算速度の低下を招く場合がある。そこで、前回算出値を適用する処理を行っても単位時間当たりのステップ数が増加しなくなる機関回転速度を設定閾値とすることにより、演算装置の性能（特性）に応じた適切な強度演算を行うことができる。【００１３】
前記周波数成分強度算出手段は、前記機関の回転に同期したトリガ信号（ＣＲＫ割り込み）の発生時期を中心とした所定期間（ＴＳ）内に得られるサンプル値を用いて前記周波数成分強度（ＳＴＦＴ）の算出を行うことが望ましい。【００１４】
この構成によれば、機関の回転に同期したトリガ信号の発生時期を中心とした所定期間内にサンプリングされたデジタル値を用いて周波数成分強度の算出が行われるので、対象とする機関の回転角度を中心とした周波数成分解析を行うことができる。その結果、トリガ信号が発生した時点を始点とする所定期間内またはトリガ信号が発生した時点を終点とする所定期間内でサンプリングされたデジタル値を用いる手法より、機関回転変動の影響を受け難くなるという効果が得られる。
【図面の簡単な説明】
【００１５】
【図１】本発明の一実施形態にかかる内燃機関及びその制御装置の構成を示す図である。
【図２】ノックセンサ出力のサンプリング及び周波数成分分析を説明するための図である。
【図３】周波数成分強度（ＳＴＦＴ）を算出する手順及びタイミング関係を説明するためのタイムチャートである。
【図４】周波数成分強度（ＳＴＦＴ）の算出期間（ＴＳ）と、ＣＲＫ割り込みの発生周期（ＣＲＭＥ）との関係を説明するためのタイムチャートである。
【図５】機関回転数（ＮＥ）と単位時間当たりの演算ステップ数との関係を示す図である。
【図６】周波数成分強度（ＳＴＦＴ）を算出する手順を説明するためのフローチャートである。
【図７】周波数成分分析による得られる強度データをスペクトル時系列マップ及び二値化スペクトル時系列マップとして示した図である。
【図８】吸気弁の着座ノイズの二値化スペクトル時系列マップを示す図である。
【図９】ノイズ除去処理を説明するために二値化スペクトル時系列マップを示す図である。
【図１０】マスタパターンマップを用いたノッキング判定を説明するために二値化スペクトル時系列マップを示す図である。
【図１１】重み付けマップの一例を示す図である。
【図１２】ノッキング判定処理のフローチャートである。
【図１３】図１２の処理で参照されるマップを示す図である。
【図１４】図１２に示す処理で実行される二値化データマップ算出処理のフローチャートである。
【図１５】図１４に示す処理で実行される二値化処理のフローチャートである。
【図１６】図１５の処理で参照されるマップを示す図である。
【図１７】図１２に示す処理で実行されるノイズ除去処理のフローチャートである。
【図１８】図１２に示す処理で実行される適合率算出処理のフローチャートである。
【図１９】図１８の処理で参照されるマップを説明するための図である。
【図２０】図１２に示す処理で実行されるノイズ学習処理のフローチャートである。
【図２１】図２０に示す処理で実行されるノイズマップ更新処理のフローチャートである。
【図２２】ノックセンサにより検出される振動の周波数成分強度を示す図である。
【発明を実施するための形態】
【００１６】
以下本発明の実施の形態を図面を参照して説明する。
図１は、本発明の一実施形態にかかる内燃機関（以下「エンジン」という）及びその制御装置の全体構成図であり、例えば４気筒のエンジン１の吸気管２の途中にはスロットル弁３が配されている。スロットル弁３にはスロットル弁開度ＴＨを検出するスロットル弁開度センサ４が連結されており、センサ４の検出信号は、電子制御ユニット（以下「ＥＣＵ」という）５に供給される。
【００１７】
燃料噴射弁６はエンジン１とスロットル弁３との間かつ吸気管２の図示しない吸気弁の少し上流側に各気筒毎に設けられており、各噴射弁は図示しない燃料ポンプに接続されていると共にＥＣＵ５に電気的に接続されて当該ＥＣＵ５からの信号により燃料噴射弁６の開弁時間が制御される。エンジン１の各気筒には、点火プラグ７が設けられており、点火プラグ７はＥＣＵ５に接続されている。ＥＣＵ５は、点火プラグ７に点火信号を供給する。
【００１８】
スロットル弁３の下流側には吸気圧ＰＢＡを検出する吸気圧センサ８、及び吸気温ＴＡを検出する吸気温センサ９が設けられている。エンジン１の本体には、エンジン冷却水温ＴＷを検出する冷却水温センサ１０及び非共振型のノックセンサ１１が装着されている。センサ８〜１１の検出信号は、ＥＣＵ５に供給される。ノックセンサ１１としては、例えば５ｋＨｚから２５ｋＨｚまでの周波数帯域の振動を検出可能なものが使用される。
【００１９】
吸気管２のスロットル弁３の上流側には、吸入空気流量ＧＡを検出する吸入空気流量センサ１３が設けられており、その検出信号はＥＣＵ５に供給される。
ＥＣＵ５には、エンジン１のクランク軸（図示せず）の回転角度を検出するクランク角度位置センサ１２が接続されており、クランク軸の回転角度に応じた信号がＥＣＵ５に供給される。クランク角度位置センサ１２は、エンジン１の特定の気筒の所定クランク角度位置でパルス（以下「ＣＹＬパルス」という）を出力する気筒判別センサ、各気筒の吸入行程開始時の上死点（ＴＤＣ）に関し所定クランク角度前のクランク角度位置で（４気筒エンジンではクランク角１８０度毎に）ＴＤＣパルスを出力するＴＤＣセンサ及びＴＤＣパルスより短い一定クランク角周期（例えば６度周期）で１パルス（以下「ＣＲＫパルス」という）を発生するＣＲＫセンサから成り、ＣＹＬパルス、ＴＤＣパルス及びＣＲＫパルスがＥＣＵ５に供給される。これらのパルスは、燃料噴射時期、点火時期等の各種タイミング制御、エンジン回転数（エンジン回転速度）ＮＥの検出に使用される。
【００２０】
エンジン１は、吸気弁（図示せず）の弁リフト量及び開角（開弁期間）を連続的に変更する第１弁作動特性可変機構と、吸気弁を駆動するカムの、クランク軸回転角度を基準とした作動位相を連続的に変更する第２弁作動特性可変機構とを有する弁作動特性可変装置２０を備えている。ＥＣＵ５は、弁作動特性可変装置２０にリフト量制御信号及び作動位相制御信号を供給し、吸気弁の作動制御を行う。第１及び第２弁作動特性可変機構の構成は、それぞれ例えば特開２００８−２５４１８号公報及び特開２０００−２２７０１３号公報に示されている。
【００２１】
ＥＣＵ５は、各種センサからの入力信号波形を整形し、電圧レベルを所定レベルに修正し、アナログ信号値をデジタル信号値に変換する等の機能を有する入力回路、中央演算処理ユニット（以下「ＣＰＵ」という）、該ＣＰＵで実行される各種演算プログラム及び演算結果等を記憶する記憶回路（メモリ）、燃料噴射弁６及び点火プラグ７に駆動信号を供給する出力回路等から構成される。
【００２２】
本実施形態では、ノックセンサ１１の出力信号の周波数成分分析を行い、その分析の結果得られる周波数成分強度に基づいてノッキング判定を行う。そこで先ず、周波数成分分析の概要を説明する。図２（ａ）は、ノックセンサ１１の出力信号波形を示し、図２（ｂ）は、同図（ａ）の期間ＴＳの波形を拡大して示す図である。本実施形態では、サンプリング周期ＴＳＭＰを２０マイクロ秒とし、連続して検出される５０個のデータを対象として、離散フーリエ変換（ＤＦＴ）による周波数成分分析を行う。その周波数成分分析の結果が図２（ｃ）に示されている。図２（ｃ）の縦軸は周波数成分強度ＳＴＦＴであり、本実施形態では、５ｋＨｚから２５ｋＨｚまでの周波数帯域における１ｋＨｚ毎の周波数（５，６，７，…，２４，２５ｋＨｚ）に対応する周波数成分強度ＳＴＦＴが、クランク角度６度毎に（エンジン１のクランク軸が６度回転する毎に）算出される。
【００２３】
図３は上記周波数成分分析を説明するためのタイムチャートであり、図３（ａ）はノックセンサ１１の出力信号を２０マイクロ秒毎にサンプリングすることにより得られる５０個の検出データＶＫＮＫ（デジタル値）が順次格納されるメモリのアドレスに対応するアドレス番号を示しており、同図（ｂ）は検出データＶＫＮＫを用いて算出されるサイン波成分強度ＤＭＦＴＳ及びコサイン波成分強度ＤＭＦＴＣが、検出データＶＫＮＫに対応してメモリに格納される状態を示している。アドレス番号は、同図（ａ）と同一であるが、１つの検出データＶＫＮＫのサンプル値に対応して、周波数５ｋＨｚから２５ｋＨｚまでの２１個の周波数に対応するサイン波成分強度ＤＭＦＴＳ及びコサイン波成分強度ＤＭＦＴＣが１ミリ秒毎に算出され、メモリに格納される。
【００２４】
図３（ｃ）は、クランク角度位置センサ１２から出力されるＣＲＫパルスを示し、本実施形態ではＣＲＫパルスの立ち下がり時期が演算実行時期の基準（以下「ＣＲＫ割り込み」という）として使用される。本実施形態では、ＣＲＫ割り込みが発生したタイミングを記憶しておき、そのタイミングを中心としたサンプリング期間ＴＳ中に得られる５０個の検出データに対応するサイン波成分強度ＤＭＦＴＳ及びコサイン波成分強度ＤＭＦＴＣを用いて各周波数（５〜２５ｋＨｚ）の周波数成分強度ＳＴＦＴ(j,i)が算出される（図３（ｄ））。ここで「ｊ」は周波数を示すインデクスパラメータであり（以下「周波数インデクス」という）、ｊ＝０，１，２，…，２０が周波数５，６，７，…，２５ｋＨｚに対応する。「ｉ」は、ＣＲＫ割り込みが発生するクランク角度ＣＡ（ピストンが上死点に位置するときの角度位置を０度とする）を示すインデクスパラメータであり（以下「クランク角インデクス」という）、ｉ＝０，１，２，…，１４が、クランク角度６，１２，１８，…，９０度に対応する。
【００２５】
図３に示す例では、インデクスパラメータｍ（＝１〜５０）が、周波数成分強度ＳＴＦＴ(j,i)の算出に適用されるデータを示す。最初のＣＲＫ割り込みタイミングでは、アドレス番号「５」の検出データ（ｍ＝２５に対応する検出データ）を中心とした５０個のデータを用いて周波数成分強度ＳＴＦＴ(j,i)が算出され、次のＣＲＫ割り込みタイミングでは、アドレス番号「４５」の検出データ（ｍ＝２５に対応する検出データ）を中心とした５０個のデータを用いて周波数成分強度ＳＴＦＴ(j,i+1)が算出される。なお、周波数成分強度ＳＴＦＴ(j,i)は、相対強度を示す無次元量である。
【００２６】
周波数成分強度ＳＴＦＴ(j,i)の算出に適用されるデータを、ｍ＝２５に対応する検出データを中心とした５０個のデータ、換言すれば、ＣＲＫ割り込みの発生タイミングを中心としたサンプリング期間ＴＳ中に得られる５０個の検出データとすることにより、対象とするクランク角度を中心とした周波数成分解析を行うことができる。その結果、ＣＲＫ割り込みが発生した時点を始点とするサンプリング期間内またはＣＲＫ割り込みが発生した時点を終点とするサンプリング期間内でサンプリングされた検出データを用いる手法より、エンジン回転変動の影響を受け難くなるという効果が得られる。
【００２７】
図４は、ＣＲＫ割り込みの発生周期ＣＲＭＥと、周波数成分強度ＳＴＦＴの算出に適用される検出データのサンプリング期間（「ＳＴＦＴ算出」と表示）ＴＳとの関係を示す図である。図４（ａ）は、エンジン回転数ＮＥが１０００ｒｐｍの例を示し、サンプリング期間ＴＳと割り込み発生周期ＣＲＭＥとが一致している。したがって、エンジン回転数ＮＥが１０００ｒｐｍを超えると、サンプリング期間ＴＳが一部重複する。図４（ｂ）は、エンジン回転数ＮＥが２０００ｒｐｍより若干高い例を示し、サンプリング期間ＴＳの重複部分が長くなっている。
【００２８】
本実施形態では、このようにサンプリング期間ＴＳが重複した状態においては、重複部分に関わるサイン波成分強度ＤＭＦＴＳの積算値及びコサイン波成分強度ＤＭＦＴＣの積算値として前回算出値を用いることにより、周波数成分強度ＳＴＦＴの算出に要する演算時間を短縮するようにしている。
【００２９】
ただし、前回算出値を使用する演算は、使用可能な重複データ数の比率が小さいときは、演算装置の特性に依存して単位時間当たりの処理ステップ数が増加する。図５の示す実線Ｌ１１は、本実施形態で使用するＣＰＵにおける、前回算出値を使用した場合の単位時間当たりの処理ステップ数ＮＳと、エンジン回転数ＮＥとの関係を示す図であり、回転数ＮＥＳＡＴ（以下「ステップ数飽和回転数ＮＥＳＡＴ」という）に達するまでは、処理ステップ数ＮＳは増加し、ＮＥ≧ＮＥＳＡＴの範囲では一定値ＮＳＳＡＴとなる。よって本実施形態では、閾値ＮＥＴＨをステップ数飽和回転数ＮＥＳＡＴに設定し、エンジン回転数ＮＥが閾値ＮＥＴＨ以上であるときは、前回算出値を使用する演算を行うようにしている。
【００３０】
図５に示す破線Ｌ１２は、ＦＦＴ（高速フーリエ変換）アルゴリズムを使用した場合の、エンジン回転数ＮＥと処理ステップ数ＮＳとの関係を示す。図５から明らかなようにＦＦＴアルゴリズムを用いる場合には、エンジン回転数ＮＥの増加に対して処理ステップ数ＮＳが直線的に増加するのに対し、ＤＦＴアルゴリズムを使用し、かつエンジン回転数ＮＥの高い領域（ＮＥ＞ＮＥＳＡＴ）で前回算出値を用いる演算を行うことにより、処理ステップ数ＮＳの増加を抑制することが可能となる。
【００３１】
図６は周波数成分強度ＳＴＦＴ算出の手順を示すフローチャートである。ステップＳ１では、検出データＶＫＮＫを取得し、ステップＳ２では、下記式（１）及び（２）によりサイン波成分強度ＤＭＦＴＳ(j,k)及びコサイン波成分強度ＤＭＦＴＣ(j,k)を算出する。ここでインデクスパラメータｋは、図３（ｂ）に示すアドレス番号であり、「１」からＮＤ（本実施形態では「５０」）までの値をとる。またΔｔは、（サンプリング周期×ＮＤ）に相当し、本実施形態では１ミリ秒である。
ＤＭＦＴＳ(j,k)＝ＶＫＮＫ(k)×
ｓｉｎ{２π×（ｊ＋５）×１０００×Δｔ×ｋ／ＮＤ}  （１）
ＤＭＦＴＣ(j,k)＝ＶＫＮＫ(k)×
ｃｏｓ{２π×（ｊ＋５）×１０００×Δｔ×ｋ／ＮＤ}  （２）
【００３２】
ステップＳ１及びＳ２の処理は、サンプリング周期（２０マイクロ秒）毎に実行される。
【００３３】
ステップＳ３では、ＣＲＫ割り込みが発生したか否かを判別し、発生していないときはステップＳ１に戻り、ＣＲＫ割り込みが発生したタイミングで、エンジン回転数ＮＥが閾値ＮＥＴＨ以上であるか否かを判別する。ＮＥ＜ＮＥＴＨであるときは、前回使用データを使用しない下記第１式（Ａ）により周波数成分強度ＳＴＦＴを算出する一方、ＮＥ≧ＮＥＴＨであるときは、前回使用データを使用する下記第２式（Ｂ）により周波数成分強度ＳＴＦＴ(j,i)を算出する。インデクスパラメータｍは、ＣＲＫ割り込みが発生したときのアドレス番号ｋが「２５」となるように設定される修正アドレス番号である（図３（ｄ）参照）。インデクスパラメータｍｐは、前回算出時に適用された修正アドレス番号であり、ｍｘは前回境界インデクスであり、ｍｙは今回境界インデクスであり、それぞれ下記式（３）及び（４）で与えられる。
【００３４】
【数１】
2010035634.tif 000003
ｍｘ＝ＮＤ−ＮＶ＋１                        （３）
ｍｙ＝ＮＶ＋１                              （４）
ここでＮＶは前回のＳＴＦＴ算出に使用したデータのうち、今回のＳＴＦＴ算出に使用できるデータの数であり、データ数ＮＶはエンジン回転数ＮＥが高くなるほど増加する。
【００３５】
第２式（Ｂ）の第１項ＴＲＭＳ１及び第３項ＴＲＭＣ１が、積算値の前回算出値に相当し、第２項ＴＲＭＳ２及び第３項ＴＲＭＣ２が新たに算出する積算値に相当する。
【００３６】
このようにサイン波成分強度の積算値の一部及びコサイン波成分強度の積算値一部を、それぞれ前回算出値（ＴＲＭＳ１，ＴＲＭＣ１）に置き換えて、周波数成分強度ＳＴＦＴ(j,i)の算出を行うことにより、周波数成分強度ＳＴＦＴを算出するときに、サイン波成分強度ＤＭＦＴＳ及びコサイン波成分強度ＤＭＦＴＣの積算演算の一部を再度行う必要がなくなり、演算速度を高めることができる。その結果、算出された周波数成分強度を用いた後処理の時間を十分に確保することが可能となる。
【００３７】
次に図７〜図１１を参照して、周波数成分強度ＳＴＦＴを用いたノッキング判定手法の概要を説明する。
上述したように周波数成分強度ＳＴＦＴは、時系列データとして図７（ａ）に示すように２次元マトリクス（以下「スペクトル時系列マップ」という）として算出される。スペクトル時系列マップは縦方向が周波数ｆ［ｋＨｚ］であり、横方向がクランク角度（燃焼行程が開始する上死点後のクランク角度）ＣＡ［ｄｅｇ］である。ここで、縦方向及び横方向のインデクスとしてそれぞれ周波数インデクスｊ（ｊ＝０〜２０）及びクランク角インデクスｉ（ｉ＝０〜１４）を用い、スペクトル時系列マップの要素を強度パラメータＫＭＡＰ(j,i)と表示する（強度パラメータＫＭＡＰは、本実施形態では周波数成分強度ＳＴＦＴ(j,i)と同じパラメータであり、スペクトル時系列マップ上の値であることからＫＭＡＰというラベルを使用する）。例えば強度パラメータＫＭＡＰ(0,0)は、左下端の「３４」、強度パラメータＫＭＡＰ(0,14)は右下端の「３１」、強度パラメータＫＭＡＰ(20,0)は左上端の「５６」、強度パラメータＫＭＡＰ(20,14)は右上端の「３０」に相当する。
【００３８】
図７（ａ）のスペクトル時系列マップを二値化すると図７（ｂ）の二値化スペクトル時系列マップが得られる。図７（ｂ）のマップの要素である二値化強度パラメータＮＫＭＡＰ(j,i)は、強度パラメータＫＭＡＰ(j,i)が「５０」以上であるとき「１」、「５０」未満であるとき「０」をとる。
【００３９】
図７（ｂ）に示す二値化スペクトル時系列マップは、図２２に破線で示すノッキングが発生したときに得られる周波数成分強度に対応するマップである。一方図２２に実線で示す吸気弁の着座ノイズ（エンジン１の吸気弁が全閉位置に到達したときに発生する振動に起因するノイズ）に対応する二値化スペクトル時系列マップは、図８に示すようになり、図７（ｂ）に示すノッキング発生時のマップと明確に相違するため、着座ノイズをノッキング発生と誤判定することがなく、正確なノッキング判定を行うことができる。
【００４０】
本実施形態ではさらに、ノッキングが発生していないと判定されたときの二値化スペクトル時系列マップに基づいて、ノイズに対応する二値化スペクトル時系列マップ（以下単に「ノイズマップ」という）を学習演算により生成し、ノイズマップ上のノイズ学習値を、検出データから算出される二値化スペクトル時系列マップ上の対応する値から差し引くことにより、ノイズの影響を除去するようにしている。
【００４１】
図９（ａ）は、学習演算により得られたノイズマップの一例を示す図であり、クランク角度６度のタイミングで６〜２５ｋＨｚの帯域でノイズ成分が表れている。図９（ｂ）は、図７（ｂ）に示す二値化スペクトル時系列マップから図９（ａ）に示すノイズ成分を減算することによりノイズ除去処理を行った後の二値化スペクトル時系列マップを示す。このようにノイズ除去処理を行うことにより、ノイズの影響を除き、より正確な判定が可能となる。
【００４２】
また本実施形態では、検出データから得られる二値化スペクトル時系列マップを、図１０（ａ）に示すマスタパターンマップと比較することにより、ノッキングが発生したか否かを判定するようにしている。マスタパターンマップは、ノッキングが発生したときに得られる典型的な二値化スペクトル時系列マップに相当するものである。
【００４３】
上記比較は具体的には以下のように行われる。二値化スペクトル時系列マップ上の二値化強度パラメータＮＫＭＡＰ(j,i)と、マスタパターンマップ上のマスタパラメータＭＭＡＰ(j,i)との積を積算することにより、強度積算値ＳＵＭＫを算出し、マスタパラメータＭＭＡＰ(j,i)そのものの積算値である基準積算値ＳＵＭＭを算出し、基準積算値ＳＵＭＭに対する強度積算値ＳＵＭＫの比率として適合率ＰＦＩＴ（＝ＳＵＭＫ／ＳＵＭＭ）を算出する。そして、適合率ＰＦＩＴが判定閾値ＳＬＶＬを超えるとノッキングが発生したと判定する。
【００４４】
図１０（ｂ）は、図９（ｂ）に示す二値化強度パラメータＮＫＭＡＰ(j,i)と、図１０（ａ）に示すマスタパターンマップ上のマスタパラメータＭＭＡＰ(j,i)との積（ＮＫＭＡＰ(j,i)×ＭＭＡＰ(j,i)）のマップを示す。この例では、図１０（ｂ）の積マップは、図９（ｂ）に示すもとのマップと全く同一となり、適合率ＰＦＩＴは、「０．８６３」となる。
【００４５】
なお、後述する実際のノッキング判定処理においては、上記強度積算値ＳＵＭＫ及び基準積算値ＳＵＭＭを算出する際にエンジン運転状態に応じた重み付けを行うようにしている。この重み付けを行うための重み付けパラメータＷＭＡＰ(j,i)が設定されている重み付けマップの一例を図１１に示す。図１１に示す重み付けマップは、周波数６ｋＨｚの成分についてクランク角１２度から７２度までのパラメータ値、周波数１０ｋＨｚの成分についてクランク角３６度から６６度までのパラメータ値、並びに周波数１１ｋＨｚ、１３ｋＨｚ、１５ｋＨｚ、及び２０ｋＨｚの成分についてクランク角１２度から３０度までのパラメータ値に対して重みが付けられるように設定されている。
【００４６】
図１１に示すような重み付けマップによる重み付けは、エンジン運転状態に依存して二値化スペクトル時系列マップの周波数に対する特性が変化することを補償するために行われるものである。エンジン運転状態に応じて重み付けを行うことにより、エンジン運転状態の変化に拘わらず正確な判定を行うことが可能となる。
【００４７】
図１２は、上述した手法によりノッキング判定を行う処理のフローチャートであり、この処理はＥＣＵ５のＣＰＵでＴＤＣパルスの発生に同期して実行される。
【００４８】
ステップＳ１１では、図１４に示す二値化データマップ算出処理を実行し、上述した二値化スペクトル時系列マップの算出を行う。ステップＳ１２では、図１７に示すノイズ除去処理を実行し、図９（ａ）に例示するようなノイズマップを用いてノイズ成分を除去する処理を行う。
【００４９】
ステップＳ１３では、図１８に示す適合率算出処理を実行し、ノイズ成分が除去された二値化スペクトル時系列マップと、マスタパターンマップとを用いて適合率ＰＦＩＴを算出する。
【００５０】
ステップＳ１４では、エンジン回転数ＮＥ及び吸気圧ＰＢＡに応じて例えば図１３に示すように設定されたＳＬＶＬマップを検索し、判定閾値ＳＬＶＬを算出する。ＳＬＶＬマップに設定されている格子点以外のエンジン回転数ＮＥ及び吸気圧ＰＢＡについては、補間演算により判定閾値ＳＬＶＬを算出する。
【００５１】
ステップＳ１５では、ステップＳ１３で算出される適合率ＰＦＩＴが判定閾値ＳＬＶＬより大きいか否かを判別し、その答が肯定（ＹＥＳ）であるときは、ノッキングが発生したと判定し、ノッキングフラグＦＫＮＯＣＫを「１」に設定する（ステップＳ１６）。
【００５２】
ステップＳ１５でＰＦＩＴ≦ＳＬＶＬであるときは、ノッキングは発生していないと判定し、ノッキングフラグＦＫＮＯＣＫを「０」に設定する（ステップＳ１７）。次いで図２０に示すノイズ学習処理を実行し、ノイズマップ（図９（ａ）参照）の更新を行う。
【００５３】
図１４は、図１２のステップＳ１１で実行される二値化データマップ算出処理のフローチャートである。
ステップＳ２１では、クランク角インデクスｉ及び周波数インデクスｊをいずれも「０」に初期化する。ステップＳ２２では、周波数インデクスｊが周波数データ数ＪＮ（本実施形態では２１）から「１」を減算した値より大きいか否かを判別する。最初はこの答は否定（ＮＯ）であるので、ステップＳ２３に進み、クランク角インデクスｉがクランク角データ数ＩＮ（本実施形態では１５）から「１」を減算した値より大きいか否かを判別する。
【００５４】
最初はステップＳ２３の答も否定（ＮＯ）であるので、強度パラメータＫＭＡＰ(j,i)を、メモリに格納されている周波数成分強度ＳＴＦＴ(j,i)に設定し（ステップＳ２５）、次いで図１５に示す二値化処理を実行する（ステップＳ２６）。ステップＳ２７では、クランク角インデクスｉを「１」だけインクリメントする。
【００５５】
図１５のステップＳ３１では、エンジン回転数ＮＥ及び吸気圧ＰＢＡに応じて図１６に示すＢＬＶＬマップを検索し、二値化閾値ＢＬＶＬを算出する。ＢＬＶＬマップは、エンジン回転数ＮＥが増加するほど二値化閾値ＢＬＶＬが増加し、かつ吸気圧ＰＢＡが増加するほど二値化閾値ＢＬＶＬが増加するように設定されている。ラインＬ１で示される設定値は、第１所定吸気圧ＰＢＡ１（例えば５３ｋＰａ（４００ｍｍＨｇ））から第２所定吸気圧ＰＢＡ２（例えば８０ｋＰａ（６００ｍｍＨｇ））の範囲で適用される。ラインＬ２及びＬ３は、それぞれ第３所定吸気圧ＰＢＡ３（例えば９３ｋＰａ（７００ｍｍＨｇ））及び第４所定吸気圧ＰＢＡ４（例えば１０７ｋＰａ（８００ｍｍＨｇ））に対応する。
【００５６】
ステップＳ３２では、強度パラメータＫＭＡＰ(j,i)が二値化閾値ＢＬＶＬより大きいか否かを判別し、その答が肯定（ＹＥＳ）であるときは、二値化強度パラメータＮＫＭＡＰ(j,i)を「１」に設定する（ステップＳ３３）。一方、ステップＳ３２でＫＭＡＰ(j,i)≦ＢＬＶＬであるときは、二値化強度パラメータＮＫＭＡＰ(j,i)を「０」に設定する（ステップＳ３４）。
【００５７】
図１４に戻り、ステップＳ２３の答が否定（ＮＯ）である間はステップＳ２４〜Ｓ２７を繰り返し実行し、クランク角インデクスｉが（ＩＮ−１）を超えると、ステップＳ２８に進み、クランク角インデクスｉを「０」に戻すとともに、周波数インデクスｊを「１」だけインクリメントし、ステップＳ２２に戻る。
【００５８】
ステップＳ２２の答が否定（ＮＯ）である間は、ステップＳ２３〜Ｓ２８を繰り返し実行し、周波数インデクスｊが（ＪＮ−１）を超えると、本処理を終了する。
【００５９】
図１４の処理により、図７（ａ）に示すスペクトル時系列マップ上の値を示す強度パラメータＫＭＡＰ(j,i)が、メモリに格納されている、周波数成分分析の結果として得られた周波数成分強度ＳＴＦＴ(j,i)に設定されるとともに、強度パラメータＫＭＡＰ(j,i)の二値化が行われ、二値化強度パラメータＮＫＭＡＰ(j,i)が算出される。すなわち、図７（ｂ）に示す二値化スペクトル時系列マップが生成される。
【００６０】
図１７は、図１２のステップＳ１２で実行されるノイズ除去処理のフローチャートである。
ステップＳ４１では、クランク角インデクスｉ及び周波数インデクスｊをともに「０」に初期化する。ステップＳ４２では、周波数インデクスｊが周波数データ数ＪＮから「１」を減算した値より大きいか否かを判別する。最初はこの答は否定（ＮＯ）であるので、ステップＳ４３に進み、クランク角インデクスｉがクランク角データ数ＩＮから「１」を減算した値より大きいか否かを判別する。
【００６１】
最初はステップＳ４３の答も否定（ＮＯ）であるので、下記式（６）により二値化強度パラメータＮＫＭＡＰ(j,i)を補正し、補正二値化強度パラメータＪＫＭＡＰ(j,i)を算出する。式（６）のＮＮＭＡＰ(j,i)は、学習処理により更新されるノイズマップ上の二値化ノイズパラメータである。
ＪＫＭＡＰ(j,i)＝ＮＫＭＡＰ(j,i)−ＮＮＭＡＰ(j,i)          （６）
【００６２】
ステップＳ４５では、補正二値化強度パラメータＪＫＭＡＰ(j,i)が負の値であるか否かを判別し、その答が否定（ＮＯ）であるときは直ちにステップＳ４７に進む。ＪＫＭＡＰ(j,i)＜０であるときは、ＪＫＭＡＰ(j,i)を「０」に設定し（ステップＳ４６）、ステップＳ４７に進む。
【００６３】
ステップＳ４７ではクランク角インデクスｉを「１」だけインクリメントし、ステップＳ４３に戻る。ステップＳ４３の答が否定（ＮＯ）である間はステップＳ４４〜Ｓ４７を繰り返し実行し、クランク角インデクスｉが（ＩＮ−１）を超えると、ステップＳ４８に進み、クランク角インデクスｉを「０」に戻すとともに、周波数インデクスｊを「１」だけインクリメントし、ステップＳ４２に戻る。ステップＳ４２の答が否定（ＮＯ）である間は、ステップＳ４３〜Ｓ４８を繰り返し実行し、周波数インデクスｊが（ＪＮ−１）を超えると、本処理を終了する。
【００６４】
図１７の処理により、ノイズが除去された補正二値化強度パラメータＪＫＭＡＰ(j,i)が得られる（図９（ｂ）参照）。
【００６５】
図１８は、図１２のステップＳ１３で実行される適合率算出処理のフローチャートである。
ステップＳ５１では、エンジン回転数ＮＥ及び吸気圧ＰＢＡに応じてマスタパターンマップを選択し、ステップＳ５２では、エンジン回転数ＮＥ及び吸気圧ＰＢＡに応じて重み付けマップを選択する。重み付けマップは、エンジン運転状態に依存して二値化スペクトル時系列マップの周波数に対する特性が変化することを補償するために設けられている。エンジン回転数ＮＥまたは吸気圧ＰＢＡ（エンジン負荷）が変化すると、燃焼室内の温度が変化し、二値化スペクトル時系列マップが変化する。したがって、エンジン回転数ＮＥ及び吸気圧ＰＢＡに応じてマスタパターンマップ及び重み付けマップを選択することにより、エンジン運転状態の変化に拘わらず正確な判定を行うことが可能となる。
【００６６】
本実施形態では、図１９に示すようにエンジン回転数ＮＥ及び吸気圧ＰＢＡにより定義される９個のエンジン運転領域に対応して、９個のマスタパターンマップ及び９個の重み付けマップが予め設定されており、ステップＳ５１では９個のマスタパターンマップうちの１つが選択され、ステップＳ５２では、９個の重み付けマップのうちの１つが選択される。図１９において低回転領域は例えばエンジン回転数ＮＥが２０００ｒｐｍ以下の領域とされ、中回転領域は２０００ｒｐｍから４０００ｒｐｍまでの領域とされ、高回転領域は４０００ｒｐｍを超える領域とされる。また低負荷領域は例えば吸気圧ＰＢＡが６７ｋＰａ（５００ｍｍＨｇ）以下の領域とされ、中負荷領域は６７ｋＰａから９３ｋＰａ（７００ｍｍＨｇ）までの領域とされ、高負荷領域は９３ｋＰａを超える領域とされる。
【００６７】
ステップＳ５３では、クランク角インデクスｉ及び周波数インデクスｊをともに「０」に初期化するとともに、強度積算値ＳＵＭＫ及び基準積算値ＳＵＭＭを「０」に初期化する。強度積算値ＳＵＭＫ及び基準積算値ＳＵＭＭは、後述するステップＳ５７で更新され、ステップＳ６０で適合率ＰＦＩＴの算出に適用される。
【００６８】
ステップＳ５４では、周波数インデクスｊが周波数データ数ＪＮから「１」を減算した値より大きいか否かを判別する。最初はこの答は否定（ＮＯ）であるので、ステップＳ５５に進み、クランク角インデクスｉがクランク角データ数ＩＮから「１」を減算した値より大きいか否かを判別する。
【００６９】
最初はステップＳ５５の答も否定（ＮＯ）であるので、ステップＳ５６に進み、下記式（７）及び（８）により、重み付けマスタパラメータＭＭＷ及び重み付け積パラメータＫＭＷを算出する。下記式のＷＭＡＰ(j,i)は、重み付けマップに設定されている重み付けパラメータである。重み付け積パラメータＫＭＷは、マスタパラメータＭＭＡＰ(j,i)と補正二値化強度パラメータＪＫＭＡＰ(j,i)との積を、重み付けパラメータＷＭＡＰ(j,i)によって重み付けしたものである。
ＭＭＷ＝ＭＭＡＰ(j,i)×ＷＭＡＰ(j,i)                        （７）
ＫＭＷ＝ＭＭＡＰ(j,i)×ＪＫＭＡＰ(j,i)×ＷＭＡＰ(j,i)       （８）
【００７０】
ステップＳ５７では、下記式（９）及び（１０）により、重み付けマスタパラメータＭＭＷ及び重み付け積パラメータＫＭＷを積算し、基準積算値ＳＵＭＭ及び強度積算値ＳＵＭＫを算出する。
ＳＵＭＭ＝ＳＵＭＭ＋ＭＭＷ                                （９）
ＳＵＭＫ＝ＳＵＭＫ＋ＫＭＷ                                （１０）
【００７１】
ステップＳ５８では、クランク角インデクスｉを「１」だけインクリメントし、ステップＳ５５に戻る。ステップＳ５５の答が否定（ＮＯ）である間はステップＳ５６〜Ｓ５８を繰り返し実行し、クランク角インデクスｉが（ＩＮ−１）を超えると、ステップＳ５９に進み、クランク角インデクスｉを「０」に戻すとともに、周波数インデクスｊを「１」だけインクリメントし、ステップＳ５４に戻る。ステップＳ５４の答が否定（ＮＯ）である間は、ステップＳ５５〜Ｓ５９を繰り返し実行し、周波数インデクスｊが（ＪＮ−１）を超えると、ステップＳ６０に進み、下記式（１１）により適合率ＰＦＩＴを算出する。
ＰＦＩＴ＝ＳＵＭＫ／ＳＵＭＭ                              （１１）
【００７２】
図２０は、図１２のステップＳ１８で実行されるノイズ学習処理のフローチャートである。
【００７３】
ステップＳ７１では、クランク角インデクスｉ、周波数インデクスｊ、加算学習パラメータＬＫ、及び減算学習パラメータＬＭをいずれも「０」に初期化する。ステップＳ７２では、周波数インデクスｊが周波数データ数ＪＮから「１」を減算した値より大きいか否かを判別する。最初はこの答は否定（ＮＯ）であるので、ステップＳ７３に進み、クランク角インデクスｉがクランク角データ数ＩＮから「１」を減算した値より大きいか否かを判別する。
【００７４】
最初はステップＳ７３の答も否定（ＮＯ）であるので、ステップＳ７４に進み、二値化強度パラメータＮＫＭＡＰ(j,i)が二値化ノイズパラメータＮＮＭＡＰ(j,i)と等しいか否かを判別する。この答が肯定（ＹＥＳ）であるときは、直ちにステップＳ８０に進む。
【００７５】
ステップＳ７４の答が否定（ＮＯ）であって、二値化強度パラメータＮＫＭＡＰ(j,i)が二値化ノイズパラメータＮＮＭＡＰ(j,i)と異なるときは、二値化強度パラメータＮＫＭＡＰ(j,i)が二値化ノイズパラメータＮＮＭＡＰ(j,i)より大きいか否かを判別する（ステップＳ７５）。この答が肯定（ＹＥＳ）であるときは、加算学習パラメータＬＫを「１」に設定し、減算学習パラメータＬＭを「０」に設定する（ステップＳ７６）。一方、ＮＫＭＡＰ(j,i)＜ＮＮＭＡＰ(j,i)であるときは、加算学習パラメータＬＫを「０」に設定し、減算学習パラメータＬＭを「１」に設定する（ステップＳ７７）。
【００７６】
ステップＳ７８では、下記式（１２）及び（１３）により、加算学習パラメータＬＫ及び減算学習パラメータＬＭを修正する。式（１２）、（１３）のＤＳＮＯＩＳＥは、例えば０．１に設定されるノイズ学習係数である。
ＬＫ＝ＤＳＮＯＩＳＥ×ＬＫ                                （１２）
ＬＭ＝ＤＳＮＯＩＳＥ×ＬＭ                                （１３）
【００７７】
ステップＳ７９では、下記式（１４）に加算学習パラメータＬＫ及び減算学習パラメータＬＭを適用し、ノイズパラメータＮＭＡＰ(j,i)を更新する。次に説明する図２１の処理で、ノイズパラメータＮＭＡＰ(j,i)を二値化することにより、二値化ノイズパラメータＮＮＭＡＰ(j,i)が算出される。
ＮＭＡＰ(j,i)＝ＮＭＡＰ(j,i)＋ＬＫ−ＬＭ                （１４）
【００７８】
ステップＳ８０ではクランク角インデクスｉを「１」だけインクリメントし、ステップＳ７３に戻る。ステップＳ７３の答が否定（ＮＯ）である間はステップＳ７４〜Ｓ８０を繰り返し実行し、クランク角インデクスｉが（ＩＮ−１）を超えると、ステップＳ８１に進み、クランク角インデクスｉを「０」に戻すとともに、周波数インデクスｊを「１」だけインクリメントし、ステップＳ７２に戻る。ステップＳ７２の答が否定（ＮＯ）である間は、ステップＳ７３〜Ｓ８１を繰り返し実行し、周波数インデクスｊが（ＪＮ−１）を超えると、ステップＳ８２に進み、図２１に示すノイズマップ更新処理を実行する。
【００７９】
図２１のステップＳ９１では、クランク角インデクスｉ及び周波数インデクスｊをともに「０」に初期化する。ステップＳ９２では、周波数インデクスｊが周波数データ数ＪＮから「１」を減算した値より大きいか否かを判別する。最初はこの答は否定（ＮＯ）であるので、ステップＳ９３に進み、クランク角インデクスｉがクランク角データ数ＩＮから「１」を減算した値より大きいか否かを判別する。
【００８０】
最初はステップＳ９３の答も否定（ＮＯ）であるので、ステップＳ９４に進み、ノイズパラメータＮＭＡＰ(j,i)がノイズ二値化閾値ＮＬＶＬ（例えば０．８）より大きいか否かを判別する。この答が肯定（ＹＥＳ）であるときは、二値化ノイズパラメータＮＮＭＡＰ(j,i)を「１」に設定する（ステップＳ９５）。一方、ＮＭＡＰ(j,i)≦ＮＬＶＬであるときは、二値化ノイズパラメータＮＮＭＡＰ(j,i)を「０」に設定する（ステップＳ９６）。
【００８１】
ステップＳ９７ではクランク角インデクスｉを「１」だけインクリメントし、ステップＳ９３に戻る。ステップＳ９３の答が否定（ＮＯ）である間はステップＳ９４〜Ｓ９７を繰り返し実行し、クランク角インデクスｉが（ＩＮ−１）を超えると、ステップＳ９８に進み、クランク角インデクスｉを「０」に戻すとともに、周波数インデクスｊを「１」だけインクリメントし、ステップＳ９２に戻る。ステップＳ９２の答が否定（ＮＯ）である間は、ステップＳ９３〜Ｓ９８を繰り返し実行し、周波数インデクスｊが（ＪＮ−１）を超えると、本処理を終了する。
【００８２】
図２０のステップＳ７５でＮＫＭＡＰ(j,i)＞ＮＮＭＡＰ(j,i)であるときは、ステップＳ７８でＬＫ＝０．１，ＬＭ＝０となり、式（１４）により、ノイズパラメータＮＭＡＰ(j,i)が「０．１」だけインクリメントされる。一方、ステップＳ７５でＮＫＭＡＰ(j,i)＜ＮＮＭＡＰ(j,i)であるときは、ステップＳ７８でＬＫ＝０，ＬＭ＝０．１となり、式（１４）により、ノイズパラメータＮＭＡＰ(j,i)が「０．１」だけデクリメントされる。そして、図２１の処理でノイズパラメータＮＭＡＰ(j,i)がノイズ二値化閾値ＮＬＶＬより大きいとき、二値化ノイズパラメータＮＮＭＡＰ(j,i)が「１」に設定される一方、ノイズパラメータＮＭＡＰ(j,i)がノイズ二値化閾値ＮＬＶＬ以下であるときは、二値化ノイズパラメータＮＮＭＡＰ(j,i)が「０」に設定される。
【００８３】
図２０及び図２１の処理により、ノッキングが発生していないと判定されたときの二値化強度パラメータＮＫＭＡＰ(j,i)に応じてノイズマップが更新され、例えば吸気弁の着座ノイズのように定常的に発生するノイズがノイズマップに反映される。その結果、ノイズの影響を除いて高精度の判定を行うことが可能となる。
【００８４】
なお、ノイズ除去処理（図１２，ステップＳ１２）を行わずに、二値化強度パラメータＮＫＭＡＰをそのまま用いて適合率ＰＦＩＴの算出を行うようにしてよい。ノイズ除去処理を行わない場合には、ノイズの影響を受ける可能性が高くなるが、図１９を参照して説明したようにノイズと区別してノッキングの判定を行うことができるからである。
【００８５】
以上詳述したように本実施形態では、クランク角６度間隔でノックセンサ１１の出力信号の周波数成分分析が行われ、その結果得られる、５ｋＨｚ〜２５ｋＨｚの周波数成分強度の時系列データであるスペクトル時系列マップが生成される。すなわち、スペクトル時系列マップの要素が二次元配列データである強度パラメータＫＭＡＰ(j,i)としてメモリに格納される。そして、強度パラメータＫＭＡＰ(j,i)を二値化することにより、二値化強度パラメータＮＫＭＡＰ(j,i)が算出され、該二値化強度パラメータＮＫＭＡＰ(j,i)に基づいてノッキングが発生したか否かが判定される。二値化強度パラメータＮＫＭＡＰ(j,i)には、エンジンの回転に伴う周波数成分分布の変化が反映されるので、二値化強度パラメータＮＫＭＡＰ(j,i)と、ノッキング発生時に特有の変化パターンに対応するマスタパターンマップ上のマスタパラメータＭＭＡＰ(j,i)とを比較することにより、ノッキングの発生を正確に判定することができる。また周波数成分分析の結果得られる強度パラメータＫＭＡＰ(j,i)を二値化することにより、データ量が減少するとともに時系列データの変化パターンが単純化されるので、メモリ容量を低減するとともに演算速度を高めることができる。
【００８６】
また二値化強度パラメータＮＫＭＡＰ(j,i)（ＪＫＭＡＰ(j,i)）が、マスタパラメータＭＭＡＰ(j,i)に近い変化パターンを示すときに、ノッキングが発生している可能性が高いので、二値化強度パラメータＮＫＭＡＰ(j,i)（ＪＫＭＡＰ(j,i)）とマスタパラメータＭＭＡＰ(j,i)との類似性（相関性）を示すパラメータを算出することにより、正確な判定を行うことができる。
【００８７】
本実施形態ではこの類似性（相関性）を示すパラメータとして適合率ＰＦＩＴを用い、適合率ＰＦＩＴが判定閾値ＳＬＶＬを超えたときにノッキングが発生したと判定される。適合率ＰＦＩＴを用いることにより、二値化強度パラメータＮＫＭＡＰ(j,i)（ＪＫＭＡＰ(j,i)）とマスタパラメータＭＭＡＰ(j,i)との類似性（相関性）を比較的簡単な演算で的確に評価し、正確な判定を行うことができる。
【００８８】
また二値化強度パラメータＮＫＭＡＰ(j,i)に基づいて、ノイズ成分の時系列データＮＮＭＡＰ(j,i)が算出され、二値化強度パラメータＮＫＭＡＰ(j,i)がノイズ成分の時系列データであるノイズパラメータＮＮＭＡＰ(j,i)により補正され、補正二値化強度パラメータＪＫＭＡＰ(j,i)に基づいてノッキング判定が行われる。したがって、上述した着座ノイズのように定常的に表れるノイズ成分を除去してより正確な判定を行うことが可能となる。
【００８９】
また二値化強度パラメータＮＫＭＡＰ(j,i)（ＪＮＫＭＡＰ(j,i)）及びマスタパラメータＭＭＡＰ(j,i)に対して、周波数に応じて設定された重み付けパラメータＷＭＡＰ(j,i)を乗算して、適合率ＰＦＩＴが算出される。ノッキングが発生したときに大きくなる周波数成分は予め判明しているので、その周波数の近傍の周波数に対応するパラメータに大きな重みを付けることにより、判定精度を高めることができる。
【００９０】
本実施形態では、ＥＣＵ５がサンプリング手段、要素強度算出手段、及び周波数成分強度算出手段を構成する。
【００９１】
なお本発明は上述した実施形態に限るものではなく、種々の変形が可能である。例えば、ノックセンサ出力のサンプリング周期や周波数成分分析を行うクランク角度間隔は上述したもの（２０マイクロ秒、６度）に限るものではなく、本発明の目的が達成される範囲内において変更可能である。また、二値化スペクトル時系列マップ（上述した実施形態では２１行×１５列のマトリクスで構成）も同様に変更可能である。
【００９２】
また上述した実施形態では、ノックセンサによる検出値の周波数分析を行う装置を示したが、本願発明はこれに限らず、エンジン運転パラメータ検出値の周波数成分分析を、エンジン回転（ＣＲＫ割り込み）に同期して行う場合に適用可能である。例えば本発明は、ＣＲＫセンサから出力されるＣＲＫパルス（割り込み）の発生周期ＣＲＭＥまたはその逆数として得られるエンジン１の瞬時回転速度の周波数成分分析などにも適用できる。
【００９３】
また上述した実施形態では、エンジン回転数ＮＥの閾値ＮＥＴＨを、単位時間当たりの演算ステップ数が飽和するステップ数飽和回転数ＮＥＳＡＴに設定したが、周波数成分強度の算出に使用するデータのサンプリング期間ＴＳが、ＣＲＫ割り込み発生周期ＣＲＭＥより短くなる回転数ＮＥＴＳ、または回転数ＮＥＴＳより若干高い（例えば１００ｒｐｍ程度高い）回転数を閾値ＮＥＴＨに設定するようにしてもよい。上述した実施形態では、回転数ＮＥＴＳは、１０００ｒｐｍである。
【００９４】
また上述した実施形態では、二値化強度パラメータＮＫＭＡＰ(j,i)及びマスタパラメータＭＭＡＰ(j,i)に重み付けパラメータＷＭＡＰ(j,i)を乗算して、適合率ＰＦＩＴを算出したが、重み付けパラメータＷＭＡＰ(j,i)を乗算せずに、すなわち重み付けを行わずに算出するようにしてもよい。
【００９５】
また上述した実施形態では、判定閾値ＳＬＶＬ、二値化閾値ＢＬＶＬ、マスタパターンマップ、及び重み付けマップをエンジン回転数ＮＥ及び吸気圧ＰＢＡに応じて、算出または選択するようにしたが、予め設定された値または１つのマップに固定しておいてもよい。
【符号の説明】
【００９６】
１  内燃機関
５  電子制御ユニット（サンプリング手段、要素強度算出手段、周波数成分強度算出手段）
１１  ノックセンサ
【図１】
2010035634.tif 000004
【図２】
2010035634.tif 000005
【図３】
2010035634.tif 000006
【図４】
2010035634.tif 000007
【図５】
2010035634.tif 000008
【図６】
2010035634.tif 000009
【図７】
2010035634.tif 000010
【図８】
2010035634.tif 000011
【図９】
2010035634.tif 000012
【図１０】
2010035634.tif 000013
【図１１】
2010035634.tif 000014
【図１２】
2010035634.tif 000015
【図１３】
2010035634.tif 000016
【図１４】
2010035634.tif 000017
【図１５】
2010035634.tif 000018
【図１６】
2010035634.tif 000019
【図１７】
2010035634.tif 000020
【図１８】
2010035634.tif 000021
【図１９】
2010035634.tif 000022
【図２０】
2010035634.tif 000023
【図２１】
2010035634.tif 000024
【図２２】
2010035634.tif 000025
【国際調査報告】
2010035634.tif 000026
2010035634.tif 000027
2010035634.tif 000028
2010035634.tif 000029