(57)【要約】【課題】操作性が良好で且つ波形のひずみのほとんどない低周波音測定を実現するために、低周波音測定器のフィルタにより生じた波形のひずみを取り除いて元の原波形を復元する波形復元装置を提供する【解決手段】波形復元装置１は、データ記録装置１１、デジタル逆フィルタ１２及び制御部１３から構成される。データ記録装置１１は、低周波音圧レベル計３により測定した測定データを記録している。デジタル逆フィルタ１２は、低周波音圧レベル計３に内蔵されているハイパスフィルタ３３と逆の周波数特性を有する逆フィルタであり、後述するデジタル逆フィルタ設計装置２により設計される。制御部１３は、データ記録装置１１及びデジタル逆フィルタ１２を制御し、データ記録装置１１からの出力にデジタル逆フィルタ処理をかけるように制御する。【選択図】　図１
2004108792.tif 000002
【特許請求の範囲】
【請求項１】
低周波音測定器により測定された測定波形データを記録しておく記録手段と、
低周波音測定器に内蔵されているハイパスフィルタと逆の特性を有するデジタル逆フィルタと、
前記記録手段から前記測定波形データを読み出し、この測定波形データに対して前記デジタル逆フィルタを適用し、原波形を復元するように、前記記録手段及び前記デジタル逆フィルタを制御する制御手段と、
を備えたことを特徴とする低周波音測定器用の波形復元装置。
【請求項２】
前記デジタル逆フィルタは、ＦＩＲ型逆フィルタであることを特徴とする請求項１記載の低周波音測定器用の波形復元装置。
【請求項３】
低周波音測定器への入力に相当する入力データを記録しておく第１の記録手段と、
低周波音測定器からの出力データを記録しておく第２の記録手段と、
前記第１の記録手段から出力される前記入力データと前記第２の記録手段から出力される前記出力データとから、前記低周波音測定器に内蔵されているハイパスフィルタの周波数応答関数を推定するハイパスフィルタ推定手段と、
前記推定されたハイパスフィルタの周波数特性と逆の特性を有するデジタル逆フィルタを設計するデジタル逆フィルタ設計手段と、
を備えたことを特徴とするデジタル逆フィルタ設計装置。
【請求項４】
前記第１の記録手段には、圧力変動測定時に前記低周波音測定器と並べて設置された差圧計の測定データが記録されていることを特徴とする請求項３記載のデジタル逆フィルタ設計装置
【請求項５】
前記第１の記録手段には、流体力学モデルからモデル計算により求められるモデル計算波形データが記録されていることを特徴とする請求項３記載のデジタル逆フィルタ設計装置。
【請求項６】
前記ハイパスフィルタ推定手段は、前記入力データ及び前記出力データを不規則振動波形とみなしてハイパスフィルタの推定を行うことを特徴とする請求項３乃至請求項５の何れか１項に記載のデジタル逆フィルタ設計装置。
【請求項７】
前記不規則振動波形とみなされた入力データ及び出力データは、規則振動波形とみなされる波形を少なくとも３つ含んでいることを特徴とする請求項６記載のデジタル逆フィルタ設計装置。
【請求項８】
前記デジタル逆フィルタ設計手段は、周波数サンプリング法を用いてＦＩＲ型逆フィルタを設計することを特徴とする請求項３乃至請求項７の何れか１項に記載のデジタル逆フィルタ設計装置。
【請求項９】
ハイパスフィルタを内蔵した低周波音測定器の測定波形データから波形のひずみを除去して原波形を復元する波形復元方法であって、
前記ハイパスフィルタと逆の特性を有するデジタル逆フィルタにより前記測定波形データをフィルタリングすることを特徴とする波形復元方法。
【請求項１０】
低周波音測定器に内蔵されているハイパスフィルタと逆の特性を有するデジタル逆フィルタを設計するデジタル逆フィルタ設計方法であって、
前記低周波音測定器への入力に相当する入力データと前記低周波音測定器からの出力データとからハイパスフィルタの周波数応答関数を推定する推定工程と、前記推定工程により推定されたハイパスフィルタの周波数特性と逆の特性を有するデジタル逆フィルタを設計する設計工程と、
を備えたことを特徴とするデジタル逆フィルタ設計方法。
【請求項１１】
低周波音測定器に内蔵されているハイパスフィルタと逆の特性を有するデジタル逆フィルタを設計するデジタル逆フィルタ設計方法であって、
低周波音測定器への入力に相当する入力データを記録する第１の記録工程と、
低周波音測定器からの出力データを記録する第２の記録工程と、
前記入力データと前記出力データとからハイパスフィルタの周波数応答関数を推定する推定工程と、
前記推定工程により推定されたハイパスフィルタの周波数特性と逆の特性を有するデジタル逆フィルタを設計する設計工程と、
を備えたことを特徴とするデジタル逆フィルタ設計方法。
【請求項１２】
前記第１の記録工程は、圧力変動測定時に前記低周波数音測定器と並べて配置された差圧計により測定した測定データを記録する工程であることを特徴とする請求項１１記載のデジタル逆フィルタ設計方法。
【請求項１３】
前記第１の記録工程は、流体力学モデルからモデル計算により求められるモデル計算波形データを記録する工程であることを特徴とする請求項１１記載のデジタル逆フィルタ設計方法。
【請求項１４】
前記推定工程は、前記入力データ及び出力データを不規則振動波形とみなしてハイパスフィルタの推定を行うことを特徴とする請求項１０乃至請求項１３の何れか１項に記載のデジタル逆フィルタ設計方法。
【請求項１５】
前記不規則波形とみなされた入力データ及び出力データは、規則振動波形とみなされる波形を少なくとも３つ含んでいることを特徴とする請求項１４記載のデジタル逆フィルタ設計方法。
【請求項１６】
前記設計工程は、周波数サンプリング法を用いてＦＩＲ型逆フィルタを設計する工程であることを特徴とする請求項１０乃至請求項１５の何れか１項に記載のデジタル逆フィルタ設計方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、低周波音測定器において波形を復元する方法に関し、詳しくは、２０Ｈｚ以下の超低周波音を含む低周波音を測定する低周波音測定器において波形を復元する方法に関する。
【０００２】
【従来の技術】
まず、測定対象である２０Ｈｚ以下の超低周波音を含む低周波音について簡単に説明する（詳しくは、非特許文献１参照）。
【０００３】
【非特許文献１】
中野有朋著，「超低周波音」，第１版，株式会社技術書院，２００２年８月２０日
【０００４】
この低周波音の発生源としては、自然的なものと人工的なものに分けることができる。自然的なものは、火山の噴火による爆発、地震による地表面の振動、海の波や津波による水面の振動、防波堤やテトラポットにあたる波、地中噴火による水中爆発及び台風や竜巻による圧力変化等によって発生する。
【０００５】
また、人工的なものは、大砲による砲撃、超音速航空機による圧力変化（ソニックブーム）、核爆発、橋梁の振動、ダムからの放水時の水流変動、新幹線のトンネル出口での圧力変化、高速列車の通過による圧力変化及び発破による衝撃等によって発生する。
【０００６】
このような低周波音は、建物の窓や戸などの建具をがたがた振動させて、騒音の問題を発生させ、また、直接人体に対して、頭痛、不眠、イライラ、耳鳴り等の悪影響をあたえるといわれている。
【０００７】
ここで、新幹線等の高速列車の通過による圧力変動について、図１５及び図１６を参照して、より具体的に説明する。
【０００８】
図１５は、高速で走行する列車周りに発生する圧力場を示す概念図である。このように、列車の先頭部では、進行方向前側にプラスの圧力場である正圧場が発生し、進行方向後側にマイナスの圧力場である負圧場が発生する。また、列車の後尾部では、進行方向前側に負圧場が発生し、進行方向後側に正圧場が発生する。更に、パンタカバー部においても、正圧場と負圧場とが発生する。
【０００９】
これらの圧力場は、定常的なものであり、音波等のようにそれ自身が変動して伝播するものではない。ところが、圧力場は列車の走行とともに移動するため、地上で観察すると低周波数の圧力変動として観測される。
【００１０】
図１６は、この列車通過時の圧力変動を示す概念図である。図１６に示すように、先頭部が通過する時に正圧・負圧、後尾部が通過する時に負圧・正圧、パンダカバー部が通過する時に正圧・負圧、の圧力変動が生じる。なお、後尾部の圧力変動が、先頭部のそれよりも小さくなっているのは、はく離の影響によるものである。
【００１１】
従来、このような圧力変動を測定する測定器として、低周波音圧レベル計や差圧計が用いられている。
【００１２】
【発明が解決しようとする課題】
図１７に、列車通過時の同一の圧力変動を、低周波音圧レベル計及び差圧計により測定した波形を示す。低周波音圧レベル計で測定した出力波形を点線で示し、差圧計で測定した出力波形を実線で示している。
【００１３】
図１７に示すように、低周波音圧レベル計の測定波形は、差圧計の測定波形と比較して、先頭部において正圧が過小・負圧が過大となり、後尾部において負圧が過小・正圧が過大となっている。
【００１４】
このように波形のひずみ現象が生じるのは、低周波音圧レベル計に内蔵されているハイパスフィルタの影響であると考えられる。
【００１５】
一方、差圧計の周波数特性は平坦であるため、差圧計を用いると波形のひずみのない測定波形が得られるが、気圧変動により０点がシフトしてしまう。このため、現地測定の際に、０点調整をする必要が生じる。
【００１６】
このような気圧変動は、差圧計に太陽光があたって温度が少し変化する等、小さな環境変化によって生じるため、頻繁に０点調整をする必要が生じ、非常に使い勝手を悪くしている。現場での測定にとって、測定器の操作性は重要であり、使い勝手が悪いと手間ばかりかかり、操作に集中せざるをえなくなり、その結果、測定意欲が減衰して測定精度へ悪影響を及ぼしてしまう。また、調整に時間がかかり、本来の測定時間が十分に確保できないといった問題が生じてしまう。
【００１７】
また、測定波形を処理する際にも０点シフトのための補正が必要になり、データ処理に時間と手間がかかってしまう。
【００１８】
以上より、操作性を考慮すれば、低周波音の測定器としては、低周波音圧レベル計を使用することが好ましい。しかし、上述した波形のひずみの問題を解決する必要がある。
【００１９】
本発明は、上記課題に鑑みてなされたものであり、操作性が良好で且つ波形のひずみのほとんどない低周波音測定を実現するために、低周波音測定器のフィルタにより生じた波形のひずみを取り除いて元の原波形を復元する波形復元方法及び波形復元装置を提供することを目的とする。
【００２０】
【課題を解決するための手段】
上記課題を解決するために、本発明に係る波形復元装置は、低周波音測定器により測定された測定波形データを記録しておく記録手段と、低周波音測定器に内蔵されているハイパスフィルタと逆の特性を有するデジタル逆フィルタと、記録手段から前記測定波形データを読み出し、この測定波形データに対して前記デジタル逆フィルタを適用し、原波形を復元するように、記録手段及びデジタル逆フィルタを制御する制御手段と、を備えたことを特徴とする。
【００２１】
また、上記課題を解決するために、本発明に係るデジタル逆フィルタ設計装置は、低周波音測定器への入力データを記録しておく第１の記録手段と、低周波音測定器からの出力データを記録しておく第２の記録手段と、第１の記録手段から出力される前記入力データと第２の記録手段から出力される出力データとから、低周波音測定器に内蔵されているハイパスフィルタの周波数応答関数を推定するハイパスフィルタ推定手段と、推定されたハイパスフィルタの周波数特性と逆の特性を有するデジタル逆フィルタを設計するデジタル逆フィルタ設計手段と、を備えたことを特徴とする。
【００２２】
また、上記課題を解決するために、本発明に係る波形復元方法は、ハイパスフィルタを内蔵した低周波音測定器の測定波形データから波形のひずみを除去して原波形を復元する波形復元方法であって、ハイパスフィルタと逆の特性を有するデジタル逆フィルタにより測定波形データをフィルタリングすることを特徴とする。
【００２３】
また、上記課題を解決するために、本発明に係るデジタル逆フィルタ設計方法は、低周波音測定器に内蔵されているハイパスフィルタと逆の特性を有するデジタル逆フィルタを設計するデジタル逆フィルタ設計方法であって、低周波音測定器への入力データと低周波音測定器からの出力データとからハイパスフィルタの周波数応答関数を推定する推定工程と、推定工程により推定されたハイパスフィルタの周波数特性と逆の特性を有するデジタル逆フィルタを設計する設計工程と、を備えたことを特徴とする。
【００２４】
【発明の実施の形態】
［第１の実施の形態］
以下、図面を参照しながら本発明の第１の実施の形態について説明する。
図１は、本実施の形態に係る波形復元装置１の概略構成を示すブロックダイヤグラムである。図１に示すように、波形復元装置１は、データ記録装置１１、デジタル逆フィルタ１２及び制御部１３から構成される。
【００２５】
データ記録装置１１は、低周波音圧レベル計３（リオン株式会社製の超低周波音圧レベル計ＸＮ−１２Ａ）により測定した測定データを記録している。デジタル逆フィルタ１２は、低周波音圧レベル計３に内蔵されているハイパスフィルタ３３と逆の周波数特性を有する逆フィルタであり、後述するデジタル逆フィルタ設計装置２により設計される。制御部１３は、データ記録装置１１及びデジタル逆フィルタ１２を制御し、データ記録装置１１からの出力にデジタル逆フィルタ処理をかけるように制御する。
【００２６】
図２は、本実施の形態に係るデジタル逆フィルタ設計装置２の概略構成を示すブロックダイヤグラムである。図２に示すように、デジタル逆フィルタ設計装置２は、差圧計データ記録装置２１、低周波音圧レベル計データ記録装置２２、ハイパスフィルタ推定部２３及びデジタル逆フィルタ設計部２４から構成される。
【００２７】
差圧計データ記録装置２１は、差圧計（株式会社ＳＳＫ製のＤＰ８−２０）により測定した測定データを記録している。低周波音圧レベル計データ記録装置２２は、低周波音圧レベル計３により測定した測定データを記録している。ハイパスフィルタ推定部２３は、後述する手法により、低周波音圧レベル計３に内蔵されているハイパスフィルタ３３の周波数特性を推定する。デジタル逆フィルタ設計部２４は、ハイパスフィルタ推定部２３により推定された周波数特性に基づき、後述する手法によりデジタル逆フィルタ１２を設計する。
【００２８】
図３は、本実施の形態に係る低周波音圧レベル計３の概略構成を示すブロックダイヤグラムである。図３に示すように、低周波音圧レベル計３は、入力側から、アッテネータ（ＡＴＴ）３１、アンプ（ＡＭＰ）３２、ハイパスフィルタ（ＨＰＦ）３３、ローパスフィルタ（ＬＰＦ）３４、アンプ（ＡＭＰ）３５から構成されている。また、カタログ値から、ハイパスフィルタは、０．２Ｈｚ以下の成分をカットし、ローパスフィルタは１０００Ｈｚ以上の成分をカットすることが分かっている。ただし、このハイパスフィルタ３３の具体的な周波数特性までは分かっていない。
【００２９】
ここで、本実施の形態の測定対象である圧力変動は、２０Ｈｚ以下に主な周波数成分がある低周波音であるので、ローパスフィルタ３４による波形のひずみは考慮せず、ハイパスフィルタ３３による波形のひずみのみを考慮すれば良い。
【００３０】
以下、本実施の形態に係る波形の復元について、（１）ハイパスフィルタ３３の周波数特性の推定、（２）デジタル逆フィルタ１２の設計、（３）波形の復元、という順番で説明する。
【００３１】
（１）ハイパスフィルタ３３の周波数特性の推定
まず、デジタル逆フィルタ設計装置２のハイパスフィルタ推定部２３におけるハイパスフィルタ３３の周波数特性の推定について説明する。
【００３２】
本実施の形態に係るハイパスフィルタ推定部２３では、差圧計データ記録装置２１からの出力を入力データ、低周波音圧レベル計データ記録装置２２からの出力を出力データとして、ハイパスフィルタ３３の周波数特性を推定する。なお、上記入力データは、実際の圧力変動測定時の低周波音圧レベル計３への入力に相当するデータである。
【００３３】
ハイパスフィルタ推定部２３においては、入力データのパワースペクトルをＰｘ（ｆ）、入出力データのクロススペクトルをＰｘｙ（ｆ）として、１つの線形システムとしての低周波音圧レベル計３の周波数応答関数を、Ｐｘｙ（ｆ）／Ｐｘ（ｆ）により求める。そして、後述するように、求めた周波数応答関数と、既知の各次数のバタワースフィルタの周波数応答を比較して、ハイパスフィルタ３３のフィルタ次数Ｌを推定する。
【００３４】
差圧計データ及び低周波音圧レベル計データは、同一の測定条件で同一の圧力変動を測定した測定データを用いた。具体的な測定対象は、高速列車通過時の圧力変動であり、低周波音圧レベル計と差圧計を並べて同時に測定した。
【００３５】
図４は、低周波音圧レベル計及び差圧計の測定位置を示す図である。図４に示すように、上り線の中心から６ｍ、８ｍ、１０ｍ、１２．５ｍ、２０ｍ、３０ｍの位置に低周波音圧レベル計及び差圧計を置いて測定を行った。
【００３６】
一般に、アナログ回路のフィルタとしては、バタワースフィルタ、チェビシェフフィルタ、楕円フィルタがよく知られているが、この中でも、バタワースフィルタが基本的なフィルタとして最もよく使用されている。そこで、本実施の形態に係る低周波音圧レベル計３のハイパスフィルタ３３として、バタワースフィルタが使用されていると仮定した。
【００３７】
ここで、バタワースフィルタを設計する場合、カットオフ周波数とフィルタ次数の２つのパラメータを与えれば良い。本実施の形態では、ハイパスフィルタ３３のカットオフ周波数は０．２Ｈｚであることが分かっているので、残りのフィルタ次数Ｌを推定すれば、ハイパスフィルタ３３の周波数特性を推定することができる。
【００３８】
図５は、６ｍの地点における差圧計による測定結果を示す図であり、図６は、６ｍの地点における低周波音圧レベル計３による測定結果を示す図である。図５の測定波形において丸印をつけたところは、測定時に強制的に０点調整したものである。
【００３９】
図５及び図６に示すように、本実施の形態に係る測定波形には、図１６に示す１列車分の波形が５つ含まれている。また、測定の際には、列車が通過している間だけデータレコーダを回しデータを記録しており、列車が通過していない間のデータは記録されていない。そして、差圧計データ記録装置２１には、図５に示す測定データが記録され、低周波音圧レベル計データ記録装置２２には、図６に示す測定データが記録されている。
【００４０】
ハイパスフィルタ推定部２３は、規則振動波形とみなされる１列車分の測定波形ごとに周波数特性の推定を行うのではなく、本実施の形態では５つの列車分の波形を、まとめて一つの不規則振動波形とみなして処理を行った。このように複数列車分の波形をまとめて処理することで、一つ一つの波形を個別に処理する場合に比べて、一回の処理で済むので、処理が簡単である。また、測定波形データのファイル数も少なくて済み、管理も容易になる。
【００４１】
また、一つの波形から推定される周波数特性は、波形が大きく変動する鋭い特性になりがちであるが、本実施の形態のように、複数の列車分の波形をまとめて処理することで、平均値のように滑らかな周波数特性を得ることができる。なお、波形の個数については、比較的滑らかな周波数特性を得るため、また、処理の簡素化のためにも、３個以上であることが望ましい。
【００４２】
なお、他の地点（８ｍ、１０ｍ、１２．５ｍ、２０ｍ、３０ｍ）についても、同様に測定・処理を行ったが、図５及び図６に示す６ｍの地点での測定波形とほとんど同じであったため、詳細な説明を省略し、６ｍの地点での測定波形についてのみ代表して説明する。
【００４３】
ハイパスフィルタ２３の周波数応答関数を求めるのは、上述したように、入力データのパワースペクトルをＰｘ（ｆ）、入出力のクロススペクトルをＰｘｙ（ｆ）とすると、Ｐｘｙ（ｆ）／Ｐｘ（ｆ）により求めることができる。
【００４４】
図７は、ハイパスフィルタ推定部２３において６ｍの地点での測定データから求めた周波数応答関数と、既知の各次数のバタワースフィルタ（１次乃至３次）の周波数特性との比較を示した図である。図７（ａ）が振幅特性の図であり、図７（ｂ）が位相特性の図である。
【００４５】
ここで、波形の復元を行う場合、振幅特性よりも位相特性が重要である。そこで、図７（ｂ）の位相特性に注目すると、１次のバタワースフィルタの位相特性が実測データから得られた位相特性に最も近い。よって、低周波音圧レベル計に用いられているハイパスフィルタは、１次のバタワースフィルタと推定される。
【００４６】
なお、この次数を推定する処理は、ハイパスフィルタ推定部２３において、最小二乗法等を利用したソフトウェア的な処理により求めても良いし、図７に示す比較図をディスプレイ等に表示して、オペレータが判断しても良い。
【００４７】
（２）デジタル逆フィルタ１２の設計
続いて、デジタル逆フィルタ設計部２４において、推定されたハイパスフィルタ３３の逆フィルタの設計を行う。本実施の形態では、周波数サンプリング法を用いたＦＩＲ型逆フィルタとして設計を行った（有本卓著，「信号・画像のディジタル処理」，産業図書，１９８０年参照）。
【００４８】
まず、周波数サンプリング法を用いてバタワースフィルタを設計するフィルタ設計法について簡単に説明する。
▲１▼入力パラメータ
入力パラメータとして、サンプリング間隔Δｔ（ｓ）、通過域周波数ｆｐ（Ｈｚ）、通過域周波数における減衰量Ａｐ（ｄＢ）、阻止域周波数ｆｒ（Ｈｚ）、阻止域周波数における減衰量Ａｒ（ｄＢ）、インパルス応答総数Ｎ（Ｎは奇数）を与える。
【００４９】
▲２▼角周波数への変換
通過域角周波数ωｐ＝２πｆｐ、阻止域角周波数ωｒ＝２πｆｒに変換する。
【００５０】
▲３▼バタワースフィルタの次数Ｍの決定
バタワースフィルタの次数Ｍは次式１で計算される。
【数１】
2004108792.tif 000003
【００５１】
▲４▼カットオフ周波数ωｃの決定
通過域から決まるカットオフ周波数ωｃｐ、及び阻止域から決まるカットオフ周波数ωｃｒは、次式２で表される。
【数２】
2004108792.tif 000004
このように、カットオフ周波数は２種類得られるが、本実施の形態では、カットオフ周波数ωｃを次式３で与えるものとする。
【数３】
2004108792.tif 000005
【００５２】
▲５▼ハイパスフィルタに変換
アナログプロトタイプのＭ次バタワースローパスフィルタＨＬ（ｓ）は、次式４で与えられる。
【数４】
2004108792.tif 000006
ここで、ＨＬ（ｓ）のｓｉ（ｉ＝１，２，…，Ｍ）を、次式５のようにすると、
【数５】
2004108792.tif 000007
【００５３】
周波数応答は、次式６のようになる。
【数６】
2004108792.tif 000008
上式で、ωをωｃ／ωに変換すると、Ｍ次バタワースハイパスフィルタＨＨ（ｊω）が、次式７のように得られる。
【数７】
2004108792.tif 000009
【００５４】
▲６▼周波数応答を単位円上でＮ個にサンプリング
角周波数ωを単位円上でＮ等分して、その分点を、ωｋ＝２πｋ／ＮΔｔ（ｋ＝０，１，…，Ｎ−１）とすると、振幅特性と位相特性は、次式８で表される。
【数８】
2004108792.tif 000010
ただし、実際の計算は、次のステップ▲７▼があるため、ｋ＝０，１，…，（Ｎ−１）／２を計算すれば良い。
【００５５】
▲７▼インパルス応答を実数にするための条件
インパルス応答を実数にするための次式９の条件を用いて、ｋ＝（Ｎ＋１）／２，…，Ｎ−１を計算する。
【数９】
2004108792.tif 000011
【００５６】
▲８▼インパルス応答の計算
インパルス応答は、次式１０の逆離散フーリエ変換を用いて計算する（ｎ＝０，１，２，…，Ｎ−１）。
【数１０】
2004108792.tif 000012
【００５７】
▲９▼フィルタリング
次式１１のたたみ込みを行う。
【数１１】
2004108792.tif 000013
ただし、実際の計算は高速フーリエ変換（ＦＦＴ）を用いて、高速化を図っている。
【００５８】
続いて、上述した周波数サンプリング法を用いて、デジタル逆フィルタ設計部２４において、ＦＩＲ型逆フィルタを設計する手法について説明する。
逆フィルタは、順フィルタ（バタワースハイパスフィルタ３３）の周波数特性Ｈ（ω）とは逆の周波数特性Ｈｉ（ω）＝１／Ｈ（ω）を持つフィルタである。したがって、逆フィルタの入力データである測定波形をｙ（ｔ）、出力データをｘ（ｔ−α）（α：任意の定数）とすると、次式１２で定義される。
【数１２】
2004108792.tif 000014
【００５９】
ここで、Ｈ（ω）、Ｈｉ（ω）を次式１３のように定義する。
【数１３】
2004108792.tif 000015
そうすると、逆フィルタの振幅と位相について、それぞれ次式１４で表される。
【数１４】
2004108792.tif 000016
【００６０】
このように、逆フィルタの振幅特性と位相特性は、順フィルタの振幅特性と位相特性から得られる。よって、逆フィルタの設計は、上述した周波数サンプリング法によるフィルタ設計法に、以下の作業を加えれば良い。
【００６１】
なお、バタワースハイパスフィルタの周波数特性Ｈ（ω）において、ω→０の時、Ｈ（ω）→０となるので、原理的に完全な復元は不可能である。しかし、列車通過時の圧力変動波形の主な周波数成分は１Ｈｚ付近にあり、ハイパスフィルタ３３のカットオフ周波数が０．２Ｈｚであることを考えれば、逆フィルタにより復元される波形には、必要とする周波数帯域の物理的情報は十分に入っていると考えられる。
【００６２】
また、上述したように、本実施の形態に係るバタワースハイパスフィルタ３３は１次であると推定されているが、一般性を持たせるために、以下の説明は、バタワースハイパスフィルタの次数をＭ次として説明する。
まず、アナログプロトタイプのＭ次バタワースローパスフィルタの伝達関数ＨＬ（ｓ）は、上述したように、一般に式４で表される。
【００６３】
このＨＬ（ｓ）の極ｓｉ（ｉ＝１，２，…，Ｍ）を、ｓｉ＝σｉ＋ｊωｉとすると、周波数特性｜Ｈｉ（ωｋ）｜、φｉ（ωｋ）は、次式１５及び１６により表される。ただし、Δｔはサンプリング間隔、ωｃはカットオフ角周波数である。
【数１５】
2004108792.tif 000017
【数１６】
2004108792.tif 000018
【００６４】
ただし、任意定数αは次式１７のように定義する。
【数１７】
2004108792.tif 000019
ここで、バタワースフィルタは非線形位相特性のため、必ずしも任意定数αを式１７のように選ぶ必要はない。実際、α＝０とおいても、結果的には同じであった。
【００６５】
しかし、線形位相特性のフィルタ設計においても、任意定数αは本来任意であり、α＝０としても良い。それを、式１７のようにおくことで完全線形位相特性が得られ、α＝０の場合よりさらに位相特性が改善される。よって、このことを考慮して、本実施の形態においても任意定数αを式１７のように定義した。
【００６６】
ところで、逆フィルタを用いた波形の復元では、Ｓ／Ｎ比の低い部分（本実施の形態では、０．２Ｈｚ以下の部分）を復元しても、ノイズばかりを復元して物理的に意味のあるデータは得られない。
【００６７】
よって、本実施の形態に係るデジタル逆フィルタ設計部２４では、周波数領域（０≦ｆ≦１／２Δｔ）を、復元領域、遮断領域、遷移領域の３つに分けて設計を行う。このように、遷移帯域を設けたのは、フィルタ設計が数学的には一種の近似であり、したがって急激な遮断は数学的に不連続をもたらし、既知のＧｉｂｂｓ現象により、カットオフ周数数付近で振幅利得にかえって誤差をまねくからである。
【００６８】
図８は、上記３つの領域に分けて設計される逆フィルタの振幅特性｜Ｈｉ（ｆ）｜の概略を示す図である。横軸は周波数ｆ、縦軸は振幅特性｜Ｈｉ（ｆ）｜を表しており、ｇ１は、遷移帯域の下限周波数、ｇ２は、遷移帯域の上限周波数である。また、遷移帯域は余弦関数で接続している。
【００６９】
このように逆フィルタを３つの領域から設計することで、Ｓ／Ｎ比の低い部分が復元されることを防止できると共に、遷移領域により測定誤差を緩和させることができる。
【００７０】
（設計例）
表１は、デジタル逆フィルタ設計部２４において、実際に逆フィルタの設計を行った際の各入力パラメータの値を示す。次数Ｍ＝１の場合だけでなく、比較のため、次数Ｍ＝２，３の場合についても設計値を示してある。なお、この入力パラメータは、図示しない入力手段から、オペレータによりデジタル逆フィルタ設計部に入力されるようにすると良い。
【表１】
2004108792.tif 000020
【００７１】
図９は、上記パラメータで設計された逆フィルタ（次数Ｍ＝１）の周波数特性を示す図である。図９（ａ）は、振幅特性を示す図であり、図９（ｂ）は位相特性を示す図である。比較のために、順フィルタ（１次バタワースハイパスフィルタ）の周波数特性も合わせて示している。図９を見ると、逆フィルタが適切に設計されていることが分かる。２次、３次の逆フィルタについては、図示を省略するが、１次の場合と同様に、適切に逆フィルタが設計されていた。
【００７２】
（３）波形の復元
このように設計されたデジタル逆フィルタを波形復元装置１のデジタル逆フィルタ１２として適用すると、低周波音圧レベル計３による測定波形から波形のひずみを除去して原波形を復元することができる。なお、本実施の形態では、デジタル逆フィルタ１２の作用をコンピュータによりソフトウェア的に処理しているが、ハードウェア的にフィルタ回路を構成してフィルタリングしても良い。
【００７３】
図１０は、本実施の形態に係る波形復元装置１により復元された波形を示す図である。なお、図１０には、差圧計による測定波形（原波形）も一緒に示してある。また、上述したように、差圧計による測定波形には０点シフトが見られるが、ここでは０点シフトを除去する補正を行った後の測定波形を示している。
【００７４】
差圧計は周波数特性が平坦であると考えられるので、復元された波形が差圧計波形に近ければ、復元が適切に行われていると判断できるが、図１０を見ると、二つの波形はほぼ一致している。
以上、詳細に説明したように、本実施の形態によれば、低周波音圧レベル計３に用いられているハイパスフィルタ３３が、１次バタワースハイパスフィルタであることが推定できた。また、このハイパスフィルタ３３による波形のひずみを除去するため、ＦＩＲ型のデジタル逆フィルタ１２を設計し、この逆フィルタを適用することにより、低周波音圧レベル計１の測定波形のひずみを取り除くことができた。
【００７５】
［第２の実施の形態］
以下、図面を参照しながら本発明の第２の実施の形態について説明する。
本実施の形態と第１の実施の形態とは、（１）ハイパスフィルタ３３の周波数特性の推定、（２）デジタル逆フィルタ１２の設計、の手法が大きく異なる。よって、第１の実施の形態と異なる事項のみを説明し、共通する事項については説明を省略する。
【００７６】
（１）ハイパスフィルタ３３の周波数特性の推定
第１の実施の形態では、デジタル逆フィルタ設計装置２の入力データとして、差圧計により測定したデータを用いている。差圧計の測定データは、差圧計と低周波音圧レベル計３とを並べて測定して得られたものであるので、列車速度、測定位置、天候の影響等の測定条件が、差圧計と低周波音圧レベル計とで同一である。よって、周波数特性を推定するための入力データとしては、この点で好ましい。
【００７７】
しかし、現場での測定によるノイズがどうしても入ってくると考えられる。
このため、第２の実施の形態では、ハイパスフィルタ推定部２３への入力データとして、簡単な流体力学モデルから計算される波形データを用いた。すなわち、図示は省略するが、第２の実施の形態に係るデジタル逆フィルタ設計装置２は、図２の差圧計データ記録装置の変わりに、後述するモデル計算波形を生成するモデル波形計算装置を備えている。
【００７８】
図１１は、この流体力学モデルを示す概念図である。ｘ軸は、列車の進行方向に平行な軸、ｙ軸は、列車の進行方向に垂直であって地面に平行な軸、ｚ軸は、列車の進行方向及び地面に垂直な軸である。図１１（ａ）は、ｘ−ｙ平面図であり、下側のｘ軸上を列車が通過するようになっている。また、上側のｘ軸は、下側のｘ軸から距離ｒ離れた位置にあり、このｘ軸上を０点として距離ｒ離れた位置での圧力変化を図示している。すなわち、ｙ軸は、列車の通過位置からの距離を表すと共に、列車の通過位置から距離ｒ離れた位置での圧力も表わしている。
【００７９】
図１１（ｂ）は、ｙ−ｚ平面図であり、ｘ軸は地平面上に位置し、ｚ軸は地面からの距離を表わしている。また、Ｐ（ｘ，ｙ，ｚ）は測定を行う観測点であり、地面からｚ１［ｍ］の高さに位置している。ｍはわき出しの中心部であり、ｈ［ｍ］は列車底面の地面からの高さ、ｂ［ｍ］は列車の断面を円形状をみなした時の半径である（後述するＡ’を用いれば、Ａ’＝πｂ２の関係）。
【００８０】
Ｕを列車速度［ｍ／ｓ］、ｍ［ｍ３／ｓ］をわき出しの強さ、ｒ１、ｒ２［ｍ］をわき出し点と観測点Ｐ（ｘ，ｙ，ｚ）との距離とし、地表面の効果を鏡像法を用いて表わすと、流れを示す速度ポテンシャルφ［ｍ２／ｓ］は、次式１８で表わされる。
【数１８】
2004108792.tif 000021
【００８１】
また、ｒ１、ｒ２は次式１９で表わされる。
【数１９】
2004108792.tif 000022
また、わき出しの強さｍは、列車断面積をＡ’［ｍ２］とすると、次式２０で表わされる。
【数２０】
2004108792.tif 000023
【００８２】
式１９の勾配をとり、流速を求め、Ｂｅｒｎｏｕｌｌｉの定理を使って整理すると、圧力係数Ｃｐは、次式２１で表わされる。
【数２１】
2004108792.tif 000024
圧力は、圧力係数Ｃｐに動圧ρＵ２／２をかければ良い。ここで、ρ［ｋｇ／ｍ３］は密度である。よって、圧力ｐは、次式２２で表わされる。
【数２２】
2004108792.tif 000025
【００８３】
以上、ｍをわき出しとみなした場合について説明した。つまり列車の先頭部を表現したものである。
一方、列車の後尾部では、ｍを吸い込みとみなして、先頭部と同様に表現できる。ただし、後尾部では、はく離の影響により、先頭部よりも圧力変動が小さくなる。このため、所定の値として後尾部圧力低減定数γを定義し、後尾部の吸い込み強さをγｍとする。この後尾部圧力低減係数γは、低周波音圧レベル計のｐｅａｋ−ｐｅａｋ値と、本実施の形態に係るモデル計算波形のｐｅａｋ−ｐｅａｋ値が合うように、所定の値に決定する。
【００８４】
また、時間の原点を適当に決め、ｔ＝ｔｉの時、先頭部のわき出し又は後尾部の吸い込みがｘ＝０にあるとし、ｘ＝Ｕ（ｔ−ｔｉ）と変換することで、時間領域のグラフになる。
【００８５】
図１１（ａ）に示すように、圧力が０になる点がわき出し点になるので、実測データから通過時圧力が０になる時刻を読み取り、そこに各列車に対応するわき出しをおけば良い。
このようにしてモデル波形計算装置において求めたモデル計算波形を図１２に示す。Ａ系、Ｂ系、Ｃ系の三種類の列車からなる５つ分の列車の波形について示している。図１２（ａ）は、上述したモデル計算により求めた波形を示す図であり、図１２（ｂ）は、比較のために低周波音圧レベル計による測定波形を示している。
【００８６】
実際の計算は、ｙ＝６．０ｍ、８．０ｍ、１０．０ｍ、１２．５ｍ、２０．０ｍ、３０ｍの６つの地点でのモデル計算を行ったが、代表してｙ＝６．０ｍの地点についてのみ結果を示している。
【００８７】
このモデル計算におけるパラメータは以下の通りである。列車断面積：Ａ’＝１２．２ｍ２（Ａ系）、１２．９ｍ２（Ｂ系）、１１．２ｍ２（Ｃ系）、列車底面高さ：ｈ＝０．５ｍ、観測点高さ：ｚ１＝１．２ｍ、後尾部圧力低減係数：γ＝０．６６１。また、上り線と下り線の間隔は、図４に示す第１の実施の形態と同様に４．３ｍである。
【００８８】
また、図１２に示す５つの列車は、左側から、上りＢ系１６両編成（２２０ｋｍ／ｈ）、下りＡ系６両編成（２４１ｋｍ／ｈ）、上りＡ系６両編成（２１７ｋｍ／ｈ）、下りＢ系１６両編成（２１４ｋｍ／ｈ）、上りＣ系１６両編成（２６５ｋｍ／ｈ）である。
【００８９】
図１３は、入力データとしてこのモデル計算波形（モデル波形計算装置）、出力データとして低周波音圧レベル計３による測定波形（低周波音圧レベル計データ記録装置２２）とした時の周波数応答関数を示す図である。図１３（ａ）は、振幅特性を示し、図１３（ｂ）は、位相特性を示している。周波数応答関数を求める方法は第１の実施の形態と同様であるので、説明を省略する。また、第１の実施の形態と同様に、ハイパスフィルタ推定部２３において、低周波音圧レベル計３に用いられているハイパスフィルタ３３の周波数特性を推定するために、図１３には、既知の各次数のバタワースフィルタ（１次乃至３次）の周波数特性を示している。
【００９０】
図１３に示すように、本実施の形態で求めた周波数応答関数は、２Ｈｚより高周波数側で振幅特性及び位相特性の振動が大きくなっている。これは、モデル計算波形にはこの周波数範囲のデータが含まれていないためと思われる。
【００９１】
そして、本実施の形態で求めた周波数応答関数が、何次のバタワースフィルタの周波数特性に近いかについてみるため、位相特性に注目すると、１次バタワースフィルタの特性が最も近い。よって、第２の実施の形態においても、第１の実施の形態と同様に、ハイパスフィルタ推定部２３において、低周波音圧レベル計３に用いられているハイパスフィルタ３３は、１次のバタワースハイパスフィルタであると推定される。
【００９２】
なお、本実施の形態では、第１の実施の形態の差圧計データ記録装置に代えて、モデル波形計算装置を用いて周波数特性の推定を行っているが、これに限定されるものではない。例えば、デジタル逆フィルタ設計装置２が、差圧計データ記録装置２１及びモデル波形計算装置の双方を備えるように構成しても良い。そして、差圧計による測定波形を入力データとして周波数特性の推定を行った後に、さらにモデル計算波形を入力データとして周波数特性の推定を行っても良い。このように二通りの方法で推定を行うことで、推定結果の信頼性をより高めることができる。
【００９３】
（２）デジタル逆フィルタ１２の設計
第１の実施の形態では、デジタル逆フィルタ設計部２４において、周波数サンプリング法を用いたＦＩＲ型逆フィルタとして設計を行っているが、本実施の形態では、双一次変換を用いたＩＩＲ型逆フィルタとして設計を行った。
まず、双一次変換を用いた設計法について簡単に説明する。
【００９４】
▲１▼入力パラメータ
入力パラメータとして、サンプリング間隔Δｔ（ｓ）、通過域周波数ｆｐ（Ｈｚ）、通過域周波数における減衰量Ａｐ（ｄＢ）、阻止域周波数ｆｒ（Ｈｚ）、阻止域周波数における減衰量Ａｒ（ｄＢ）を与える。
【００９５】
▲２▼周波数の正規化
上で与えた通過域周波数ｆｐ、阻止域周波数ｆｒをサンプリング周波数ｆｓで正規化し、角周波数に変換する。すなわち、サンプリング周波数ｆｓ＝１／Δｔ、デジタルフィルタ通過域角周波数ωｄｐ＝２πｆｐ／ｆｓ、ωｄｒ＝２πｆｒ／ｆｓと変換する。
【００９６】
▲３▼バタワースフィルタの次数Ｌの決定
バタワースフィルタの次数Ｌは、次式２３で決定される。
【数２３】
2004108792.tif 000026
【００９７】
▲４▼カットオフ周波数ωｃの決定
通過域から決まるカットオフ周波数ωｃｐ、及び阻止域から決まるカットオフ周波数ωｃｒは、第１の実施の形態と同様に式２で表される。そして、本実施の形態でも同様に、カットオフ周波数ωｃを式３で与えるものとする。
【００９８】
▲５▼アナログプロトタイプをデジタルプロトタイプに変換
Ｍ次バタワースローパスフィルタのアナログプロトタイプＨＬａは、次式２４で表わされる。
【数２４】
2004108792.tif 000027
ただし、ｓｉ（ｉ＝１，２，…，Ｍ）は、ＨＬａ（ｓ）の安定なＭ個の極である。
続いて、式２４において、ｓを次式２５のように変換する。
【数２５】
2004108792.tif 000028
【００９９】
すると、次式２６に示すような、Ｍ次バタワースローパスフィルタのデジタルプロトタイプＨＬｄが得られる。
【数２６】
2004108792.tif 000029
【０１００】
▲６▼ハイパスフィルタに変換
デジタルプロトタイプのＭ次バタワースローパスフィルタを、ωｃをカットオフ周波数に持つハイパスフィルタに変換する。すなわち、上式２６のｚ−１を次式２７のように変換する。
【数２７】
2004108792.tif 000030
以上の変換を適用して、式２８のように整理し、各係数ａｉ、ｂｉを求める。
【数２８】
2004108792.tif 000031
【０１０１】
▲７▼フィルタリング
Ｍ次バタワースデジタルフィルタを次式２９でフィルタリングする。
【数２９】
2004108792.tif 000032
ただし、ｙ（ｎ−ｉ）＝０、ｎ−ｉ＜０である。
続いて、本実施の形態に係るＩＩＲ型逆フィルタを設計する手法について説明する。
【０１０２】
上述したように、順フィルタであるハイパスフィルタの周波数特性Ｈ（ｚ）は、式２８で表わされる。よって、逆フィルタの周波数特性Ｈｉ（ｚ）は、次式３０のようになる。
【数３０】
2004108792.tif 000033
【０１０３】
つまり、順フィルタの零点が、逆フィルタの極になる。フィルタの次数Ｍ＝６次まで逆フィルタの極を調べると、ｚ−１＝１（Ｍ重根）となっている。デジタルフィルタの安定条件は、「極がすべてｚ平面状の単位円の外部にある」ことである。したがって、上式３０のままでは安定しない。
【０１０４】
そこで、逆フィルタの安定性を確保するため、極を単位円の外に移動させる。すなわち、本実施の形態では、逆フィルタの分母を、ｃ（ｚ−１−ｒ・１）Ｍ、１＜ｒ、と変形し、これを展開して式３０の形に整理し、各係数を求める。このため、本実施の形態では、入力パラメータとしてｒが追加されることになる。
【０１０５】
（設計例）
表２は、本実施の形態に係る方法により、実際に逆フィルタの設計を行った際の各入力パラメータの値を示す。次数Ｍ＝１の場合だけでなく、比較のため、次数Ｍ＝２，３の場合についても設計値を示してある。
【表２】
2004108792.tif 000034
表３は、上記の入力パラメータを与えた時の出力結果を示す表である。
【表３】
2004108792.tif 000035
【０１０６】
（３）波形の復元
図１４は、本実施の形態に係る波形復元装置１により復元された波形を示す図である。また、図１４には、差圧計による測定波形も一緒に示してある。
図１０に示す第１の実施の形態に係るＦＩＲ型逆フィルタの場合と比較すると、ＦＩＲ型逆フィルタによる復元波形が安定しているのに対して、ＩＩＲ型逆フィルタによる復元波形には大きなうねりのような振動が含まれている。よって、逆フィルタとしてはＦＩＲ型のほうが適していると考えられる。
【０１０７】
一方、本実施の形態のモデル計算波形のように、ノイズを含まない波形の復元の場合には、ＩＩＲ型逆フィルタのほうが適していることがわかっている。
よって、低周波音圧レベル計３による実測波形を復元するといった、ノイズを含む測定データを復元する場合にはＦＩＲ型逆フィルタが適しているが、ノイズの少ない環境での測定データであれば、ＩＩＲ型フィルタも適用できると考えられる。
【０１０８】
以上、詳細に説明したように、本実施の形態によれば、第１の実施の形態と異なる構成により、第１の実施の形態における効果と同様の効果が得られる。
なお、上述したように、第１の実施の形態では、入力データを差圧計の測定データ、デジタル逆フィルタをＦＩＲ型逆フィルタとする組み合わせであり、また、第２の実施の形態では、入力データをモデル計算波形、デジタル逆フィルタをＩＩＲ型逆フィルタとする組み合わせとなっているが、本発明の実施の形態は、もちろんこの組み合わせに限定されるものではない。例えば、入力データを差圧計の測定データ、デジタル逆フィルタをＩＩＲ型逆フィルタとする組み合わせ、入力データをモデル計算波形、デジタル逆フィルタをＦＩＲ型逆フィルタとする組み合わせであっても良い。
【０１０９】
【発明の効果】
以上詳細に説明したように、本発明によれば、低周波音測定器のフィルタにより生じた波形のひずみを取り除いて元の原波形を復元する波形復元装置及び波形復元方法を提供することができ、操作性が良好で且つ波形のひずみのほとんどない低周波音測定を実現することができる。
【図面の簡単な説明】
【図１】第１の実施の形態に係る波形復元装置の概略構成を示すブロックダイヤグラムである。
【図２】第１の実施の形態に係るデジタル逆フィルタ設計装置の概略構成を示すブロックダイヤグラムである。
【図３】第１の実施の形態に係る低周波音圧レベル計の概略構成を示すブロックダイヤグラムである。
【図４】第１の実施の形態に係る差圧計と低周波音圧レベル計の測定位置を示す図である。
【図５】第１の実施の形態に係る６ｍの地点での差圧計による測定波形を示す図である。
【図６】第１の実施の形態に係る６ｍの地点での低周波音圧レベル計による測定波形を示す図である。
【図７】第１の実施の形態に係る差圧計の測定波形から求めた周波数特性を示す図である。
【図８】第１の実施の形態に係る逆フィルタの周波数特性の概略を示す図である。
【図９】第１の実施の形態に係る逆フィルタの周波数特性を示す図である。
【図１０】第１の実施の形態に係る波形復元装置により復元された波形を示す図である。
【図１１】第２の実施の形態に係る流体力学モデルを示す概念図である。
【図１２】第２の実施の形態に係るモデル計算により求めた波形を示す図である。
【図１３】第２の実施の形態に係るモデル計算波形から求めた周波数特性を示す図である。
【図１４】第２の実施の形態に係る波形復元装置により復元された波形を示す図である。
【図１５】高速で走行する列車周りに発生する圧力場を示す概念図である。
【図１６】列車通過時の圧力変動を示す概念図である。
【図１７】列車通過時の圧力変動を低周波音圧レベル計及び差圧計により測定した測定波形を示す図である。
【符号の説明】
１　　波形復元装置
２　　デジタル逆フィルタ設計装置
３　　低周波音圧レベル計
１１　データ記録装置
１２　デジタル逆フィルタ
１３　制御部
２１　差圧計データ記録装置
２２　低周波音圧レベル計データ記録装置
２３　ハイパスフィルタ推定部
２４　デジタル逆フィルタ設計部
３１　アッテネータ
３２　アンプ
３３　ハイパスフィルタ
３４　ローパスフィルタ
３５　アンプ
【図１】
2004108792.tif 000036
【図２】
2004108792.tif 000037
【図３】
2004108792.tif 000038
【図４】
2004108792.tif 000039
【図５】
2004108792.tif 000040
【図６】
2004108792.tif 000041
【図７】
2004108792.tif 000042
【図８】
2004108792.tif 000043
【図９】
2004108792.tif 000044
【図１０】
2004108792.tif 000045
【図１１】
2004108792.tif 000046
【図１２】
2004108792.tif 000047
【図１３】
2004108792.tif 000048
【図１４】
2004108792.tif 000049
【図１５】
2004108792.tif 000050
【図１６】
2004108792.tif 000051
【図１７】
2004108792.tif 000052