周期的な波形要素を含む音波の変動強度を推定する音質評価システムにおいて、前記音波の振幅包絡を周波数解析して最もスペクトルの大きい基本周波数を算出する周波数解析手段と、当該算出された基本周波数の逆数を算出して前記波形要素の周期を求め、当該波形要素の周期を用いて波形要素の振幅が最大となる時刻であるオンセット時刻、振幅が最小となる時刻であるオフセット時刻、スペクトル重心、音響パワー、正規化された音響パワーのうち少なくとも一つを含む要素特徴量を求める要素特徴量算出手段と、当該算出された要素特徴量を用いて変動強度を推定する変動強度推定手段と、を備えたことを特徴とする音質評価システム。請求項1における音質評価システムにおいて、前記算出された波形要素の1/4周期ごとに振幅包絡の傾きを求め、当該傾きの符号が変化する時刻間の長さである時間領域を求めて当該時間領域内で振幅が最大となる時刻を探索し、当該時刻における振幅の逸脱量である振幅逸脱量を算出する振幅逸脱量算出手段を備え、前記変動強度推定手段が、当該振幅逸脱量および前記要素特徴量を用いて変動強度を推定するものである音質評価システム。請求項2における音質評価システムにおいて、前記探索された時刻ごとの間隔を抽出し、各時間間隔の逸脱量である時間逸脱量を算出する時間逸脱量算出手段を備え、前記変動強度推定手段が、当該時間逸脱量および前記要素特徴量を用いて変動強度を推定するものである音質評価システム。周期的な波形要素を含む音波の変動強度を推定する音質評価方法において、前記音波の振幅包絡を周波数解析して最もスペクトルの大きい基本周波数を算出するステップと、当該算出された基本周波数の逆数から前記波形要素の周期を求め、当該波形要素の周期を用いて波形要素の振幅が最大となる時刻であるオンセット時刻、振幅が最小となる時刻であるオフセット時刻、スペクトル重心、音響パワー、正規化された音響パワーのうち少なくとも一つを含む要素特徴量を求めるステップと、当該算出された要素特徴量を用いて変動強度を推定するステップと、を備えたことを特徴とする音質評価方法。請求項4における音質評価方法において、前記算出された波形要素の1/4周期ごとに振幅包絡の傾きを求め、当該傾きの符号が変化する時刻間の長である時間領域を求めて当該時間領域内で振幅が最大となる時刻を探索し、当該時刻における振幅の逸脱量である振幅逸脱量を算出するステップを備え、当該振幅逸脱量および前記要素特徴量を用いて変動強度を推定するようにした音質評価方法。請求項5における音質評価方法において、前記探索された時刻ごとの間隔を抽出し、各時間間隔の逸脱量である時間逸脱量を算出するステップを備え、当該時間逸脱量および前記要素特徴量を用いて変動強度を推定するようにした音質評価方法。このような排気音を発する場合、一般的には、その音の大きさラウドネスや音の鋭さシャープネスなどを大きくすればよいが、このように音の大きさや鋭さを大きくすると、逆に周囲に対する騒音が大きくなってしまう。一方、エンジンの排気音などのような振幅変動音を発するとともに、その変調周波数を4Hz~8Hz程度にすると、最も歩行者に刺激を与えることができる。しかしながら、従来の変動強度の推定方法を用いて最適な振幅変動音を評価する場合、発せられる衝撃音の波形形状が異なっている場合であっても、これを評価することができなかった。すなわち、従来における変動強度の推定方法では、その音のラウドネスの時間履歴にどれだけ変動成分が含まれているかを調べるだけであるため、断続的に発せられるそれぞれの音のエンベロープの立ち上がりが速い場合と遅い場合を比較した場合、実際には、立ち上がりの速い音の方が変動感が強く感じられるにもかかわらず、これらの違いをうまく評価することができなかった。また、4サイクルエンジン音の場合、吸気、圧縮、爆発、排気の4つのサイクルを繰り返すため、排気音もこれらの4つのサイクルに合わせた波形要素を含むことになるが、実際には、爆発音の最大振幅がずれたり、あるいは、発音時刻がずれたりする場合がある。これらの違いは、人間に変動感の違いとして感じられるが、従来の方法では、これを大きく評価に組み込むことができなかった。そこで、本発明は、上記課題に着目して、周期的な波形要素からなる音の音質を評価する場合、その波形要素の性質を考慮して、実際の変動感に近づけられるようにした音質評価システムを提供することを目的とする。すなわち、本発明は上記課題を解決するために、周期的な波形要素を含む音波の変動強度を推定する音質評価システムにおいて、前記音波の振幅包絡を周波数解析して最もスペクトルの大きい基本周波数を算出する周波数解析手段と、当該算出された基本周波数の逆数を算出して前記波形要素の周期を求め、当該波形要素の周期を用いて波形要素の振幅が最大となる時刻であるオンセット時刻、振幅が最小となる時刻であるオフセット時刻、スペクトル重心、音響パワー、正規化された音響パワーのうち少なくとも一つを含む要素特徴量を求める要素特徴量算出手段と、当該算出された要素特徴量を用いて変動強度を推定する変動強度推定手段とを備えるようにしたものである。このようにすれば、周期的な波形要素の特徴を考慮することによって、振幅変動の値に波形要素の値を組み込むことができ、より人間が聴いたときの変動感に近づけることができるようになる。また、このような発明において、波形要素の1/4周期ごとに振幅包絡の傾きを求め、当該傾きの符号が変化した時刻間の長さである時間領域を求めて当該時間領域内で振幅が最大となる時刻を探索し、当該時刻における振幅の逸脱量である振幅逸脱量を算出する振幅逸脱量算出手段を備え、当該振幅逸脱量および前記要素特徴量を用いて変動強度を推定するようにする。このようにすれば、エンジンの不規則な爆発音の大きさのずれなども考慮して変動強度を求めることができるようになる。さらには、前記探索された時刻ごとの間隔を抽出し、各時間間隔の逸脱量である時間逸脱量を算出する時間逸脱量算出手段を備え、当該時間逸脱量および前記要素特徴量を用いて変動強度を推定するようにする。このようにすれば、エンジンの不規則な爆発音の時間のずれなども考慮して変動強度を求めることができるようになる。本発明によれば、周期的な波形要素を含む音波の変動強度を推定する音質評価システムにおいて、前記音波の振幅包絡を周波数解析して最もスペクトルの大きい基本周波数を算出する周波数解析手段と、当該算出された基本周波数の逆数から前記波形要素の周期を求め、当該算出された要素周期を用いて波形要素の振幅が最大となるオンセット時刻、振幅が最小となるオフセット時刻、スペクトル重心、音響パワー、正規化された音響パワーのうち少なくとも一つを含む要素特徴量を求める要素特徴量算出手段と、当該算出された要素特徴量を用いて変動強度を推定する変動強度推定手段とを備えるようにしたので、振幅変動の値に波形要素の値を組み込むことができ、より人間が聴いたとき変動感に近づけることができるようになる。以下、本発明の一実施の形態について図面を参照しながら説明する。この実施の形態における音質評価システム100は、パーソナルコンピューターなどによって構成されるものであって、マイクを介して入力された音の振幅変動を算出して出力するようにしたものである。より具体的には、音波の入力を受け付ける入力受付手段1と、その受け付けた音の振幅包絡を算出する振幅包絡算出手段2と、その算出された振幅包絡を周波数解析して基本周波数を算出する周波数解析手段3とを備えている。そして、特徴的には、その基本周波数の逆数から波形要素の周期を求め、その周期ごとに波形要素を推定して、その波形要素のオンセット時刻、オフセット時刻、スペクトル重心、音響パワー、正規化された音響パワーなどを求めて変動強度を求めるようにしたものである。以下、本実施の形態における音質評価システム100について詳細に説明する。まず、入力受付手段1は、評価の対象となる音源からの音波の入力を受け付ける。ここで入力される音波としては、バイクの排気音を想定する。この排気音は、エンジン音を構成する音の一つであり、シリンダー内で燃料と空気を混合させて爆発させる衝撃波の一種である。通常、エンジンは1秒間に何回も爆発や排気などを繰り返すため、その排気音も断続的なものになる。一般に、インジェクターを搭載したエンジンでは、燃料と空気をコンピューターで電子的に制御して混合させるのに対して、キャブレターを搭載したエンジンは、吸気口内の気圧が下がることによって燃料と空気を混合させるようになっており、このようなエンジンでは、その排気音の発音時刻の間隔や爆発音の大きさが不均一になって、図2に示すような複数の波形要素を含む波形形状となる。図2aのうち、縦方向の太い破線で示された符号11は波形要素の発音時刻、横線で示される符号12は発音時刻における振幅、符号13は発音時刻の間隔InterOnsetInterval以下、IOIと称するである。図2aにより、単一の爆発音である波形要素の最大振幅および発音時刻は不均一であり、等間隔なIOIからの逸脱や爆発音の最大振幅の逸脱が含まれていることがわかる。それらの逸脱量は、知覚印象に大きく影響を与えるため、変動感を音響波形から推定する場合は、それらを考慮する必要がある。また、排気音は、波形要素である爆発音の繰り返しで構成されているため、その繰り返し回数が変動感に影響を及ぼすと考えられる。また、排気音は車種によっても大きく変わるため、図3に示すような排気音の波形要素から、スペクトル重心や音響パワー、正規化された音響パワー、立ち上がり時間および立ち下がり時間の長さ、立ち上がり時間および立ち下がり時間における傾きやその高さなどを求め、変動強度を推定するものとする。振幅包絡算出手段2は、この入力された波形を平滑化して振幅包絡を算出する。この振幅包絡を算出する場合は、まず、入力された音響波形を二乗してパワー関数を算出し、その後、所定の窓幅の移動平均を数回繰り返して算出する。すると、図2bに示すような平滑化された振幅包絡が得られる。次に、周波数解析手段3は、このようにして得られた振幅包絡を周波数解析する図2c。この周波数解析を行う場合は、全波整流によって得られた振幅包絡に対してFastFourierTransformを行い、得られた周波数から最もスペクトルの大きな周波数である基本周波数を抽出する。この得られた基本周波数は、入力された音響波形に含まれる波形要素の周波数とほぼ等しい周波数となる。一次変動算出手段4は、このようにして得られた値を用いて変動成分を抽出する。この変動成分を抽出する場合は、周波数解析手段3によって解析された周波数に対して、1Hz~20Hzを通す重み付き帯域通過フィルタを通し、そのフィルタに通した波形に対し逆フーリエ変換InverseFastFourierTransformを行う。そして、逆フーリエ変換した後の波形の実効値RMS値を求めて、そのRMS値を一次変動パラメタの値とする。二次変動算出手段5は、各波形要素の最大振幅やIOIがどれだけ逸脱しているかを算出する。この二次変動の値を算出する場合、各波形要素のオンセット時刻およびオフセット時刻を推定する。このオンセット時刻やオフセット時刻の推定においては、同様に、音響波形を二乗して得られたパワー関数に所定の窓幅の移動平均を数回行い、得られた波形の基本周波数を求めた後にそれの逆数を算出して波形要素の周期を求める。次に、その波形要素の1/4周期ごとにその傾きを最小二乗法を用いて推定し、その傾きの符号が正から負に変化した場合の長さ、すなわち1/2周期の長さを持つ時間領域を求め、その領域内で振幅が最大となる時刻をオンセット時刻とし、また、傾きの符号が負から正に変化した場合の長さ、すなわち1/2周期の長さを持つ時間領域を求め、その領域内で振幅が最小となる時刻をオフセット時刻とする。次に、得られたオンセット時刻に対して、逸脱がない所定の状態のオンセット時刻からのずれを抽出し、時間逸脱量ODkとした。ここでkは波形要素のID、Nは波形要素の総数を示すものとする。なお、k1の波形要素におけるオンセット時刻を基準とし、各波形要素におけるODkを算出している。また、オンセット時刻における音圧レベルを音圧レベルの最大値からのずれとして振幅逸脱量ADkとした。そして、得られた各逸脱量xkに対して移動平均を行い、傾向曲線を求め、次の<数1>~<数5>を用いて表1における特徴量pijを求める。このように算出される一次変動と二次変動の特徴量をまとめると図8のようになり、排気音の変動成分や、時間逸脱量や振幅逸脱量における標準偏差や傾向からの差分量の和、傾向に対する全体の変動幅、傾向の階差の輪、逸脱がない状態における値からの誤差の和などを求めることができる。三次変動算出手段6は、音響波形から波形要素を切り出した後に、各波形要素のスペクトル重心、音響パワー、正規化音響パワー、立ち上がり時間および立ち下がり時間の長さ、立ち上がり時間および立ち下がり時間における傾きやその高さを算出する。この三次変動の値を算出する場合、まず、二次変動の算出の場合と同様に、オンセット時刻およびオフセット時刻を抽出し、これらの時刻を用いて波形要素を切り出す。そして、その切り出された波形要素から、スペクトル重心、音響パワー、正規化音響パワー、立ち上がり時間および立ち下がり時間の長さ、立ち上がり時間および立ち下がり時間における傾きやその高さの9種類の特徴量を算出する。スペクトル重心を算出する場合は<数6>、音響パワーの算出をする場合は<数7>、正規化音響パワーを算出する場合は<数8>、立ち上がり時間における高さの算出を行う場合は<数9>、立ち下がり時間における高さを算出する場合は<数10>を用いる。また、各傾きは最小二乗法を用いて得られた一次回帰式の傾きであり、その傾きを用いて各高さを算出する。そして、1つの波形要素に対して算出した各特徴におけるそれらの平均、変動幅、標準偏差および絶対値の総和を計算し、図9に示すように、9種類の特徴×4=36個の特徴量を抽出する。変動強度推定手段7は、このように算出された一次変動の特徴量から三次変動の特徴量までを用いて最終的な変動強度を推定する。この変動強度を推定する場合、抽出した一次変動の特徴量から三次変動までの特徴量について主成分分析を行い、それらの特徴量を無相関となる特徴量に修正する。次に、実験によって得られた主観評価値と修正後の特徴量を用いて重回帰分析を行うことで変動感を推定する重回帰式を求めるようにする。次に、このように構成された音質評価システム100を用いて一次変動の特徴量から三次変動までの特徴量を抽出し、入力された音の変動強度を求める方法について説明する。まず、一次変動の特徴量を抽出する場合は、図4に示すように、入力された音響波形に対してステップS1全波整流を行って振幅包絡を算出しステップS2、周波数解析であるFFTを行うステップS3。次に、1Hz~20Hzを通す重み付き帯域通過フィルタを通しステップS4、そのフィルタに通した波形に対して逆フーリエ変換を行いステップS5、IFFT後の波形の実効値RMS値である一次変動の特徴量を求めるステップS6。次に、二次変動の特徴量を求める場合は、図5に示すように、入力された音響波形ステップT1を二乗して得られたパワー関数にステップT2、例えば、窓幅20msの移動平均を3回行って振幅包絡を算出しステップT3、その得られた振幅包絡の波形を周波数解析して基本周波数を算出するステップT4。そして、その基本周波数の逆数から波形要素の周期を求めステップT5。次に、排気音の1/4周期ごとにその傾きを最小二乗法を用いて推定し、その傾きの符号が正から負に変化した場合の長さ、すなわち1/2周期の長さを持つ時間領域を求め、その領域内で振幅が最大となる時刻をオンセット時刻として抽出し、また、傾きの符号が負から正に変化した場合の長さ、すなわち1/2周期の長さを持つ時間領域を求め、その領域内で振幅が最小となる時刻をオフセット時刻として抽出するステップT6。次に、得られたオンセット時刻に対して逸脱がない状態のオンセット時刻からのずれを抽出し、時間逸脱量を算出するステップT7。また、オンセット時刻における音圧レベルを、音響波形における音圧レベルの最大値からのずれとして音響波形から抽出し、振幅逸脱量を抽出するステップT8。そして、このように得られた逸脱量に対して移動平均を行って傾向曲線を求め、上記<数1>~<数5>を用いて二次変動の特徴量を求めるステップT9。三次変動の特徴量を求める場合は、図6に示すように、二次変動の特徴量の算出の場合と同様に、各波形要素のオンセット時刻およびオフセット時刻を抽出しステップU1、その各時刻を用いて音響波形から波形要素を切り出すステップU2。そして、その切り出された波形要素におけるスペクトル重心、音響パワー、正規化された音響パワー、立ち上がり時間の長さや立ち下がり時間の長さ、立ち上がり時間および立ち下がり時間における傾きやその高さの9種類の特徴量を<数6>~<数10>を用いて算出するステップU3。そして、このように算出された各特徴量における平均、変動幅、標準偏差および絶対値の総和を算出し、三次変動の特徴量を抽出するステップU4。次に、このように一次変動から三次変動の特徴量を算出した後、図7に示すように、各特徴量を主成分分析ステップSTU1を行うことによって無相関となる特徴量に修正し、重回帰式による推定を行った後ステップSTU2、変動強度の推定を行うステップSTU3。このように上記実施の形態によれば、周期的な波形要素を含む音波の変動強度を推定する音質評価システム100において、当該音波の振幅包絡を周波数解析して最もスペクトルの大きい基本周波数を算出する周波数解析手段3と、当該算出された基本周波数の逆数から前記波形要素の周期を求め、当該要素周期を用いて波形要素のオンセット時刻、オフセット時刻、スペクトル重心、音響パワー、正規化された音響パワーのうち少なくとも一つを含む特徴量を求める三次変動算出手段6と、当該算出された要素特徴量を用いて変動強度を推定する変動強度推定手段7とを備えるようにしたので、振幅変動の値に波形要素の値を組み込むことができ、より人間が聴いた際の変動感に近づけることができるようになる。なお、上記実施の形態では、振幅変動を有する音波の変動強度を推定する場合について説明したが、例えば、救急自動車のサイレンなどのように周波数変動を有する音波の変動強度も推定する場合にも適用することができる。このような場合としては、周波数変動音を周波数解析して基本周波数を算出し、その基本周波数を有する時間関数を上記実施の形態における振幅包絡とみなして同様の処理をすることができる。この場合、スペクトル重心は求めずに、音響パワー、正規化音響パワー、立ち上がり時間および立ち下がり時間の長さ、立ち上がり時間および立ち下がり時間における傾きやその高さの合計8個のパラメータを計算する。なお、音響パワーは基本周波数の時間関数から1周期の波形の面積を用いて求めるものとする。一次変動と二次変動を含む排気音に対する変動感の主観評価実験を行うためには、様々な二次変動を持つ排気音が必要となるが、収録した現実の排気音が持つ二次変動は網羅的でないという問題がある。そこで、収録した排気音から時間および振幅の逸脱を抽出した後に、得られた各逸脱を強調または抑制することにより、様々な二次変動を持つ排気音を合成する。排気音は、爆発音の繰り返しによって構成されるので、記録した排気音から爆発音を抽出し、それを時間軸上でつなぎ合わせることで排気音を合成する。しかし、そのような方法で合成した排気音はオートバイの音色を持っているため、変動感のみだけでなく、別の印象や感覚を感じてしまう可能性がある。そこで、オートバイの音色をもたない排気音を合成するために、記録した排気音から抽出した爆発音をそのまま用いるのではなく、その音の振幅包絡をもつピンクノイズを用いるようにした。合成した音を用いて、防音室にて正常な聴力を持つ健聴者5名に対して主観評価実験を行った。用いた刺激音は上述の排気音とする。なお、用いた二次変動の傾向は、キャブレターを搭載したエンジン駆動車3台、インジェクターを搭載したエンジン駆動車3台の計6台に含まれるものであり、逸脱量を強調または抑制する割合は、0、150、300の3通りとし、合成した音は、時間および振幅逸脱量において強調する割合を組み合わせた9通りの方法によって強調された音である。各オートバイの排気音は、量子化ビット数16bit、サンプリング周波数48kHzで録音された排気音である。主観評価実験では、シェッフェの一対比較法の浦変法を用いた。聴取者には、2つの音を提示し、最初に聞いた音に比べて後に聞いた音がどの程度変動しているかを7段階評価で回答させた。抽出した特徴量と、実験によって得られた主観評価値に対して重回帰分析を行い、FS推定で用いる重回帰式を求める。重回帰分析では、目的変数を主観評価値、説明変数を抽出したパラメタとする。重回帰分析において、複数の説明変数間に強い相関があると、偏回帰係数の推定量が不安定になるため、求められた予測精度の信頼性が低くなる。そこで、抽出した特徴に対して主成分分析を行い、それらの特徴を無相関な変数に修正し、それらを用いて重回帰分析を行う。なお、新たに用いる説明変数は累積寄与率がほぼ100までの第5主成分までとする。回帰分析によって得られた回帰方程式の予測精度は決定係数R2で表され、決定係数が08以上でかなり良い精度、決定係数が05以上且つ08未満で良い精度であると言える。重回帰分析では、説明変数を増やせば増やすほど決定係数が高くなると言われており、そのような決定係数の欠点を補うために定義されている自由度調整済み決定係数adjustedR2を算出する。1次変動の特徴のみで変動感を推定した場合のR2と、抽出した全ての特徴を考慮に入れた場合のadjustR2を比較した。その結果、抽出したすべての特徴を考慮に入れた場合の方が、一次変動の特徴のみで変動感を推定した場合よりも予測精度が高いことが確認された。