JP4764485B2 - 周波数成分分析装置 - Google Patents

Description

本発明は、内燃機関の運転パラメータの検出値についての周波数成分分析を、機関回転に同期して行う周波数成分分析装置に関する。

特許文献１には、内燃機関に装着されたノックセンサの周波数成分強度を、アダマール変換を用いて近似的に算出し、算出された周波数成分強度に基づいてノッキング判定を行うノッキング検出装置が示されている。

また特許文献２には、離散フーリエ変換（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)は、通常は下記式（Ａ）により算出される。ここでインデクスパラメータｍは、ＣＲＫ割り込みが発生したときのアドレス番号ｋが「２５」となるように設定される修正アドレス番号である（図３（ｄ）参照）。

式（Ａ）の算出には、サイン波成分強度ＤＭＦＴＳ及びコサイン波成分強度ＤＭＦＴＣの積算値ＳＵＭＤＳ及びＳＵＭＤＣを算出する必要がある。例えば周波数５ｋＨｚの積算値ＳＵＭＤＳ（5kHz)は、図７に示すテーブルを用いて下記式（３）により算出される。式（３）では、量子化サイン波信号ＱＳＩＮの値が「０」となる項は省かれている。また積算値ＳＵＭＤＣ（5kHz)は、図８に示すテーブルを用いて下記式（４）により算出される。式（４）では、量子化コサイン波信号ＱＣＯＳの値が「０」となる項は省かれている。
ＳＵＭＤＳ（5kHz)＝ＶＫＮＫ(2)＋ＶＫＮＫ(3)＋ＶＫＮＫ(4)＋ＶＫＮＫ(5)
−ＶＫＮＫ(7)−ＶＫＮＫ(8)＋……
−ＶＫＮＫ(49)−ＶＫＮＫ(50) （３）
ＳＵＭＤＣ（5kHz)＝ＶＫＮＫ(1)＋ＶＫＮＫ(2)−ＶＫＮＫ(5)−ＶＫＮＫ(6)
−ＶＫＮＫ(7)＋ＶＫＮＫ(10)＋……
−ＶＫＮＫ(47)＋ＶＫＮＫ(50) （４）

他の周波数に対応する積算値ＳＵＭＤＳ，ＳＵＭＤＣも同様に算出することができる。このように量子化サイン波信号ＱＳＩＮ及び量子化コサイン波信号ＱＣＯＳを用いることにより、積算値ＳＵＭＤＳ，ＳＵＭＤＣの算出に乗算演算が不要となり、周波数成分強度ＳＴＦＴの算出に必要な演算量を大幅に低減することができる。その結果、ＣＰＵにおける演算ステップ数を量子化を行わない場合に比べて６６％程度低減することができ、算出された周波数成分強度ＳＴＦＴを用いた後処理の時間を十分に確保することが可能となる。また３値化正弦波信号を用いることにより、アダマール変換を適用する場合に比べて周波数成分強度ＳＴＦＴの精度を高めることができる。

次に３値化に適用する閾値ＴＨ３ＰＳ，ＴＨ３ＰＣの設定について説明する。なお、上述したようにマイナス側の閾値ＴＨ３ＭＳ，ＴＨ３ＭＣは、プラス側の閾値ＴＨ３ＰＳ，ＴＨ３ＰＣに負号を付したものである。図９は、閾値ＴＨ３ＰＳと、量子化誤差の大きさを示す誤差率ＲＥＲＲＳとの関係を示す図であり、同図（ａ）が周波数６ｋＨｚに対応し、同図（ｂ）が周波数７ｋＨｚに対応する。誤差率ＲＥＲＲＳが最小となる閾値ＴＨ３ＰＳの値は、６ｋＨｚでは「０．５８９」であり、７ｋＨｚでは「０．４８１」となり、サイン波信号の周波数に依存して変化する。

これは、一定のサンプリング周期ＴＳＭＰ（２０μｓ）でサンプリングすると、図１０に示すようにサイン波信号の周波数に依存してサンプリング位相が変化するためである。本実施形態では、閾値ＴＨ３ＰＳ，ＴＨ３ＰＣは、誤差率ＲＥＲＲＳ，ＲＥＲＲＣが最小となる値に設定され、信号周波数に応じて図１１に示すように設定されている。このように閾値ＴＨ３ＰＳ，ＴＨ３ＰＣ（ＴＨ３ＭＳ，ＴＨ３ＭＣ）を設定することにより、量子化誤差の影響を最小限に抑制することができる。

なお、上記誤差率ＲＥＲＲＳ，ＲＥＲＲＣは、下記式（５）〜（１０）を用いて算出される。下記式において「ＮＦ」は、周波数成分強度を算出する周波数の数（本実施形態では「２１」）であり、「ｆ」は、誤差率を算出する対象周波数（ｋＨｚ）である。

