"音響信号の周波数信号を分析する周波数分析手段と、前記周波数信号の位相の時間変化を近似する位相曲線を算出する位相曲線算出手段と、前記位相曲線と前記周波数信号の位相との誤差を算出する誤差算出手段と、前記誤差に基づいて、前記音響信号が周期音の信号か否かを識別する音響信号識別手段とを備える音識別装置。前記音響信号識別手段は、前記誤差の算出に用いた前記周波数信号の位相に対応する音響信号を、前記誤差が所定の閾値以上の場合に非周期音の信号として識別し、前記誤差が前記所定の閾値未満の場合に周期音の信号として識別する請求項1記載の音識別装置。前記音響信号識別手段は、所定の時間幅に含まれる前記誤差の算出に用いた前記周波数信号の位相に対応する音響信号を、当該所定の時間幅に含まれる前記誤差の和又は平均値が所定の閾値以上の場合に非周期音の信号として識別し、前記誤差の和又は平均値が前記所定の閾値未満の場合に周期音の信号として識別する請求項1記載の音識別装置。さらに、所定数の前記位相との差が小さくなるように、前記所定数の前記位相とは異なる他の前記位相に±2π×mラジアンmは自然数を加算することにより、他の前記位相を補正する位相補正手段を備える請求項1記載の音識別装置。さらに、互いに異なる角度範囲ごとに、当該角度範囲内に収まるように前記位相に±2π×mラジアンmは自然数を加算することにより、前記位相を補正する位相補正手段を備え、前記位相曲線算出手段は、前記角度範囲ごとに、前記位相曲線を算出し、前記誤差算出手段は、前記角度範囲ごとに、前記誤差を算出し、前記位相補正手段は、さらに、前記誤差が最小となるときの角度範囲を選択し、前記音響信号識別手段は、選択された前記角度範囲における前記誤差に基づいて、前記音響信号が周期音の信号か否かを識別する請求項1記載の音識別装置。前記音響信号は、混合音の音響信号であり、前記音響信号識別手段は、前記誤差が所定の閾値未満の場合に、前記誤差の算出に用いた前記周波数信号の位相に対応する音響信号を、エンジン音の信号として識別する請求項1記載の音識別装置。前記周波数分析手段は、各々が音響信号の入力を受け付ける互いに離間して配置された複数のマイクロホンで受け付けられた複数の音響信号のそれぞれについて周波数信号を分析し、前記音識別装置は、さらに、前記音響信号識別手段により少なくとも1つのマイクロホンから入力された音響信号が周期音の信号として識別された場合に、前記複数のマイクロホンで受け付けられた複数の前記音響信号の到達時間差に基づいて、前記周期音の音源方向を検知する方向検知部を備える請求項1記載の音識別装置。前記位相曲線は、位相の値を変数とする2次多項式で表される曲線である請求項1記載の音識別装置。前記誤差算出手段は、前記周波数信号の位相と、当該周波数信号の時刻と同時刻における前記位相曲線の値との差を、前記誤差として算出する請求項1記載の音識別装置。音響信号の周波数信号を分析するステップと、前記周波数信号の位相の時間変化を近似する位相曲線を算出するステップと、前記位相曲線と前記周波数信号の位相との誤差を算出するステップと、前記誤差に基づいて、前記音響信号が周期音の信号か否かを識別するステップとを含む音識別方法。音響信号の周波数信号を分析するステップと、前記周波数信号の位相の時間変化を近似する位相曲線を算出するステップと、前記位相曲線と前記周波数信号の位相との誤差を算出するステップと、前記誤差に基づいて、前記音響信号が周期音の信号か否かを識別するステップとをコンピュータに実行させるためのプログラム。第1の従来技術では、スペクトルサブトラクション法SS法によって雑音を除去している。SS法とは、一定時間の音を周波数分析し、得られた周波数ごとのパワーを雑音として減算することで特定の時間の音を抽出する方法である。そのため、事前に雑音を推定する必要がある。定常的なパワーを有する音が周囲の騒音下にある場合は、雑音を推定することができるため、雑音を除去することが可能となる。しかしながら、風雑音など非定常雑音はパワーが時間的に変化するため、このような非定常雑音に対しては、SS法は頑強性が低く、風雑音と車両音とを精度良く判別して車両音を検知することはできない。第2の従来技術は、エンジン音などの音は周期音周波数はほぼ一定であり、位相が時刻に対して一定の値になる性質を利用し、周期音を識別するものである。車両が定常走行同じ速度での走行しており、エンジンの回転数が一定エンジン音の周波数が時刻に対して一定の値になるの場合は、周期音を識別することができる。一方、車両の加速又は減速によりエンジンの回転数が変化する場合は、識別精度を向上させるために周波数の時間変化に対してさらに対策を行う必要がある。特に、自車両の死角位置に存在する車両死角車両を検知するアプリケーション等では、より安全に運転を支援するためにも、衝突すると大きな事故につながりやすい加速している車両を正確に検知することは重要な課題となる。本発明は、従来の課題を解決するもので、エンジン音や音声などの周期音と、風雑音、雨音、暗騒音などの非周期音とを区別して、周期音もしくは非周期音の周波数信号を判定する音識別装置に関し、特にエンジン音や音声などのように周波数が時間とともに変化する周期音を正確に判定する音識別装置を提供することを目的とする。上記目的を達成するために、本発明のある局面に係る音識別装置は、音響信号の周波数信号を分析する周波数分析手段と、前記周波数信号の位相の時間変化を近似する位相曲線を算出する位相曲線算出手段と、前記位相曲線と前記周波数信号の位相との誤差を算出する誤差算出手段と、前記誤差に基づいて、前記音響信号が周期音の信号か否かを識別する音響信号識別手段とを備える。周波数が時間とともに変化する場合には位相も時間とともに変化する。このような位相の時間変化を位相曲線として表し、位相曲線との誤差に基づいて周期音の信号か否かを識別することにより、エンジン音や音声などの周期音と、風雑音、雨音、暗騒音などの非周期音とを区別して、周期音もしくは非周期音の周波数信号を判定することができる。特にエンジン音や音声などのように周波数が時間とともに変化する周期音を正確に判定することができる。好ましくは、上記音識別装置は、さらに、所定数の前記位相との差が小さくなるように、前記所定数の前記位相とは異なる他の前記位相に±2π×mラジアンmは自然数を加算することにより、他の前記位相を補正する位相補正手段を備える。これにより、他の時刻の位相と大きく外れた位相を補正することができ、周期音の識別を精度良く行うことが可能となる。また、上記音識別装置は、さらに、互いに異なる角度範囲ごとに、当該角度範囲内に収まるように前記位相に±2π×mラジアンmは自然数を加算することにより、前記位相を補正する位相補正手段を備え、前記位相曲線算出手段は、前記角度範囲ごとに、前記位相曲線を算出し、前記誤差算出手段は、前記角度範囲ごとに、前記誤差を算出し、前記位相補正手段は、さらに、前記誤差が最小となるときの角度範囲を選択し、前記音響信号識別手段は、選択された前記角度範囲における前記誤差に基づいて、前記音響信号が周期音の信号か否かを識別してもよい。これにより、他の時刻の位相と大きく外れた位相を補正することができ、周期音の識別を精度良く行うことが可能となる。さらに好ましくは、前記周波数分析手段は、各々が音響信号の入力を受け付ける互いに離間して配置された複数のマイクロホンで受け付けられた複数の音響信号のそれぞれについて周波数信号を分析し、前記音識別装置は、さらに、前記音響信号識別手段により少なくとも1つのマイクロホンから入力された音響信号が周期音の信号として識別された場合に、前記複数のマイクロホンで受け付けられた複数の前記音響信号の到達時間差に基づいて、前記周期音の音源方向を検知する方向検知部を備える。周期音が識別された場合に、音響信号の到達時間差をもとに車両の方向を検知している。このため、雑音の影響を受けずに車両の方向を正確に検知することができる。さらに好ましくは、前記位相曲線は、位相の値を変数とする2次多項式で表される曲線である。音響信号の周波数の変化が1次方程式で表すことができる場合には、位相の変化は2次方程式で表すことができる。このため、位相曲線を2次多項式で表される曲線とすることにより、位相変化を精度良く表すことができる。なお、本発明は、このような特徴的な手段を備える音識別装置として実現することができるだけでなく、音識別装置に含まれる特徴的な手段をステップとする音識別方法として実現したり、音識別方法に含まれる特徴的なステップをコンピュータに実行させるプログラムとして実現したりすることもできる。そして、そのようなプログラムは、CD-ROMCompactDiscReadOnlyMemory等の記録媒体やインターネット等の通信ネットワークを介して流通させることができるのは言うまでもない。本発明の音識別装置等によれば、エンジン音や音声などの周期音と、風雑音、雨音、暗騒音などの非周期音とを区別して、周期音もしくは非周期音の周波数信号を判定できる。特にエンジン音などのように周波数が時間とともに変化する周期音を正確に判定する音識別装置等を提供することができる。本発明の特徴は、エンジン音や音声などの周期音の周波数の時間変化の特徴に着目したことである。本願発明者らは、音の発生メカニズムや実際に集音した音声データを分析した。その結果、エンジン音や音声などの周期音の時間-周波数空間での周波数の時間変化が線形区分で近似できるという新たな知見を発見した。この知見により線形区分近似された音の位相の時間変化は曲線でモデル化できることを見出した。これにより、周波数が時間変化する場合においても正確に周期音を識別することができる。なお、本願発明における周期音とは、位相が一定の音又は位相の変化が連続的な音を示す。ここで、図1を用いて本発明で用いる位相の定義を行う。図1aには、入力したエンジン音の例が模式的に示されている。横軸は時間を表しており、縦軸は振幅を表している。ここではエンジンの回転数が時刻に対して一定であり、エンジン音の周波数が変化しない場合の例が示されている。また、図1bには、フーリエ変換を用いて周波数分析を行う場合の基底波形である周波数fの正弦波ここではエンジン音の周波数と同じ値を所定の周波数fとしているが示されている。横軸と縦軸は図1aと同じである。この基底波形と入力した混合音との畳み込み処理を行うことで周波数信号位相を求める。この例では、基底波形を時間軸方向に移動させずに固定し、入力したエンジン音と畳み込み処理を行うことで、時刻ごとの周波数信号位相を求めている。この処理で求めた結果を図1cに示す。横軸は時間を表しており、縦軸は位相を表している。この例では、エンジンの回転数が時刻に対して一定であり、入力したエンジン音の周波数が時刻に対して一定である。このため、所定の周波数fでの位相は加速度的に増加又は減少はしていない。この例では回転数が一定であるエンジン音の周波数と同じ値を所定の周波数fとしたが、エンジン音の周波数よりも小さい値を所定の周波数fとした場合には位相は一次関数的に増加する。また、エンジン音の周波数よりも大きい値を所定の周波数fとした場合には、位相は一次関数的に減少する。いずれの場合も、所定の周波数fでの位相は、加速度的に増加又は減少はしていない。なお、音声信号分野や高速フーリエ変換FFTなどでは基底波形を時間軸方向にずらしながら畳み込みを行うのが一般的である。この基底波形を時間軸方向にずらしながら畳み込みを行う場合は後に位相を補正することで本発明で定義する位相へと変換することが可能である。以下、図を用いて説明する。図2は位相を説明する図である。図2aには、入力したエンジン音の例が模式的に示されている。横軸は時間を表しており、縦軸は振幅を表している。また、図2bには、フーリエ変換を用いて周波数分析を行う場合の基底波形である周波数fの正弦波ここではエンジン音の周波数と同じ値を所定の周波数fとしているが示されている。横軸と縦軸は図2aと同じである。この基底波形と入力した混合音との畳み込み処理を行うことで周波数信号位相を求める。この例では、基底波形を時間軸方向に移動させながら、入力したエンジン音と畳み込み処理を行うことで、時刻ごとの周波数信号位相を求めている。この処理で求めた結果を図2cに示す。横軸は時間を表しており、縦軸は位相を表している。入力したエンジン音は周波数fであるため、周波数fでの位相のパターンは、1/fの時刻の周期で規則的に繰り返されることとなる。そこで算出された位相ψtから規則的に繰り返される位相を補正ψ ́t=mod2πψt-2πftfは分析周波数すること図2dに示すような位相が得られる。つまり、位相補正を行うことにより、図1cに示す、本発明で定義される位相へと変換することが可能となる。次に、エンジン回転数に伴うエンジン音の時刻に対する周波数の変化について説明する。図3は、本実施の形態の特徴を説明する概念図である。図3は自動車のエンジン音を後述するDFT分析部2402において分析したスペクトログラムである。縦軸は周波数を、横軸は時間をそれぞれ示しており、色の濃度は周波数信号のパワーの大きさを示している。濃い色黒い色はパワーが大きいことを示している。図3は、風などの雑音をなるべく除去したデータであり、色の濃い部分黒っぽい部分がおおむねエンジン音を示している。一般的にこのようにエンジン音は、時間とともに回転数が変化したデータであり、スペクトログラムから時間の経過と共に周波数が変化していることが分かる。エンジンは、所定数のシリンダーがピストン運動を行うことで駆動系を回転させている。そして、車両から発せられるエンジン音は、このエンジンの回転に依存した音と、エンジンの回転には依存しない固定振動音や非周期音とからなる。特に車両の外部から検知できる主な音は、エンジンの回転に依存した周期音であり、本実施の形態では、このエンジンの回転に依存する周期音をエンジン音として抽出する。図3の点線の円501、502及び503に示すように、エンジン音は、回転数が変化することで周波数が部分的に時刻に応じて変化している。ここで、周波数の変化に着目すると、周波数がランダムに変化したり、離散的に飛んだりすることはほとんどなく、所定の時間間隔でみると線形モデルで近似できることが分かる。したがって、エンジン音は下記の式1で示すような、の区分線形で近似できることがわかる。具体的には、所定の区間で見た場合、時刻tにおける周波数fは、初期値f0から時刻tに比例比例係数Aして増減する線分で線形近似できると考えられる。例えば、加速時の車両では、一般的にエンジンの回転数がほぼ線形的に増加する。加速時の周波数変化を示す区間Bでは右肩上がりに周波数が増加していることが分かる。この区間ではエンジン回転数は増加しており車両は加速している。すなわち、このエンジン音の周波数は、傾きAがプラスの区分線形で近似できる。また、減速時の車両では、エンジンの回転数が線形的に減少する。減速時の周波数変化を示す区間Aでは右肩下がりに周波数が減少していることが分かる。よって、このエンジン音の周波数は、傾きAはマイナスの区分線形で近似できる。さらに、定常走行時の車両では、エンジンの回転数は変化しない。定常走行時の周波数変化を示す区間Cではほぼ一定の周波数で推移していることが分かる。この区間ではエンジン回転数は一定であり車両は定常走行している。すなわち、このエンジン音の周波数は、傾きAが0の区分線形で近似できる。そして、周波数fを上記式1で表した場合、時刻tにおける位相ψは、と表すことができる。ここで、右辺の第3項のψ0は初期位相であり、第2項2πf0tは時刻tに比例して角周波数2πf0tだけ位相が進むことを示している。そして、第1項πAt2は、位相が二次曲線で近似できることを示している。前述のように、エンジン音などの周期音では、位相の時間変化を曲線でモデル化できる。一方、風雑音などの非周期音では、位相の時間変化は周期性を有さずランダムであり、二次曲線で近似することができない。本願発明者らは、この位相の時間変化の違いに着目し、エンジン音などの周期音であって、かつ、エンジン音などのように音の周期が変化する音と、風雑音、雨音、暗騒音などの非周期音とを区別して、周期音もしくは非周期音の周波数信号を判定できることを見出した。特に、死角車両を検知するアプリケーション等では、加速時の車両を瞬時に検知することが可能となる。ここで、エンジン回転数の増減とエンジン音の位相との関係について分析する。図4aは、区間Cにおける、エンジン回転数が一定のときのエンジン音を模式的に示した図である。ここではエンジン音の周波数をfとする。図4bは基底波形を示す図である。ここでは基底波形の周波数をエンジン音の周波数fと同じ値にしている。図4cは、基底波形に対する位相を示す図である。エンジン回転数が一定であるエンジン音は図1aに示す正弦波のように一定の周期を有する。このため、図4cに示すように、所定の周波数fでの位相は、時間変化に対して加速度的に増加又は減少はしない。なお、対象とする音が一定の周波数であり、基底波形の周波数が低い場合、位相は徐々に遅れることとなる。しかし、減少量は一定となるため、位相の形状は線形的に減少することとなる。一方、対象とする音が一定の周波数であり、基底波形の周波数が高い場合、位相は徐々に早くなる。しかし、その増加量は一定となるため、位相の形状は線形的に増加することとなる。図5aは、区間Bにおける、エンジン回転数が増加して車両が加速するときのエンジン音を模式的に示した図である。このときエンジン音の周波数は時間とともに増加する。図5bは基底波形を示す図である。例えば、基底波形の周波数はfとする。図5cは、基底波形に対する位相を示す図である。エンジン音は正弦波のように周期性を有しつつ、徐々に周期が高くなる波形を有することから、図5cに示すように、基底波形に対する位相は、時間変化に対して加速度的に増加する。図6aは、区間Aにおける、エンジン回転数が減少して車両が減速するときのエンジン音を模式的に示した図である。このときエンジン音の周波数は時間とともに減少する。図6bは基底波形を示す図である。例えば、基底波形の周波数はfとする。図6cは、基底波形に対する位相を示す図である。エンジン音は正弦波のように周期性を有しつつ、徐々に周期が低くなる波形を有することから、図6cに示すように、基底波形に対する位相は、時間変化に対して加速度的に減少する。以下、本発明の実施の形態について、図面を参照しながら説明する。実施の形態1実施の形態1に係る雑音除去装置について説明する。図7及び図8は、本発明の実施の形態1における雑音除去装置の構成を示すブロック図である。図7において、雑音除去装置1500は、マイクロホン2400と、DFT分析部2402と、雑音除去処理部1504とを含む。DFT分析部2402は、請求の範囲に示す周波数分析手段に対応する。マイクロホン2400は、外部からの混合音2401を集音する。この混合音は車両のエンジン音と風雑音とから構成されている。DFT分析部2402は、入力された混合音2401に対してフーリエ変換処理を施し、混合音2401の周波数信号を複数の周波数帯域のそれぞれについて求める。なお、DFT分析部2402が行うフーリエ変換処理の代わりに、高速フーリエ変換、離散コサイン変換、又はウェーブレット変換などの別の周波数変換方法による周波数変換を用いても良い。DFT分析部2402により求められた周波数信号は、M個の周波数帯域を有する。このM個の周波数帯域を指定する番号を記号jj=1~Mで表す。雑音除去処理部1504は、位相補正部1501jj=1~Mと、抽出音判定部1502jj=1~Mと、音抽出部1503jj=1~Mとを含む。つまり、位相補正部、抽出音判定部及び音抽出部は、周波数帯域ごとに設けられている。位相補正部1501jj=1~Mは、請求の範囲に示す位相補正手段に対応する。音抽出部1503jj=1~Mは、請求の範囲に示す音響信号識別手段に対応する。位相補正部1501jj=1~Mは、M個の位相補正部から構成され、j番目の位相補正部1501jは、周波数帯域jについての処理を実行する。本明細書中において、同様の参照符号の記載を行う処理部についても、同様である。位相補正部1501jj=1~Mは、DFT分析部2402が求めた周波数帯域jの周波数信号に対して、時刻tの周波数信号の位相をψtラジアンとするときに、ψ ́t=mod2πψt-2πftfは分析周波数に位相を補正する。抽出音判定部1502jj=1~Mは、所定の時間幅において、分析の対象とする時刻の位相補正された周波数信号から位相の時間変化を近似する位相曲線近似曲線を算出し、算出された位相曲線と分析対象とする時刻の位相との誤差を算出する。このとき位相距離位相曲線と分析対象とする時刻の位相との誤差は補正後の位相ψ ́tを用いて計算する。最後に、音抽出部1503jj=1~Mは、抽出音判定部1502jj=1~Mが計算した誤差位相距離をもとに、誤差が第2の閾値以下になる周波数信号を抽出音として抽出する。これらの処理を、所定の時間幅を時間方向に移動させながら行うことにより、時間-周波数領域ごとに抽出音の周波数信号2408を取り出すことができる。図8に、抽出音判定部1502jj=1~Mの構成を示すブロック図を示す。抽出音判定部1502jj=1~Mは、周波数信号選択部1600jj=1~Mと、位相距離判定部1601jj=1~Mと、位相曲線算出部1602jj=1~Mとから構成される。位相曲線算出部1602jj=1~Mは、請求の範囲に示す誤差算出手段に対応する。周波数信号選択部1600jj=1~Mは、所定の時間幅において、位相補正部1501jj=1~Mが位相補正した周波数信号の中から位相曲線の算出及び位相距離を計算するのに用いる周波数信号を選択する。位相曲線算出部1602jj=1~Mは、周波数信号選択部1600jj=1~Mが選択した周波数信号の補正された位相ψ ́tを用いて時間経過に伴い位相が変化する位相形状を二次曲線として計算する。そして、位相距離判定部1601jj=1~Mは、位相曲線算出部1602jj=1~Mが算出した位相曲線と分析対象とする時刻の補正後の位相との位相距離を判定する。なお、本発明の必須の構成要件は、図7に示したDFT分析部2402と、音抽出部1503と、図8に示した位相距離判定部1601jと、位相曲線算出部1602jとである。DFT分析部2402が、図1cに示した本発明で定義する位相を直接導き出すことができるのであれば、位相補正部1501jは不要である。また、マイクロホン2400は、必須の構成要素ではない。次に、以上のように構成された雑音除去装置1500の動作について説明する。以下では、j番目の周波数帯域について説明を行うが、他の周波数帯域についても同様の処理が行われる。ここでは、周波数帯域の中心周波数と分析周波数とが一致する場合を例にして説明を行う。分析周波数とは、位相距離を求めるψ ́t=mod2πψt-2πftにおける周波数fである。雑音除去装置1500は、この周波数fに抽出音が存在するか否かを判定する。他の方法として、周波数帯域を含む複数の周波数を分析周波数として抽出音の判定を行ってもよい。この場合は、中心周波数の周辺の周波数に抽出音が存在するか否かを判定することができる。図9及び図10は、雑音除去装置1500の動作手順を示すフローチャートである。マイクロホン2400は、外部からの混合音2401を集音し、集音した混合音を、DFT分析部2402に出力するS200。DFT分析部2402は、混合音2401を受け付けて、混合音2401に対してフーリエ変換処理を施し、混合音2401の周波数信号を周波数帯域jごとに求めるステップS300。次に、位相補正部1501jは、DFT分析部2402が求めた周波数信号に対して、周波数帯域jごとに、位相ψを、位相ψ ́t=mod2πψt-2πftfは分析周波数に変換することで位相補正を行うステップS1700j。ここで、時刻tの周波数信号の位相をψtラジアンとする。ここで、本発明において位相を用いる理由及び位相補正を行う方法例を説明する。図3は、自動車のエンジン音をDFT分析部2402において分析したスペクトログラムである。縦軸は周波数、横軸は時間をそれぞれ示し、色の濃度は周波数信号のパワーの大きさを示している。濃い色はパワーが大きいことを示している。図3は、風などの雑音をなるべく除去したデータであり、色の濃い部分がおおむねエンジン音を示している。分析に使用したエンジン音は、時間とともに回転数が変化したデータであり、スペクトログラムから時間の経過と共に周波数が変化していることが分かる。図11は、DFT分析におけるパワーと位相について説明する図である。図11aは、図3と同様に、自動車のエンジン音をDFT分析したスペクトログラムである。図11bは時刻t1から所定の時間窓幅のハニング窓を用いて複素空間上に周波数信号601を表したものである。周波数f1、f2、f3等、各周波数のパワーと位相が算出される。周波数信号601の長さがパワーを示し、周波数信号601と実軸とのなす角が位相を示している。そして、図11aにt1、t2、t3と示すように、時間シフトを行いながら各時刻における周波数信号を求めることとなる。なお、一般的に、スペクトログラムは各時刻における各周波数のパワーを示すのみであり、位相については省略されている。図3や図11aに示すスペクトログラムも同様に、DFT分析したパワーの大きさのみを表示したものである。図11cに、図11aにおける所定の周波数例えば周波数f4での時間方向の位相の変動を示す。横軸は時間を示す。縦軸は周波数信号の位相を示しており、0~2πラジアンの間の値で示される。図11dに、図11aにおける、所定の周波数例えば周波数f4のパワーの時間変動を示す。横軸は時間軸である。縦軸は、周波数信号の大きさパワーを表している。周波数信号の位相ψt及び大きさパワーPtは、周波数信号の実部をxtと表すこととして、周波数信号の虚部をytと表すこととすると、及びである。ここで、記号tは周波数信号の時刻を表している。図12を用いて、風などの雑音がある場合の自動車のエンジン音を説明する。図12aは、図3と同様、自動車のエンジン音をDFT分析部したスペクトログラムである。縦軸は周波数、横軸は時間をそれぞれ示しており、色の濃度は周波数信号のパワーの大きさを示している。しかし、図3と異なり、風などの雑音が含まれているため、エンジン音以外の周波数にも色の濃い部分時刻t1、t2などが存在し、エンジン音なのか風雑音なのかパワーのみではまったく分からない状態となっている。図12bは、時刻t2におけるエンジン音部分のある周波数f4の所定時間のパワーの推移を示すグラフである。風雑音の影響でパワーが乱れていることが分かる。図12cは、時刻t3におけるエンジン音がない部分である周波数f4の所定時間のパワーの推移を示すグラフである。非定常のパワーが存在することが分かる。また、図12b及び図12cを比較しても、パワーのみでは風雑音なのか、エンジン音が存在するのかまったく区別がつかないことが分かる。そこで、本発明では、エンジン音を抽出すべく、位相の時間変化を用いる。まず、エンジン音の位相特性について説明する。エンジンは、所定数のシリンダーがピストン運動を行うことで駆動系を回転させている。そして、車両から発せられるエンジン音は、このエンジンの回転に依存した音と、エンジンの回転には依存しない固定振動音又は非周期音とからなる。特に、車両の外部から検知できる主な音は、エンジンの回転に依存した周期音である。本発明では、このエンジンの回転に依存する周期音をエンジン音として抽出する。図3の点線の円501、502、503に示すように、エンジン音は、回転数が変化することで周波数が変化していることが分かる。ここで周波数の変化に着目すると、周波数がランダムに変化したり、離散的に飛んだりすることはほとんどなく、所定の時間間隔でみるとほぼ時刻に比例して周波数が変化していることが分かる。したがって、エンジン音は上記式1で示すような区分線形で近似できる。具体的には所定の時間区間で見た場合、時刻tにおける周波数fは、初期値f0から時刻tに比例比例係数Aして増減する線分で線形近似できると考えられる。そして、周波数fを上記式1で表した場合、時刻tにおける位相ψは、上記式2として表すことができる。次に、位相の時間変化の近似処理を容易にするための位相補正処理について説明を行う。図11cに示されている周波数信号の位相ψtを位相ψ ́t=mod2πψt-2πftfは分析周波数に変換することで位相補正を行う。初めに、位相補正部1501jは、基準の時刻を決定する。図13aは、図11cと同じ内容のものであり、この例では、図13aの黒丸印の時刻t0を基準の時刻に決定している。次に、位相補正部1501jは、位相を補正する周波数信号の複数の時刻を決定する。この例では、図13aの5個の白丸印の時刻t1、t2、t3、t4、t5を、位相を補正する周波数信号の時刻に決定している。ここで、基準の時刻t0における周波数信号の位相をと表すこととして、位相を補正する5個の時刻における周波数信号の位相をと表すことにする。これらの補正する前の位相を図13aにおいて×印で示してある。また、対応する時刻の周波数信号の大きさはで表すことができる。次に、図14に、時刻t2における周波数信号の位相を補正する方法を示す。図14aと図13aとは、同じ内容である。また、図14bは、1/ffは分析周波数の時間間隔で、かつ等角速度で、0~2πラジアンまで規則的に変化する位相を実線で表している。ここで、補正したあとの位相をと表すことにする。図14bにおいて、基準の時刻t0と時刻t2との位相を比較する。時刻t2の位相は、時刻t0の位相よりだけ大きい。そこで、図14aにおいて、基準の時刻t0の位相ψt0との時間差に起因する位相のずれを補正するために、時刻t2の位相ψt2からΔψを差し引いてψ ́t2を求める。これが位相補正後の時刻t2の位相である。このとき、時刻t0の位相は基準の時刻における位相であるので位相補正後も同じ値となる。具体的には、位相補正後の位相をにより求める。位相補正したあとの周波数信号の位相を図13bに×印で示す。図13bの表示の方法は図13aと同様であるため説明を省略する。再度、図9を参照して、抽出音判定部1502jは、位相補正部1501jが求めた補正後の位相情報を用いて、位相の形状を算出する。そして分析の対象とする時刻の周波数信号と、分析の対象とする時刻とは異なる複数の時刻における周波数信号との位相距離誤差を求めるステップS1701j。図10は、抽出音の周波数信号を判定する処理ステップS1701jの動作手順を示すフローチャートである。初めに周波数信号選択部1600jは、位相補正部1501jが求めた所定の時間幅における位相補正された周波数信号から、位相曲線算出部1602jが位相の形状を計算する際に用いる周波数信号を選択するステップS1800j。ここでは、分析の対象とする時刻をt0として、時刻t0と時刻t1、t2、t3、t4、t5とにおける周波数信号の位相から位相の形状を算出する。このとき、位相曲線を求めるときに用いた周波数信号時刻t0~t5の6個の周波数信号は所定の値以上の数から構成されている。このことは、位相距離を求めるために選択された周波数信号の数が少ない場合に、位相の時間変化の規則性を判定することが困難になるからである。ここでの所定の時間幅の時間長は、抽出音の位相の時間変化の性質に基づいて決定することとしてもよい。次に、位相曲線算出部1602jは、位相曲線を算出するステップS1801j。位相曲線は例えば以下の二次多項式式11で近似して算出することとする。図15は位相曲線の算出処理を説明するための図である。図15に示すように、所定の数の点から二次曲線を算出することができる。本発明では二次曲線を重回帰曲線として算出する。具体的には、各時刻tii=0,1,2,3,4,5おける補正後の位相をψ ́tiとした場合、二次曲線Ψtの各係数A2、A1、A0はそれぞれ、と表せる。なお、各係数はである。再度、図10を参照して、位相距離判定部1601jは、位相曲線算出部1602jが算出した形状と分析対象とする時刻の補正後の位相との位相距離を計算するステップS1802j。この例では、位相距離誤差E0は位相の差分誤差であり、で求める。なお、分析対象とする点を除いて形状を算出し、算出した形状と分析対象とする点との位相の差を計算することとしてもよい。この計算方法によると、分析対象とする点に算出した形状から著しく外れるノイズが含まれる場合、より精度良く形状を近似することが可能となる。なお、本例では分析対象である時刻t0に対し、時刻t1、t2、t3、t4、t5の位相から位相の形状を算出している。例えば時刻t2を分析対象とする場合つまり時刻t0 ́とした場合、あらたに時刻t1 ́、t2 ́、t3 ́、t4 ́、t5 ́の位相から位相曲線を算出して誤差を算出してもよいし、既に算出されたt0、t1、t2、t3、t4、t5による位相曲線から誤差を算出することとしてもよい。すなわち、すでに算出された位相曲線を用いた誤差は、となる。この方法によると、位相曲線の算出回数が減るため、計算量の削減の効果を奏する。さらに、分析対象を所定の区間とし、誤差の平均によって、分析対象区間中の全ての周波数信号が誤差か否かを区別することとしてもよい。例えば誤差の平均は以下の式22であらわすことができる。なお、分析対象とする区間は周囲の状況に応じて可変とすることも可能である。例えば、車両が急激に加速又は減速することが多い交差点付近では分析対象とする区間を短くし、加速又は減速が比較的少ない場合は分析対象とする区間を長くすることとしてもよい。再度、図9を参照して、音抽出部1503jは、位相距離誤差が閾値以下である分析の対象とする周波数信号の各々を抽出音として抽出するステップS1702j。図16は、位相距離を求める所定の時間幅96msにおける、混合音の周波数信号の位相補正された位相ψ ́tを模式的に示した図である。横軸は時間tを示し、縦軸は位相補正された位相ψ ́tを示す。黒丸印は分析の対象とする周波数信号の位相を示し、白丸印は位相曲線を求めるのに用いられた周波数信号の位相を示す。太い破線1101は算出された位相曲線である。位相補正された各点をもとに位相曲線として二次曲線が算出されていることが分かる。細い破線1102は誤差の閾値例えば20度とするを示す。つまり、上側の破線1102は、破線1101を上に閾値分だけシフトさせたものであり、下側の破線1102は、破線1101を下に閾値分だけシフトさせたものである。分析の対象とする周波数信号の位相が2つの破線1102内に収まっていれば、その周波数信号は抽出音周期音の周波数信号であると判定され、2つの破線1102内に収まっていなければ、その周波数信号は雑音の周波数信号であると判定される。図16aにおいて、黒丸印で示される分析対象とされる周波数信号の位相は、位相の二次曲線との誤差が閾値未満である。このため、音抽出部1503jは、その周波数信号を抽出音の周波数信号として抽出する。図16bにおいて、黒丸印で示される分析対象とされる周波数信号の位相のそれぞれは、位相の二次曲線との誤差が閾値以上である。このため、音抽出部1503jは、これらの周波数信号を抽出音の周波数信号として抽出せずに、雑音として除去する。図17は、本実施の形態に示す手法によるエンジン音抽出処理を説明する図である。上記式1に示すようにエンジン音を区分線形で近似した場合、位相は上記式11に示すように二次曲線で近似できる。図17aは、図5に示したのと同じスペクトログラムである。図17b~図17eは、図17aに四角印で示す4つの領域における周波数信号を示すグラフである。4つの領域のそれぞれは、1つの周波数帯域を有する領域である。図17b~図17eに示すグラフにおいて、横軸は時間を示し、縦軸は位相を示す。白い丸印は実際の分析された周波数信号を示し、太い破線は算出された近似曲線を示す。また、細い破線は、抽出音と雑音との閾値を示す。図17bは、エンジンの回転数が下がっている、すなわち時間-周波数空間で周波数の時間変化が負の傾きの一次式で近似できるエンジン音部分の補正後の位相を示したグラフである。位相曲線は上に凸の形状をしていることが分かる。そして分析された各周波数信号はほぼ閾値以内に収まっていることが分かる。図17cは、エンジンの回転数が上がっている、すなわち時間-周波数空間で周波数の時間変化が正の傾きの一次式で近似できるエンジン音部分の補正後の位相を示したグラフである。位相曲線は下に凸の形状をしていることが分かる。そして分析された各周波数信号はほぼ閾値以内に収まっていることが分かる。図17dは、エンジンの回転数が一定、すなわち時間-周波数空間で周波数が変化しない二次の係数が0で近似できるエンジン音部分の補正後の位相を示したグラフである。位相曲線は二次の項が0となり、一次直線の形状をしていることが分かる。そして分析された各周波数信号はほぼ閾値以内に収まっていることが分かる。このことから、二次曲線による表現は周波数が変化しないエンジン音も含めて識別できることが分かる。図17eは、風雑音部分の補正後の位相を示したグラフである。風雑音の周波数信号の位相は、ばらばらであるため、二次の近似曲線を算出したとしても、その曲線からの誤差が大きく、閾値以内である信号部分はほとんどないことが分かる。このように、算出された曲線と曲線からの誤差をもとに風雑音とエンジン音を区別することができる。図18は、位相曲線からの誤差を説明する図である。横軸はエンジン音、雨音、風雑音の各音響信号を示している。縦軸は、本手法による位相曲線からの誤差の平均及び分散を示している。つまり、縦軸における線分の幅がとり得る誤差の範囲を示し、ひし形が平均値を示している。例えば、エンジン音の場合、誤差の範囲は、1度から18度の間であり、誤差の平均は10度である。分析条件は以下となる。8kHzでサンプリングされた各音声に対し、256点32msで周波数分析を行い、768点96msを区間として位相曲線の算出を行った。そして、位相曲線からの誤差の平均と分散を算出した。図18より、エンジン音は平均値が10度と位相曲線からの誤差が小さいのに対し、雨音は68度、風雑音は48度と位相の位相曲線との誤差が大きいことが分かる。このようにエンジン音のような周期音と、風雑音のような非周期音とは位相曲線からの誤差に大きな違いがあることが分かる。本例の場合、例えば閾値を20度等に設定し、閾値以下をエンジン音として適切に抽出できることとなる。図19は音識別を説明する図である。各グラフの横軸は時間を示し、縦軸は周波数を示す。図19aは風雑音とエンジン音が混合した音を周波数分析したスペクトログラムである。色の濃さはパワーの大きさを表しており、色が濃いほどパワーが大きいことを示している。分析条件は以下となる。8kHzでサンプリングされた音声に対し、512点で周波数分析を行い、1536点を区間として位相曲線の算出を行った。そして位相曲線からの誤差の閾値を20度としてエンジン音の抽出を行った。図19bは本実施の形態における手法で風雑音とエンジン音を識別したグラフである。黒い部分がエンジン音として抽出された部分である。図19aでは風などの影響により雑音が混合されているため、いったいどの部分がエンジン音であるか抽出することは困難である。しかしながら、本実施の形態における方法でエンジン音を抽出した場合、適切にエンジン音が抽出できていることを示している。特にエンジンの回転数が急激に増加する部分や、減少する部分、定常音ともに抽出できていることが分かる。なお、位相補正部1501jは、位相補正において以下に説明する位相補正の処理を、さらに行うこととしてもよい。以下に説明する位相補正を行う場合には、位相曲線の算出や位相曲線との誤差算出等の処理が付随して行われる。このため、位相補正部1501jは、抽出音判定部1502jでの計算結果を随時参照しながら処理を行うものとする。図20は、さらに実施される位相補正を説明する図である。図20のグラフはいずれもエンジン音の一部を周波数分析したグラフであり、横軸は時間を、縦軸は位相をそれぞれ示す。各白い丸印は位相補正部1501iで位相補正された周波数信号である。図20aにおいて、白い丸印で示される周波数信号の位相を用いて位相曲線を算出すると太い破線で示す曲線が算出される。細い破線は誤差閾値である。算出された位相曲線との誤差を算出すると、各周波数信号と曲線との誤差は大きく、閾値から大きく外れた点が多いことが分かる。ここで時刻t6、t7、t8、t9の周波数信号の位相に着目すると、他の時刻の位相と大きく外れていることが分かる。これは、位相が0~2πの周期でトーラス状になっていることに起因する。そこで、トーラス状に起因する現象を考慮して位相曲線を算出することとしてもよい。これにより、他の時刻の位相と大きく外れた位相を補正することができ、音の識別を精度良く行うことが可能となる。例えば、前、後、又は前後のN個の位相を用いて位相を補正することとしてもよい。図20bにおいて、例えば、時刻t1からt5例えば、N=5の位相の平均を算出する。平均位相がψ=2π×10/360であったとする。次に時刻t6の位相がψ6=2π×170/360であったとする。しかし位相はトーラスになっているため、ψ6=2π×170/360±2πの可能性がある。なお、実際は±2π×mmは自然数の可能性があるが、ここではm=1の場合のみを考慮することとする。なお、周波数が大きく変化する場合、位相変化も大きくなるため、分析する音に応じて考慮するmを可変としてもよい。また、平均を算出するための位相の選択時刻は、時刻t1~t5に限定されるものではなく、任意の時刻を用いることが可能である。次に、時刻t6の位相ψ6を平均位相ψとの誤差が小さい値に補正する。図20bの場合、ψ6=2π×170/360-2πとなる。同様に時刻t7の位相を時刻t2~t5の位相と時刻t6の補正後の位相とを用いて補正する。本例の場合ψ7=ψ7-2πと補正される。同様の処理を時刻t8、t9と行っていく。図20cに補正後の位相を示す。時刻t6、t7、t8、t9の位相が補正されていることが分かる。この補正後の位相情報を用いて位相曲線を算出した場合、太い破線で示される曲線が算出される。図20cの場合、曲線とその閾値内に各周波数信号が含まれるため、エンジン音として適切に抽出されることとなる。なお、位相の補正手法はこれに限ったものではない。例えば、まず位相曲線を算出し、算出された形状との誤差が大きい各点に±2πの位相補正を行うこととしてもよい。また、位相が取り得る角度範囲を補正することとしてもよい。以下、図を用いて説明する。図21は位相の補正処理を説明する図である。図21のグラフはいずれも縦軸が位相、横軸が時間をそれぞれ示す。白い丸印は各時刻における周波数信号の位相を示す。図21aは0から2πを角度範囲とした場合の周波数信号の位相を示す。各位相を元に位相曲線が算出され黒い曲線で示している。図21cは曲線との誤差を元に位相を補正している。具体的には時刻t1の位相に2πを加算する補正がされている。また時刻t8の位相に-2πを加算する補正がされている。一方、図21bは-πからπを角度範囲とした場合の周波数信号の位相を示す。図21aと同様に、各位相を元に位相曲線が算出され黒い曲線で示している。図21dは曲線との誤差を元に位相を補正している。具体的には時刻t10の位相に-2πを加算する補正がされている。図21cの角度範囲の場合の曲線との誤差と、図21dの角度範囲の場合の曲線との誤差を比較した場合、図21cの角度範囲の場合の曲線との誤差が小さくなる。そこで図21cの角度範囲を用いた位相曲線を用いる。このように、角度範囲の制御を行い位相曲線を算出することとしてもよい。これにより、他の時刻の位相と大きく外れた位相を補正することができ、より音の識別を精度良く行うことが可能となる。以上説明したように、本実施の形態によると、時間-周波数領域ごとにエンジン音や音声などの周期音と、風雑音、雨音、暗騒音などの非周期音とを区別して、周期音もしくは非周期音の周波数信号を判定できる。特にエンジン音などのように時間-周波数空間で周波数が時間とともに変化する周期音を正確に判定することができる。特に死角車両を検知するアプリケーション等において、衝突すると大きな事故につながりやすい加速している車両を正確に検知することができる。実施の形態2次に、実施の形態2に係る車両検知装置について説明する。実施の形態2に係る車両検知装置は、複数のマイクロホンから入力される各々の混合音から、エンジン音抽出音の周波数信号を判定し、音の到達時間差より車両の到達方向を算出し、運転者に接近車両の存在を知らせるものである。図22及び図23は、本発明の実施の形態2における車両検知装置の構成を示すブロック図である。図22において、車両検知装置4100は、マイクロホン41071と、マイクロホン41072と、DFT分析部1100と、車両検知処理部4101と、方向検知部4108とを含む。車両検知処理部4101は、位相補正部4102jj=1~Mと、抽出音判定部4103jj=1~Mと、音抽出部4104jj=1~Mとを含む。また、図23において、抽出音判定部4103jj=1~Mは、位相距離判定部4200jj=1~Mと、位相曲線算出部4201jj=1~Mと、周波数信号選択部4202jj=1~Mとを含む。図22において、マイクロホン41071は、外部からの混合音24011を集音する。マイクロホン41072は、外部からの混合音24012を集音する。この例では、マイクロホン41071とマイクロホン41072はそれぞれ自車両の左前と右前のバンパーに設置されている。これらの混合音の各々は、例えば8kHzでサンプリングされた車両のエンジン音と風雑音とから構成されている。なお、サンプリング周波数は8kHzに限定されるものではない。DFT分析部1100は、入力された混合音24011と混合音24012の各々に対して離散フーリエ変換処理を施し、混合音24011と混合音2の周波数信号を求める。ここで、DFTの時間窓幅は256点38msである。以下では、DFT分析部1100から求められた周波数帯域の個数をMとして、それらの周波数帯域を指定する番号を記号jj=1~Mで表すこととする。この例では、車両のエンジン音が存在する10Hz~500Hzの周波数帯域を10Hz間隔ごとに分割してM=50、周波数信号を求める。位相補正部4102jj=1~Mは、DFT分析部1100が求めた周波数帯域jj=1~Mの周波数信号に対して、時刻tの周波数信号の位相をψtラジアンとするときに、ψ ́ ́t=mod2πψt-2πf ́tf ́は周波数帯域の周波数に位相を補正する。なお、この例は、ψtを分析周波数で補正するのではなく、周波数信号を求めた周波数帯域の周波数f ́で補正を行う。抽出音判定部4103jj=1~Mは、所定の時間幅において、分析の対象とする時刻の位相補正された周波数信号から位相曲線を算出し、算出された位相曲線をもとに抽出音の判定を行う。このとき、位相距離を求めるときに用いた周波数信号の数は第1のしきい値以上の数から構成されている。ここでは所定の時間幅を96msとする。なお、このとき位相距離は補正後の位相ψ ́ ́tを用いて計算する。抽出音判定部4103jj=1~Mが実施する処理は、実施の形態1に示した抽出音判定部1502jj=1~Mが実施する処理と同様である。このため、その詳細な説明は繰り返さない。図23に、抽出音判定部4103jj=1~Mの構成を示すブロック図を示す。抽出音判定部4103jj=1~Mは、位相距離判定部4200jj=1~Mと、位相曲線算出部4201jj=1~Mと、周波数信号選択部1600jj=1~Mとから構成される。周波数信号選択部4202jj=1~Mは、所定の時間幅において、位相補正部4102jj=1~Mが位相補正した周波数信号の中から位相曲線の算出及び位相距離誤差を計算するのに用いる周波数信号を選択する。周波数信号選択部4202jj=1~Mが実施する処理は、実施の形態1に示した周波数信号選択部1600jj=1~Mが実施する処理と同様である。このため、その詳細な説明は繰り返さない。位相曲線算出部4201jj=1~Mは、周波数信号の補正された位相ψ ́ ́tを用いて時間経過に伴い位相が変化する位相形状を曲線として計算する。位相曲線算出部4201jj=1~Mが実施する処理は、実施の形態1に示した位相曲線算出部1602jj=1~Mが実施する処理と同様である。このため、その詳細な説明は繰り返さない。そして位相距離判定部4200jj=1~Mは、位相曲線算出部4201jj=1~Mが算出した位相曲線との位相距離が第2のしきい値以下か否かを判定する。具体的には位相曲線を算出する区間を768点96msとして位相曲線を算出し、位相距離を算出する。位相距離判定部4200jj=1~Mによる位相曲線の算出方法と位相距離誤差の算出方法とは、実施の形態1に示した位相距離判定部1601jj=1~Mのそれらと同様である。このため、その詳細な説明は繰り返さない。次に、音抽出部4104jj=1~Mは、抽出音判定部4103jj=1~Mが判定した位相距離をもとにエンジン音を抽出する。具体的には、誤差の閾値を20度とし、閾値以下をエンジン音として抽出する。音抽出部4104jj=1~Mが実施する処理は、実施の形態1に示した音抽出部1503jj=1~Mが実施する処理と同様である。このため、その詳細な説明は繰り返さない。なお、音抽出部4104jj=1~Mは、さらに、エンジン音が抽出されたときに、抽出音検知フラグ4105を出力する。再度、図22を参照して、方向検知部4108は、音抽出部4104jj=1~Mから抽出音検知フラグ4105が出力されたときに、音抽出部4104jj=1~Mで抽出されたエンジン音の時間-周波数領域に対し、車両の存在する方向を特定する。例えば、その領域におけるエンジン音の到達時間差をもとに車両の方向を検知する。例えばいずれか1つのマイクロホンにおいてエンジン音が抽出された場合、両マイクロホンを用いて車両の存在する方向を特定する。風雑音は両マイクロホンに対して均一ではなく、一方のマイクロホンにのみ風雑音が存在し、他方には存在しない場合もあるからである。なお、両マイクロホンにおいてエンジン音が抽出された場合にのみ、方向を特定することとしてもよい。マイクロホン41071及びマイクロホン41072の間隔をdmとする。エンジン音が自車両に対して方角θラジアンから検出されるとする。マイク間での到達時間差をΔtsとし、音速をcm/sとすると、方向θラジアンは以下の式23で表すことができる。最後に、車両検知装置4100に接続された提示部4106は、方向検知部4108で検知された車両の方向を運転者に知らせる。例えば、提示部4106は、どちらの方向から車両が来ているのかをディスプレイに表示するようにしても良い。車両検知装置4100及び提示部4106は、これらの処理を、所定の時間幅を時間方向に移動させながら行う。次に、以上のように構成された車両検知装置4100の動作について説明する。以下では、j番目の周波数帯域周波数帯域の周波数はf ́について説明を行う。図24及び図25は、車両検知装置4100の動作手順を示すフローチャートである。初めに、マイクロホン41071及び41072は、外部からの混合音2401をそれぞれ集音し、集音した混合音を、DFT分析部2402に出力するステップS201。DFT分析部1100は、混合音24011と混合音24012を受け付けて、混合音24011と混合音24012のそれぞれに対して離散フーリエ変換処理を施し、混合音24011と混合音24012の周波数信号を求めるステップS300。次に、位相補正部4102jは、DFT分析部1100が求めた周波数帯域j周波数f ́の周波数信号に対して、時刻tの周波数信号の位相をψtラジアンとするときに、位相ψtを位相ψ ́ ́t=mod2πψt-2πf ́tf ́は周波数帯域の周波数に変換することで位相補正を行うステップS4300j。次に、抽出音判定部4103j位相距離判定部4200jは、混合音混合音24011、混合音24012ごとに、所定の時間幅における第1のしきい値以上の数から構成されている位相補正された周波数信号第1のしきい値は、所定の時間幅における時刻の周波数信号の80の数であるの位相ψ ́ ́tを用いて、分析周波数fを設定して、設定された分析周波数fを用いて位相距離を求めるステップS4301j。図25を用いて、ステップS4301jの処理について詳細に説明する。初めに周波数信号選択部4202jは、位相補正部4102jが求めた所定の時間幅における位相補正された周波数信号から、位相曲線算出部4201jが位相の形状の計算に用いる周波数信号を選択するステップS1800j。そして、位相曲線算出部4201jが、位相曲線を算出するステップS1801j。次に、位相距離判定部4200jは、位相曲線算出部4201jが算出した形状と分析対象とする時刻の補正後の位相との位相距離を計算するステップS1802j。再度、図24を参照して、音抽出部4104jは、位相距離が第2のしきい値以下になる所定の時間幅における周波数信号をエンジン音の周波数信号に判定するステップS4302j。方向検知部4108は、音抽出部4104jで抽出されたエンジン音の時間-周波数領域に対し、車両の存在する方向を特定し、提示部4106は、方向検知部4108で検知された車両の方向を運転者に知らせるステップS4303。以上説明したように、実施の形態2に係る車両検知装置によると、エンジン音が抽出された場合に、エンジン音の到達時間差をもとに車両の方向を検知している。このため、雑音の影響を受けずに車両の方向を正確に検知することができる。以上、本発明の実施の形態に係る雑音除去装置及び車両検知装置について説明したが、本発明は、これらの実施の形態に限定されるものではない。例えば、上述の実施の形態ではエンジン音の抽出を例に挙げて説明を行ったが、本発明が抽出対象とする音はエンジン音に限定されるものではなく、例えば、人間もしくは動物の音声又はモーター音などのように周期音であれば本発明を適用可能である。また、音抽出部は、周波数信号ごとに周期音か雑音かを判断したが、所定の時間幅ごとにその時間幅に含まれる周波数信号が周期音か雑音かを判断するようにしても良い。例えば、図16を参照して、音抽出部は、所定の時間幅ごとに、その時間幅に含まれる周波数信号の位相に対する、位相曲線算出部が求めた二次曲線との誤差が閾値未満となる位相の割合が所定の割合以上となる場合に、その時間幅に含まれる周波数信号の全てが周期音と判断し、上記割合が所定の割合未満となる場合に、その時間幅に含まれる周波数信号の全てが雑音と判断しても良い。また、上記の各装置は、具体的には、マイクロプロセッサ、ROM、RAM、ハードディスクドライブ、ディスプレイユニット、キーボード、マウスなどから構成されるコンピュータシステムとして構成されても良い。RAM又はハードディスクドライブには、コンピュータプログラムが記憶されている。マイクロプロセッサが、コンピュータプログラムに従って動作することにより、各装置は、その機能を達成する。ここでコンピュータプログラムは、所定の機能を達成するために、コンピュータに対する指令を示す命令コードが複数個組み合わされて構成されたものである。さらに、上記の各装置を構成する構成要素の一部又は全部は、1個のシステムLSILargeScaleIntegration:大規模集積回路から構成されているとしても良い。システムLSIは、複数の構成部を1個のチップ上に集積して製造された超多機能LSIであり、具体的には、マイクロプロセッサ、ROM、RAMなどを含んで構成されるコンピュータシステムである。RAMには、コンピュータプログラムが記憶されている。マイクロプロセッサが、コンピュータプログラムに従って動作することにより、システムLSIは、その機能を達成する。さらにまた、上記の各装置を構成する構成要素の一部又は全部は、各装置に脱着可能なICカード又は単体のモジュールから構成されているとしても良い。ICカード又はモジュールは、マイクロプロセッサ、ROM、RAMなどから構成されるコンピュータシステムである。ICカード又はモジュールは、上記の超多機能LSIを含むとしても良い。マイクロプロセッサが、コンピュータプログラムに従って動作することにより、ICカード又はモジュールは、その機能を達成する。このICカード又はこのモジュールは、耐タンパ性を有するとしても良い。また、本発明は、上記に示す方法であるとしても良い。また、これらの方法をコンピュータにより実現するコンピュータプログラムであるとしても良いし、前記コンピュータプログラムからなるデジタル信号であるとしても良い。さらに、本発明は、上記コンピュータプログラム又は上記デジタル信号をコンピュータ読み取り可能な記録媒体、例えば、フレキシブルディスク、ハードディスク、CD-ROM、MO、DVD、DVD-ROM、DVD-RAM、BDBlurayDisc登録商標、半導体メモリなどに記録したものとしても良い。また、これらの記録媒体に記録されている上記デジタル信号であるとしても良い。また、本発明は、上記コンピュータプログラム又は上記デジタル信号を、電気通信回線、無線又は有線通信回線、インターネットを代表とするネットワーク、データ放送等を経由して伝送するものとしても良い。また、本発明は、マイクロプロセッサとメモリを備えたコンピュータシステムであって、上記メモリは、上記コンピュータプログラムを記憶しており、上記マイクロプロセッサは、上記コンピュータプログラムに従って動作するとしても良い。また、上記プログラム又は上記デジタル信号を上記記録媒体に記録して移送することにより、又は上記プログラム又は上記デジタル信号を上記ネットワーク等を経由して移送することにより、独立した他のコンピュータシステムにより実施するとしても良い。さらに、上記実施の形態及び上記変形例をそれぞれ組み合わせるとしても良い。今回開示された実施の形態は全ての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなくて請求の範囲によって示され、請求の範囲と均等の意味及び範囲内での全ての変更が含まれることが意図される。"