図１２は周波数成分強度ＳＴＦＴ算出の手順を示すフローチャートである。ステップＳ１では、検出データＶＫＮＫを取得し、ステップＳ２では、前記式（１）及び（２）によりサイン波成分強度ＤＭＦＴＳ(j,k)及びコサイン波成分強度ＤＭＦＴＣ(j,k)を算出する。その際上述したように図７及び図８に示すＱＳＩＮテーブル及びＱＣＯＳテーブルを使用する。

ステップＳ１及びＳ２の処理は、サンプリング周期（２０マイクロ秒）毎に実行される。

ステップＳ３では、ＣＲＫ割り込みが発生したか否かを判別し、発生していないときはステップＳ１に戻り、ＣＲＫ割り込みが発生したタイミングで、エンジン回転数ＮＥが回転数閾値ＮＥＴＨ以上であるか否かを判別する。ＮＥ＜ＮＥＴＨであるときは、前回使用データを使用しない前記式（Ａ）により周波数成分強度ＳＴＦＴを算出する一方、ＮＥ≧ＮＥＴＨであるときは、前回使用データを使用する下記式（Ｂ）により周波数成分強度ＳＴＦＴ(j,i)を算出する。インデクスパラメータｍｐは、前回算出時に適用された修正アドレス番号であり、ｍｘは前回境界インデクスであり、ｍｙは今回境界インデクスであり、それぞれ下記式（１１）及び（１２）で与えられる。式（Ｂ）の演算においても図７及び図８に示すテーブルが使用され、演算量が低減される。

ｍｘ＝ＮＤ−ＮＶ＋１ （１１）
ｍｙ＝ＮＶ＋１ （１２）
ここでＮＶは前回のＳＴＦＴ算出に使用したデータのうち、今回のＳＴＦＴ算出に使用できるデータの数であり、データ数ＮＶはエンジン回転数ＮＥが高くなるほど増加する。

式（Ｂ）の第１項ＴＲＭＳ１及び第３項ＴＲＭＣ１が、積算値の前回算出値に相当し、第２項ＴＲＭＳ２及び第４項ＴＲＭＣ２が新たに算出する積算値に相当する。

このようにサイン波成分強度の積算値の一部及びコサイン波成分強度の積算値一部を、それぞれ前回算出値（ＴＲＭＳ１，ＴＲＭＣ１）に置き換えて、周波数成分強度ＳＴＦＴ(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)を乗算せずに、すなわち重み付けを行わずに算出するようにしてもよい。

また上述した実施形態では、判定閾値ＳＬＶＬ、二値化閾値ＢＬＶＬ、マスタパターンマップ、及び重み付けマップをエンジン回転数ＮＥ及び吸気圧ＰＢＡに応じて、算出または選択するようにしたが、予め設定された値または１つのマップに固定しておいてもよい。

１ 内燃機関
５ 電子制御ユニット（サンプリング手段、要素強度算出手段、量子化手段、周波数成分強度算出手段）
１１ ノックセンサ

Claims (3)

1. 内燃機関の運転パラメータの検出値についての周波数成分分析を、前記機関の回転に同期して行う周波数成分分析装置において、
   前記運転パラメータを所定時間間隔でサンプリングし、サンプル値をデジタル値に変換するサンプリング手段と、
   前記検出値に含まれる複数の周波数成分に対応する、第１要素の強度及び前記第１要素に対して位相が９０度ずれた第２要素の強度を、所定数のサンプル値について算出する要素強度算出手段と、
   前記第１要素強度及び前記第２要素強度を用いて前記複数周波数の周波数成分強度を、前記機関の回転に同期して算出する周波数成分強度算出手段とを備え、
   前記要素強度算出手段は、前記第１及び第２要素に対応する正弦波信号をそれぞれＮ値化（Ｎは３または５）し、Ｎ値化正弦波信号を算出する量子化手段を有し、
   前記要素強度算出手段は、前記運転パラメータの検出値と、前記Ｎ値化正弦波信号とを用いて前記第１及び第２要素の強度を算出することを特徴とする周波数成分分析装置。
2. 前記量子化手段は、前記Ｎ値化する際に適用する閾値を前記正弦波信号の周波数に応じて変更することを特徴とする請求項１に記載の周波数成分分析装置。
3. 前記量子化手段は、前記機関の回転数が所定回転数以下であるときに５値化を行い、前記機関の回転数が前記所定回転数より高いときに３値化を行うことを特徴とする請求項１または２に記載の周波数成分分析装置。

Images