(57)【要約】  未知の信号に対して、入力信号を計測する。そしてその値から、直接的にその入力信号中に含まれる正弦波のパラメータ（周波数、振幅、位相）を推定する技術を提供する。  解析信号の支配方程式を立て、さらに任意の荷重関数を導入することによって、積分形式に置き換える。それを解くことによって、正弦波パラメータが推定される。その結果、従来の近似による推定ではなく、より厳密に正弦波のパラメータを得ることができるようになった。
2008136443.tif 000002
【特許請求の範囲】
【請求項１】
解析の対象である信号中の正弦波のパラメータを推定する方法において、
前記信号の微分を含む支配方程式を立てる工程と、
前記微分を含む支配方程式に任意の荷重関数を導入し、有限区間の積分形式に置き換える工程と、
前記積分形式に変換した支配方程式を解く工程と、
を含むことを特徴とする正弦波パラメータ推定方法。
【請求項２】
解析の対象である信号中の正弦波のパラメータを推定する方法において、
任意の荷重関数を導入した有限区間の積分形式で、前記信号の微分を含む支配方程式を構成する工程と、
前記積分形式の支配方程式を解く工程と、
を含むことを特徴とする正弦波パラメータ推定方法。
【請求項３】
請求項１又は２記載の正弦波パラメータ推定方法において、
前記荷重関数は、前記有限期間の整数分の１の周期の三角関数であることを特徴とする正弦波パラメータ推定方法。
【請求項４】
請求項１又は２記載の正弦波パラメータ推定方法において、
前記積分形式の支配方程式を解く工程は、
異なる３個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、
前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、
を含むことを特徴とする正弦波パラメータ推定方法。
【請求項５】
請求項１又は２記載の正弦波パラメータ推定方法において、
前記積分形式の支配方程式を解く工程は、
異なる２個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、
前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、
を含むことを特徴とする正弦波パラメータ推定方法。
【請求項６】
請求項１又は２記載の正弦波パラメータ推定方法において、
前記積分形式の支配方程式を解く工程は、
３個以上の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、
前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、
を含むことを特徴とする正弦波パラメータ推定方法。
【請求項７】
請求項１又は２記載の正弦波パラメータ推定方法において、
前記積分形式の支配方程式を解く工程は、
１個の荷重周波数と２個以上の異なる窓関数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、
前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、
を含むことを特徴とする正弦波パラメータ推定方法。
【請求項８】
請求項４〜７のいずれかに記載の正弦波パラメータ推定方法において、
前記パラメータを推定する工程は、
前記推定された周波数と、積分境界値項と、その推定された周波数の信号の振幅と位相とを正弦波のパラメータとして求める工程、
を含むことを特徴とする正弦波パラメータ推定方法。ここで、積分境界値項とは積分区間の信号やその微分の両端の値の差をいう。
【請求項９】
請求項４〜７のいずれかに記載の正弦波パラメータ推定方法において、
前記パラメータを推定する工程は、
前記推定された周波数に関して、位相と振幅とをパラメータとして含む正弦波モデルを構築し、この正弦波モデルを解析対象となる解析信号にフィッティングさせることによって、最もフィットする正弦波モデルの振幅と位相とを求める工程、
を含むことを特徴とする正弦波パラメータ推定方法。
【請求項１０】
請求項１〜９のいずれかに記載の正弦波パラメータ推定方法において、
前記積分形式に置き換える工程又は処理は、
前記積分の積分境界を、前記信号を標本化するタイミングの中間に設定することを特徴とする正弦波パラメータ推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響センシング、音響信号処理分野に関し、特に信号波形中に含まれる正弦波パラメータを推定する方法に関する。【０００２】
さらに本発明は、音響や光など、広く波動場のセンシングや信号処理、通信やデータ記録における復調や復号化の技術においては、正弦波パラメータの推定は極めて重要な技術事項であるので、それらの分野にも強く関連する。
【背景技術】
【０００３】
音声信号などの信号波形解析においては、その信号に含まれる正弦波のパラメータを推定することは基本的、且つ重要な事項である。【０００４】
したがって、従来から信号中の正弦波パラメータ推定技術として種々のものが提案されている。広く利用されている技術は、ＦＦＴ（高速フーリエ変換）を利用したものであり、特にその得られたＦＦＴスペクトルを２次補間して正弦波パラメータを推定する方法は簡単で精度も高いので広く利用されている。【０００５】
先行技術文献
正弦波パラメータ推定と関連性があると思われるいくつかの先行技術を挙げる。【０００６】
下記特許文献１には、正弦波の振幅と周波数とを推定する方法等が開示されている。ここに開示されている技術は、サンプリング信号に基づき、これらの信号間を補間すること等によって正弦波の推定を行うことを特徴とする。【０００７】
また、下記特許文献２には、未知の信号を推定するスペクトル分析装置が開示されている。ここに開示されているスペクトル分析装置によれば、入力信号の周波数や位相が推定できるとされている。【０００８】
また、下記特許文献３には、入力音声の正弦波を抽出し、その正弦波を利用して音声合成する技術が開示されている。【０００９】
また、下記特許文献４には、ノイズが音声信号に混入した場合、そのノイズの期間における正弦波を推定し、入力信号を補間することによってノイズキャンセルを行う装置が開示されている。【００１０】
また、下記非特許文献１においては、ＦＦＴの２次補間による正弦波パラメータの推定の基準に関する研究が開示されている。特に、その設計パラメータ、窓関数や、窓の長さ、零詰めを行う際の零詰めの量等について研究がなされている。【００１１】
【特許文献１】特開平０７−０８３９６４号公報
【特許文献２】特開２０００−２５８４７８号公報
【特許文献３】特開２００４−１０９８０９号公報
【特許文献２】特開２００６−１３５８０６号公報
【非特許文献１】安部, スミス,  ”ＦＦＴの２次補間に基づく正弦波パラメータ推定法の設計基準”,信学技報, vol.104 No.304, pp.7-12, Nov. 2004.
【発明の開示】
【発明が解決しようとする課題】
【００１２】
このように、未知の信号に対して正弦波のパラメータを調べる従来の手法は、補間等のテクニックを用いた近似であり、精度の向上にも一定の限界があると考えられる。【００１３】
そこで、近似ではなく、直接的に正弦波のパラメータ（振幅、位相）を推定できる手法が望ましいことは言うまでもない。【００１４】
理論的には、ナイキストの定理から、サンプリング周期が十分に高ければそのサンプリングした信号の値から（量子化誤差を除けば）元の正弦波信号は完全に再現できるはずであり、そのような手法が望まれている。【００１５】
本願発明は、このような背景に鑑みなされたものであり、補間等の近似ではなく、入力信号を計測した値から、直接的にその入力信号中に含まれる正弦波のパラメータを推定する技術を実現することを目的とする。
【課題を解決するための手段】
【００１６】
本願発明は、上記課題を解決するために、微分方程式の有限区間荷重積分に基づく正弦波パラメータの直接推定を行うことをその解決原理とするものである。【００１７】
また、本特許においては、正弦波のパラメータとは、周波数と振幅と位相を意味する。なお、正弦波のパラメータという場合、周波数と振幅と位相とに加えて直流バイアスを意味する場合もある。【００１８】
本特許においては、この直流バイアスを直接パラメータとして扱ってはいないが、周波数が直流という特別の場合の振幅として推定することができる。【００１９】
具体的には、上記課題を解決するため以下のような手段を採用している。【００２０】
（１）本発明は、上記課題を解決するために、解析の対象である信号中の正弦波のパラメータを推定する方法において、前記信号の微分を含む支配方程式を立てる工程と、前記微分を含む支配方程式に任意の荷重関数を導入し、有限区間の積分形式に置き換える工程と、前記積分形式に変換した方程式を解く工程と、を含むことを特徴とする正弦波パラメータ推定方法である。【００２１】
（２）本発明は、解析の対象である信号中の正弦波のパラメータを推定する方法において、任意の荷重関数を導入した有限区間の積分形式で、前記信号の微分を含む支配方程式を構成する工程と、前記積分形式の支配方程式を解く工程と、を含むことを特徴とする正弦波パラメータ推定方法。【００２２】
（３）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記荷重関数は、前記有限期間の整数分の１の周期の三角関数であることを特徴とする正弦波パラメータ推定方法である。【００２３】
整数分の１にとることによって、後述するように、未知数や計数の個数を近似無しで減じ、かつ、必要な個数は残すことが可能である。【００２４】
（４）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、異なる３個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、を含むことを特徴とする正弦波パラメータ推定方法である。【００２５】
（５）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、異なる２個の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、を含むことを特徴とする正弦波パラメータ推定方法である。【００２６】
（６）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、３個以上の荷重周波数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、を含むことを特徴とする正弦波パラメータ推定方法である。【００２７】
（７）また、本発明は、上記（１）又は（２）記載の正弦波パラメータ推定方法において、前記積分形式の支配方程式を解く工程は、１個の荷重周波数と２個以上の異なる窓関数の荷重関数を導入し、それらに関する連立方程式群を構成する工程と、前記連立方程式群を解くことによって前記信号中の正弦波の少なくとも周波数を含むパラメータを推定する工程と、を含むことを特徴とする正弦波パラメータ推定方法である。【００２８】
（８）また、本発明は、上記（４）〜（７）のいずれかに記載の正弦波パラメータ推定方法において、前記パラメータを推定する工程は、前記推定された周波数と、積分境界値項と、その推定された周波数の信号の振幅と位相とを正弦波のパラメータとして求める工程、を含むことを特徴とする正弦波パラメータ推定方法。ここで、積分境界値項とは積分区間の信号やその微分の両端の値の差をいう。【００２９】
（９）また、本発明は、上記（４）〜（７）のいずれかに記載の正弦波パラメータ推定方法において、前記パラメータを推定する工程は、前記推定された周波数に関して、位相と振幅とをパラメータとして含む正弦波モデルを構築し、この正弦波モデルを解析対象となる解析信号にフィッティングさせることによって、最もフィットする正弦波モデルの振幅と位相とを求める工程、を含むことを特徴とする正弦波パラメータ推定方法である。【００３０】
（１０）また、本発明は、上記（１）〜（９）のいずれかに記載の正弦波パラメータ推定方法において、前記積分形式に置き換える工程又は処理は、前記積分の積分境界を、前記信号を標本化するタイミングの中間に設定することを特徴とする正弦波パラメータ推定方法である。【００３１】
参考発明の手段
以下、本件発明と関連する参考発明の手段を説明する。【００３２】
上で述べた正弦波パラメータ推定方法は、以下のようなハードウェアを用いて実施することが好適である。【００３３】
（１１）本参考発明は、上記（１）〜（１０）のいずれかに記載の正弦波パラメータ推定方法の前記荷重積分を支援する装置において、前記解析信号に前記荷重を乗算するＮ個の乗算手段と、前記Ｎ個の乗算手段毎にそれぞれ設けられ、前記乗算手段の各出力信号をＮ周期遅延させるＮ個の遅延手段と、前記Ｎ個の遅延手段毎にそれぞれ設けられ、遅延後の出力信号と、遅延前の信号との差を求め、その差を累積値に累積するＮ個の加算手段と、前記加算手段毎にそれぞれ設けられ、前記加算手段の出力信号を累積するＮ個の累積手段と、を含むことを特徴とする正弦波パラメータ推定方法支援装置である。ここで、前記Ｎは正の整数である。【００３４】
（１２）さらに、上で述べた正弦波パラメータ推定方法を実施する為に、任意の荷重関数を導入した有限区間の積分形式で、前記信号の支配方程式を構成する手段と、前記方程式を解く手段と、を含むことを特徴とする正弦波パラメータ推定支援装置を用いることが好ましい。
【発明の効果】
【００３５】
以上述べたように、本発明によれば、解析信号中に含まれる正弦波のパラメータを直接推定することができる。
【図面の簡単な説明】
【００３６】
【図１】本実施例における反復法２の処理の流れのフローチャートである。
【図２】本実施例におけるシミュレーションした周波数の推定結果を表すグラフである。
【図３】本実施例におけるｓｉｎｃ関数の積分関数（正弦積分）と、複素正弦波荷重のｓｉｎｃ積分関数（ω＝２π／４）のグラフである。
【図４】本実施例における積分範囲を標本化メッシュに対して標本点の中間位置までシフトした場合の荷重分布と誤差分布のグラフである。
【図５】本実施例における移動有限フーリエ変換を実行する構成のブロック図である。
【図６】本実施例において、方程式の他の解法の窓関数と荷重関数の例を示すグラフである。
【符号の説明】
【００３７】
１０  乗算手段
２０  遅延部
２２  単位遅延手段
３０  累積部
３２  加算手段
【発明を実施するための最良の形態】
【００３８】
第１  解析信号の周波数等の推定
ａ）支配方程式
解析の対象となる信号を、解析信号と呼ぶ。まず、支配方程式について説明する。【００３９】
解析信号に対して下記式（１）
【数１】
2008136443.tif 000003
の一般解は、
【数２】
2008136443.tif 000004
となる。ここで、
【数３】
2008136443.tif 000005
であるから、その結果、下記式（２）
【数４】
2008136443.tif 000006
が得られる。すなわち、上で述べた式（１）（数１）が支配方程式である。複素正弦波の周波数は、この式中のａの虚部に相当し、その振幅変化率はａの実部に表現される。この支配方程式が、請求の範囲の「微分を含む支配方程式」の好適な一例に相当する。【００４０】
ｂ）荷重積分観測方程式
観測期間を［−Ｔ／２，Ｔ／２］（すなわり、総時間はＴ）とし、その期間中、信号周波数は一定であると仮定する。すなわり、上記ａ）で述べた微分方程式が一様に成立していると仮定する。この条件は、任意の荷重関数
【数５】
2008136443.tif 000007
を導入すると、下記の同値関係式（３）によって、積分形式に置き換えることができる。【００４１】
【数６】
2008136443.tif 000008
この支配方程式が、請求の範囲の「積分形式に変換した支配方程式」の好適な一例に相当する。この式（３）において、ｗ（ｔ）は、具体的には、区間［−Ｔ／２，Ｔ／２］で完備な関数系
【数７】
2008136443.tif 000009
を荷重関数として用いればよい。【００４２】
この式（３）の下で微分を消去するために、式（３）の後段の積分を実行すると、
【数８】
2008136443.tif 000010
となる。ここで、
【数９】
2008136443.tif 000011
と選ぶと、これらは区間［−Ｔ／２，Ｔ／２］で完備な直交系をなし、
【数１０】
2008136443.tif 000012
のように、荷重の微分による荷重積分が原関数の荷重積分に比例する。さらに、
【数１１】
2008136443.tif 000013
のように、全ての
【数１２】
2008136443.tif 000014
に関する積分境界値項が符号を除いて同一となる。ここで、積分境界値項とは積分区間の信号やその微分の両端の値の差を言う。以下、同様である。さてこの結果、
【数１３】
2008136443.tif 000015
及び、荷重積分を、下記式（４）に示すように
【数１４】
2008136443.tif 000016
と置くことによって、題意の荷重積分式は
【数１５】
2008136443.tif 000017
を未知変数とする線形方程式（下記式（５））
【数１６】
2008136443.tif 000018
に帰着させることができる。ａを任意の複素数としても未知数は４個であり、方程式は上記の実部と虚部の２個である。したがって、上記の線形方程式は、任意の２個の周波数で連立することによって解くことが可能になる。【００４３】
ｃ）ＤＣ成分との連立法
さて、具体的にω＝０の式と上式とを連立させると、
【数１７】
2008136443.tif 000019
より、下記の式（６）
【数１８】
2008136443.tif 000020
と解かれる。【００４４】
ここで注意すべきことは、ｆ（ｔ）は解析の対象である解析信号だということである。したがって、ｆ（ｔ）は虚部を有し、
【数１９】
2008136443.tif 000021
にはこのｆ（ｔ）のヒルベルト変換のフーリエ変換が加わり、
【数２０】
2008136443.tif 000022
も虚部を有する。解析信号化が上記の荷重積分で同時になされると考えることも可能ではあるが、それは何を解析信号であると認識するかの問題でもある。【００４５】
なお、ヒルベルト変換のフーリエ変換が加わっても
【数２１】
2008136443.tif 000023
は２倍されるだけであろうかという疑問が生じるかもしれないが、正しいヒルベルト変換ではないので２倍にはならない点にも注意すべきである。【００４６】
ｄ）近接２周波数成分を用いた連立法
今、利用する２個の周波数をω１、ω２とおき、それらの周波数におけるｓｔをｓ１，ｓ２とおくと、連立方程式は
【数２２】
2008136443.tif 000024
となる。したがって、
【数２３】
2008136443.tif 000025
となり、よって、下記（７）式、（８）式
【数２４】
2008136443.tif 000026
を得る。【００４７】
ｅ）振幅と位相の推定
さて、推定された周波数ａとして
【数２５】
2008136443.tif 000027
とおくと（ここでＡは複素振幅を表す）、
【数２６】
2008136443.tif 000028
より下記（９）式
【数２７】
2008136443.tif 000029
が求められる。ここで、ａＴ＝２ｎπのときはｓｉｎ（ａＴ／２）＝０となって適用不可能となるが、この場合は、単一周波数スペクトルとなって、振幅も位相もこのスペクトル値から決定することができる。【００４８】
第２  実信号の周波数等の推定
ａ）支配方程式
解析の対象となる実信号に関して、まず上の記載と同様に、支配方程式についての説明を行う。実信号に対して下記式（１０）
【数２８】
2008136443.tif 000030
の一般解は、周知のように、下記式（１１）
【数２９】
2008136443.tif 000031
となる。ここで、ＡやＢは符号も含めて任意である。すなわち、上で述べた式（１０）（数２８）が支配方程式であり、上記式（１１）の一般解の形の意味でのａの推定が周波数を推定することに相当する。この支配方程式が、請求の範囲の「微分を含む支配方程式」の好適な一例に相当する。ここでａの符号は決定しない、又は意味を有しないと言い換えてもよい。【００４９】
ｂ）荷重積分観測方程式
観測期間を、上の記載と同様に、［−Ｔ／２，Ｔ／２］（すなわち、総時間はＴ）とし、その期間中、信号周波数は一定、すなわち、上記ａ）で述べた微分方程式（支配方程式）が一様に成立していると仮定する。この条件は、上述したのと同様に、任意の荷重関数
【数３０】
2008136443.tif 000032
を導入すると、下記の同値関係式（１２）によって、積分形式に置き換えることができる。【００５０】
【数３１】
2008136443.tif 000033
この支配方程式が、請求の範囲の「積分形式に変換した支配方程式」の好適な一例に相当する。この式（１２）（数３１）において、上述した荷重関数は、具体的には、区間［−Ｔ／２，Ｔ／２］で完備な関数系
【数３２】
2008136443.tif 000034
を用いればよい。一般に、
【数３３】
2008136443.tif 000035
であるので、前節で述べたように、
【数３４】
2008136443.tif 000036
とおき、
【数３５】
2008136443.tif 000037
とおくと、前節で述べたのと同様の論法によって、同値条件を構成する周波数ごとの方程式として、下記式（１３）
【数３６】
2008136443.tif 000038
が得られる。【００５１】
なお、支配方程式を荷重関数を用いて積分形式にすることは上で述べたような機械的な作業である。特に、荷重関数の周期を積分区間の整数分の一にとっていることも、単純な作業の実現に寄与している。そのため、コンピュータ等の上で実施することが好適である。このようにして構成した式に、実際に計測したデータを導入すれば、正弦波のパラメータを変数として含む方程式が完成する。【００５２】
この処理に際して、まず支配方程式を決めて、それから荷重関数を導入して変形する処理をコンピュータの上で行わせるの好適であるし、また、積分形式の支配方程式の形を予め準備しておき、その係数を、荷重関数に応じて設定して方程式を構成することも好ましい。【００５３】
請求の範囲において方程式を「立てる」としている動作は、上述のように方程式を構成する処理を意味し、コンピュータで実行することが迅速な処理の為には好ましい。【００５４】
また、荷重関数の周期は、対象となる信号（正弦波）の周期と区別する為に、特に、荷重周期と呼ぶ場合がある。同様に、荷重関数の周波数は、対象となる信号の周波数と区別する為に、特に、荷重周波数と呼ぶ場合がある。【００５５】
ｃ）ＤＣ成分を含む連立法
さて、ω＝０に対する下記の関係を用いると、
【数３７】
2008136443.tif 000039
積分境界値が１個消去され、関係式は
【数３８】
2008136443.tif 000040
となる。この方程式は複素変数の方程式であり、
【数３９】
2008136443.tif 000041
の実部に注意し、上記方程式の実部のみを取り出すと、
【数４０】
2008136443.tif 000042
が得られる。よって、以下の式（１４）のように、
【数４１】
2008136443.tif 000043
と求められる。また、積分境界値は、
【数４２】
2008136443.tif 000044
のように書くことができる（式（１５）及び式（１６））。【００５６】
ｄ）ＤＣ成分を含む３周波数を用いた連立法
荷重積分観測方程式を解くもう一つの連立法は、上記ＤＣに加えてさらに２個の周波数ω１、ω２を連立する方法である。それぞれの周波数におけるｓｔをｓ１、ｓ２とおくと、これらを連立することによって
【数４３】
2008136443.tif 000045
であり、よってこの前段の第１式にｓ２ω２を乗じ、後段の第２式にｓ１ω１を乗じて差し引くことによって、
【数４４】
2008136443.tif 000046
のように、残りの積分境界値も消去され、よって、下記式（１７）
【数４５】
2008136443.tif 000047
が得られる。対称性が良く、従来の周波数成分の補間のほうほうと対応が取りやすいという利点を有する。【００５７】
ｅ）近接３周波数成分を用いた連立法
３個の周波数をω１、ω２、ω３とおき、そのときのｓｔをそれぞれｓ１、ｓ２、ｓ３とおくと、これらを連立することによって、
【数４６】
2008136443.tif 000048
を得る。この方程式の解は下記式（１８）、式（１９）、式（２０）のようになる。【００５８】
【数４７】
2008136443.tif 000049
ただし、
【数４８】
2008136443.tif 000050
である（式（２１））。なお、ここで、２×２と３×３の行列の逆行列の公式を用いている（下記数４９参照）。【００５９】
【数４９】
2008136443.tif 000051
また、念のため、ω３＝０、ｓ３＝１の場合について求めてみると
【数５０】
2008136443.tif 000052
となり、前節の結果と一致することが理解されよう。【００６０】
ｆ−１）振幅と位相の推定
さて、推定された周波数と、同時に推定された積分境界値項
【数５１】
2008136443.tif 000053
から、振幅・位相を求めることが考えられる。すなわち、ａ＞０を推定された周波数として、
【数５２】
2008136443.tif 000054
とおく。すると、
【数５３】
2008136443.tif 000055
であるので、
【数５４】
2008136443.tif 000056
すなわち、下記式（２２）
【数５５】
2008136443.tif 000057
を求めることができ、この結果から、振幅と位相が求められる。なお、数５３の導出には、下記三角関数の和と積の公式が用いられている。【００６１】
【数５６】
2008136443.tif 000058
さて、ａＴ＝２ｎπのときには、ｓｉｎ（ａＴ／２）＝０となってこの方法を適用することが困難であるが、この場合は、ωＴ＝２ｎπの中の単一周波数のみにスペクトルが立つ（表れる）ので、周波数も位相もこのスペクトルから決定可能である。【００６２】
したがって、アルゴリズムとしては、切り換えるための判断ステップ、又は、切り換えまで統合した推定アルゴリズムを用いることが好ましい。そのようなスペクトルのみが立つ場合を検出することはスペクトルを検査すれば比較的容易に判断できるので、その判断に基づきアルゴリズムを切り換えるのは当業者であれば容易に実施することができる。【００６３】
ｆ−２ａ）振幅と位相の他の推定方式
また、他の方式として、推定に用いた周波数成分
【数５７】
2008136443.tif 000059
へ正弦波モデル
【数５８】
2008136443.tif 000060
をフィッティングする方法もある。この方法は、上述した単一の周波数のみにスペクトルが立つ場合でも、条件判断をしてアルゴリズムを切り換える必要がないので便利である。また、全ての周波数成分の影響を受ける積分境界値を用いないため、複数の周波数が混合する信号を取り扱う場合にも高い精度を維持できると考えられる。【００６４】
すなわち、標記（でいいですよね）モデルの最小二乗評価関数は
【数５９】
2008136443.tif 000061
を展開すると、
【数６０】
2008136443.tif 000062
であり、Ａで微分して零とおくと、
【数６１】
2008136443.tif 000063
よって、
【数６２】
2008136443.tif 000064
となる。この式では、Ａは常に実数であるが、正の値だけでなく負にもなりえる点には注意すべきである。次に、φで微分して零とおくと、
【数６３】
2008136443.tif 000065
すなわち、
【数６４】
2008136443.tif 000066
すなわち、
【数６５】
2008136443.tif 000067
を得る。この式では、位相が［−π／２，π／２］の範囲でしか求まらないが、上記Ａの符号を正に固定し、符号を位相φに含ませると、全範囲で位相を推定することが可能である。【００６５】
ｆ−２ｂ）上記ｆ−２の別解
複素振幅をＡで表し、最小二乗評価関数
【数６６】
2008136443.tif 000068
をＡ、Ａ＊で微分して零とおくと、
【数６７】
2008136443.tif 000069
よって、
【数６８】
2008136443.tif 000070
となる。【００６６】
ｇ）複素正弦波の同時推定法（複数推定法）
観測信号（「解析信号」、「実信号」とも呼び、解析の対象となる信号を意味する）に複数の正弦波が含まれる場合について説明する。この観測信号に含まれる正弦波の個数Ｎは既知であるとする。このＮ個は予め知ることができる場合も多い。例えば、離散スペクトル分布のピークの個数から予めＮを求めることができる。【００６７】
このＮの算出もコンピュータ上で行うことが好ましい。スペクトル分布からピークを算出し、その個数をカウントすることはコンピュータによる一般的な信号処理として知られているので、そのような技術を用いればよい。【００６８】
さて、このＮ個の周波数を、
【数６９】
2008136443.tif 000071
とおく。この記法は、本節ｇ）においてのみ採用する記法であり、他の節では他の記法を採用する場合がある。また、正弦波を
【数７０】
2008136443.tif 000072
とおくと、
【数７１】
2008136443.tif 000073
と記述できる。これらによって表される各信号が満たすべき方程式は、
【数７２】
2008136443.tif 000074
である。よって、これに対応する荷重積分式は、
【数７３】
2008136443.tif 000075
である。この方程式は複素方程式であり、１個のωに対して２個の方程式を与える。一方、観測式としては全ての正弦波の混合信号に対する荷重積分
【数７４】
2008136443.tif 000076
しか得られない。この式も複素方程式であり、同様に２個以上のωの組に対して２個の方程式を与える。これに対して、未知数は、１個の正弦波に対して、周波数
【数７５】
2008136443.tif 000077
の１個と、実数である２種類の積分境界値である。したがって、この境界積分値の１＋１＝２個を足し、計３個であり、これに周波数Ｎ個を乗じて、合計３Ｎ個の未知数がある。ここで、未知数を周波数と振幅と位相で考えた場合も同様の結果になる。【００６９】
さて、これらの未知数と荷重積分は一意に結ばれるため、ここで、２個の方程式が立てられる。【００７０】
よって、推定に用いるωの個数をＭとすると、これらを未知数として解が得られるための必要条件は、下記式（２３）
【数７６】
2008136443.tif 000078
となる。Ｎ＝１の場合は、この式（２３）から、Ｍ＞３／２と求まる。すなわち、２個の周波数が必要で、前節の２周波数連立法の結果と一致することが理解されよう。【００７１】
なお、これらの方程式の構築は、コンピュータ上で実行することが好ましい。所定のωから方程式を得ることは上述のように機械的な作業であるので、汎用的なコンピュータ上でも容易に実行することが可能である。【００７２】
本特許において、方程式を「立てる」と述べている部分があるが、本特許で「立てる」と述べている作業は、上述のごとく機械的に、決まった形の方程式を「準備する」作業であるので、汎用的なコンピュータで十分に実行することが可能である。【００７３】
さて、標本値からのフーリエ係数の計算によってわずかに加わる誤差を無視すると、有限フーリエ変換は離散フーリエ変換で高速に実行することができる。このとき、例えば、２５６点のＦＦＴを使用することを前提とすれば、Ｍ（推定に用いるωの個数）＝１２７とすると（ここで、説明を簡単にするため、実の直流成分と最高周波数成分とは除く）、Ｎ＜（２Ｍ＋２）／３＝２５７／３＝８５．３３が得られる。したがって、計算上、８５個までの正弦波のパラメータを推定することができる。同様に、１２８点のＦＦＴでは４２個、６４点では２１個、３２点では１０個、１６点では５個、８点では２個、４点では１個の正弦波のパラメータを推定することができる。
【実施例】
【００７４】
以下、具体的な実施例を説明する。【００７５】
本発明の実施においては、直接的な解法は困難であると考えられるので、推定精度を必要なだけ向上させるための反復法を用いるのが実用的であると考えられる。【００７６】
「解法１」（反復法１）
ステップ１
まず、単一周波数を想定して全ての隣接３周波数荷重積分値から、周波数を推定する。【００７７】
ステップ２
次に、主要な周波数成分の周波数と位相とを推定する。【００７８】
ステップ３
上記の主要な周波数成分（周波数の摂動を含む）による荷重積分値の分布を算出する。【００７９】
ステップ４
観測した荷重積分値の分布との誤差によって推定値を更新する。【００８０】
ステップ５
上記のステップ３とステップ４を必要回数繰り返す。更新する量が一定値未満になった場合に所定の精度が得られたものと判断し、終了する。【００８１】
「解法２」（反復法２）
この方法の処理の流れのフローチャートが図１に示されている。【００８２】
ステップ１
まず、単一周波数を想定して全ての隣接３周波数荷重積分値から、周波数と振幅と位相とを推定する。この処理は、図１中、Ｓ１−１に相当する。【００８３】
ステップ２
次に、ピークを検出することによって、主要な周波数成分を取り出し（図１中ステップＳ１−２に相当する）、これら周波数成分の周波数と振幅と位相とを推定する（図１中のＳ１−３に相当する）。【００８４】
ステップ３
上記の主要な周波数成分の推定に用いた荷重積分値から他の主要な周波数成分の影響を差し引く。この処理は、図１中、Ｓ１−４に相当する。【００８５】
ステップ４
修正された荷重積分値から、新たに周波数と振幅と位相とを推定する。この処理は図１中、Ｓ１−５に相当する。【００８６】
ステップ５
上記のステップ２の後半（Ｓ１−３）と、ステップ３（Ｓ１−４）と、ステップ４（Ｓ１−５）を必要回数繰り返す。差し引く量が一定値未満になった場合に所定の精度が得られたものと判断し、終了する。この繰り返し判断処理は図１中、Ｓ１−６に相当する。【００８７】
実際にこのアルゴリズムをコンピュータ上でシミュレーションした結果が図２に示されている。図２は、Ｎ＝３２による単一並びに複数正弦波信号に対する周波数の推定結果を表すグラフである。このグラフにおいて、横軸は与えられた周波数（右上がり４５度の周波数成分である）であり、縦軸が推定された周波数である。【００８８】
図２（ａ）は１個の周波数を推定する場合であり、図２（ｂ）は、２個の周波数を推定する場合、そして、図２（４）は、４個の周波数を推定する場合である。【００８９】
全ての隣接する３荷重積分で単一周波数を仮定して独立に推定し、ｄｅｔＡの値が所定の第１しきい値以上のものを大きな黒点で表し、第１しきい値より小さな第２しきい値以上のものを小さな黒点で示している。複数の周波数成分が交差する付近ではうねり状の誤差が生じている。【００９０】
ｈ）基本値系列を用いる荷重積分の高精度化
一般にナイキスト（Ｎｙｑｕｉｓｔ）条件を満たして間隔Δｔで標本化された標本値の時系列
【数７７】
2008136443.tif 000079
は、ｓｉｎｃ関数補間された連続関数
【数７８】
2008136443.tif 000080
と等価である。よって区間［−Ｔ／２，Ｔ／２］の荷重積分は、
【数７９】
2008136443.tif 000081
となる。ただし、ωＴ＝２ｎπ（ｎ：整数）のように記述される。上式（数８２）中の積分のｔ＞０をグラフで表したものが、図３に示されている。【００９１】
この図３においては、ｓｉｎｃ関数の積分関数（正弦積分）と、複素正弦波荷重のｓｉｎｃ積分関数（ω＝２π／４）のグラフが示されている。実線がＳｉ（ｔ）であり、破線がＳｉｕ（ｔ）の実部であり、一点鎖線はＳｉｕ（ｔ）の虚部である。【００９２】
すなわち、ｔ＜０でほぼ零、ｔ＝０では１／２、ｔ＞０では、ほぼ１の値をとる。しかし、原点付近では波打があり、そのｔ＞０における１との交差、及び、ｔ＜０における０との交差は、ｔの整数点のほぼ中間点に位置する。このことから、積分範囲を標本点の中間に取ることによって精度を格段に向上させることが可能である。【００９３】
なお、本特許では、解析の対象である信号は、所定のサンプリング周期でサンプリング（標本化）することを前提としている。【００９４】
正弦積分の整数点上での値は、原点（ｔ＝０）からそれぞれ、
０．００００００（ｔ＝０）
０．５８９４９０（ｔ＝１）
０．４５１４１２（ｔ＝２）
０．５３３０９４（ｔ＝４）
０．４７４９７０（ｔ＝５）
０．５２０１０８（ｔ＝６）
０．４８３２０６（ｔ＝７）
となる。中間点での値は、ｔ＝０．５からそれぞれ、
０．４３６３２８（ｔ＝０．５）
０．５１１９６１（ｔ＝１．５）
０．４９５２３７（ｔ＝２．５）
０．５０２５１９（ｔ＝３．５）
０．４９８４５２（ｔ＝４．５）
０．５０１０４７（ｔ＝５．５）
である。【００９５】
このように、積分範囲を標本化メッシュに対して標本点の中間位置までシフトした場合の荷重分布と誤差分布のグラフが図４に示されている。【００９６】
図４（１）は、積分境界を標本化メッシュの中間に取った場合の誤差分布を表すグラフである。その横軸は、サンプリング点を表し、縦軸は誤差を表す。また、図４（２）は、積分境界を標本化メッシュの中間に取った場合の荷重分布を表すグラフである。その横軸は、サンプリング点を表し、縦軸は荷重を表す。【００９７】
この図４に示すように、荷重分布のばらつきは大きく減少し、境界に隣接する標本点に小さな誤差が残るのみとなる。この荷重は、境界の外側隣接点で０．５−０．４３６３２８＝０．０６３６７２となり、内側隣接点で０．５＋０．４３６３２８＝０．９３６３２８となる。上記の補償は−１、０及びＮ−１、Ｎの標本値に該当し、適用には前後に計２画素（サンプル）の拡張が必要である。補償については次に述べる。【００９８】
補償処理
さて、このような積分範囲のシフト（積分範囲を標本点の中間位置までずらす）の影響を補償するために、以下のような処理を行う。まず、積分領域の標本化メッシュに対するシフト量をεとすると、積分境界値項は
【数８０】
2008136443.tif 000082
のように位相シフト項
【数８１】
2008136443.tif 000083
が乗じられる。よって、標本化間隔Δｔを用いて、補償量は、具体的には、
【数８２】
2008136443.tif 000084
となることが理解されよう。【００９９】
ｈ）推定法の性質についての検討
以上、種々の推定法についてそれぞれ検討を行ったが、推定法は以下のような性質を有していると想像される。【０１００】
１．無雑音、単一周波数の正弦波で、積分によるフーリエ係数を用いる場合は、厳密解が得られる。【０１０１】
２．標本値を用いた離散フーリエ変換を利用しても、標本化密度を十分に高めることによって、厳密解にいくらでも近づけることが可能である。【０１０２】
３．無雑音、含まれる周波数の個数が既知、積分によるフーリエ係数を観測地に用いる場合は厳密解が得られる。【０１０３】
第３  正弦波周波数推定の応用
以上述べたように、本実施の形態によれば、信号中に含まれる正弦波のパラメータを推定（周波数、振幅、位相）できる。このような正弦波の推定は、音声信号処理の基本処理とも言うべきものであり、その応用範囲は極めて広い。応用できると考えられる分野、用途、目的としては以下のようなものが考えられる。【０１０４】
・分析区間によらない厳密な周波数推定を測定。【０１０５】
・調和関係にない複数正弦波の周波数の測定
・一般調和解析
・音声、楽器音、異常音の時間周波数分析
・音の立ち上がり、立ち下がりの分析
・時間応答による残響分析
・周波数精度と周波数分解能の明確な概念的な分離
・周波数変調の復調、位相変調の微分復調
・ドップラ検出、反射波による振動検出、音源の振動検出
本発明は、音声信号の解析だけでなく、その他各種信号に応用できることは言うまでもない。また、本特許においては、正弦波のパラメータとは、周波数と振幅と位相を意味する。【０１０６】
なお、正弦波のパラメータという場合、周波数と振幅と位相とに加えて直流バイアスを意味する場合もある。【０１０７】
本特許においては、この直流バイアスを直接パラメータとして扱ってはいないが、周波数が直流という特別の場合の振幅として推定することができる。【０１０８】
第４  移動フーリエ変換
時間積分型音源恒等式に必要な荷重積分、すなわち有限フーリエ変換には、窓関数が不要という顕著な特徴がある。これを数値計算法やハードウェア信号処理に利用することを検討する。【０１０９】
上述した特徴は重なりのある密な短時間フーリエ変換の計算の場合に特に顕著に発揮される特徴である。すなわち、この計算は、各標本点で複素指数関数
【数８３】
2008136443.tif 000085
を乗じ、これらを区間幅Ｎで移動平均することと備えられる。荷重が「１」で一定であるため、この移動平均は、各標本点で、区間の前端と後端の差を累積メモリに加算するだけで実行可能である。これを、ここでは、移動（有限）フーリエ変換と称することにする。【０１１０】
このときに生じる分析区間に応じた初期位相のずれは、後で補償すれば問題はない。【０１１１】
４．１  演算量の評価
移動フーリエ変換の乗算回数と加算回数（複素数）は、各標本点でＮ＋Ｎ＝２Ｎ回である。全データ長さをＭとすると、重なりなしのＦＦＴで乗算回数は
【数８４】
2008136443.tif 000086
であり、完全に密な移動演算とすると、
【数８５】
2008136443.tif 000087
である。これに対して上記の移動フーリエ変換法の乗算回数は、
【数８６】
2008136443.tif 000088
であるので、ＦＦＴをそのまま用いるよりも少ない値となる。もし、窓関数が必要であれば、窓付きでＦＦＴを完全に密に適用するには、乗算回数
【数８７】
2008136443.tif 000089
が必要である。これに対して上記の移動フーリエ変換法の場合は、
【数８８】
2008136443.tif 000090
であり、やはり少なくなる。【０１１２】
次に、密な窓関数付きＦＦＴと、窓関数なしの移動フーリエ変換法との比較では、乗算回数では
【数８９】
2008136443.tif 000091
倍の開きがある。Ｎ＝２５６の場合は、上記値は、２５６×８＝２０４８となり、２０４８倍の開きとなる。また、窓関数なしのＦＦＴと移動フーリエ変換法とでは
【数９０】
2008136443.tif 000092
倍の開きがある。もし、Ｎ＝２５６の場合は、上記値は８となり、８倍の開きとなる。【０１１３】
すなわち、従来法（例えばフーリエ位相法）で密に遅延推定を行う場合に比べて３桁の高速化が可能となると期待される。また、音源恒等式を密にＦＦＴを用いて行う場合と比較しても、少なくとも１桁の高速化が期待される。【０１１４】
４．２  移動平均の漸化式の計算法
一般に平均化の荷重分布を
【数９１】
2008136443.tif 000093
とすると、区間［−Ｔ／２、Ｔ／２］の移動平均は、
【数９２】
2008136443.tif 000094
であり、
【数９３】
2008136443.tif 000095
である。ところが、一方、
【数９４】
2008136443.tif 000096
である。ここで、
【数９５】
2008136443.tif 000097
は区間内及びその両端で連続かつ微分可能である。したがって、
【数９６】
2008136443.tif 000098
では、
【数９７】
2008136443.tif 000099
となる。すなわち、
【数９８】
2008136443.tif 000100
のような漸化式による計算が可能である。この漸化式の右辺の中括弧は単に積分の境界における信号値の差であり、区間内の値を用いずに容易に計算できる。【０１１５】
一方、離散値を用いる場合は、上記の漸化式は厳密な式となり、厳密解が得られる。その理由は、区間［０，Ｎ−１］の
【数９９】
2008136443.tif 000101
の移動平均
【数１００】
2008136443.tif 000102
は
【数１０１】
2008136443.tif 000103
であり、
【数１０２】
2008136443.tif 000104
とすると、
【数１０３】
2008136443.tif 000105
であり、
【数１０４】
2008136443.tif 000106
ならば
【数１０５】
2008136443.tif 000107
となるからである。なお、数値計算誤差の累積を防ぐ手段を講じることが好ましい。具体的には、漸化式の直流応答を有限なものとする等の手段を高じることが好ましい。このような手段は従来から知られている工夫であるから、当業者であれば容易に実行することができる。【０１１６】
４．３  移動有限フーリエ変換の実施の一例
さて、これまで説明してきた事項に基づき、移動有限フーリエ変換を実行する構成のブロック図が図５に示されている。【０１１７】
この図５に示すように、入力データ（観測信号）
【数１０６】
2008136443.tif 000108
に
【数１０７】
2008136443.tif 000109
を乗じてＮ系統、Ｎ段の遅延部に入力する。【０１１８】
図５におけるＮ個の乗算手段１０は、入力した各データに上述した荷重を乗じる。乗算した結果は、遅延部２０に供給される。この遅延部２０は、単位遅延手段２２をＮ段接続して構成されている。【０１１９】
遅延部２０の先頭と終端の差が累積部３０入力となり、各時刻の累積値が各時刻のＳＴＦＴ出力を与える。ここで遅延部２０の先頭と終端の差と、それまでの累積値との加算が加算手段３２によって行われる。遅延部２０の先頭のデータすなわち入力しようとするデータは符号がそのままで加算手段３２に供給され、遅延部２０の終端のデータは符号が反転されてから加算手段３２に供給される。また、累積部３０の出力信号がさらに加算手段３２に供給される。このような構成によって、加算手段３２は先頭のデータから終端のデータを減算して「差」を求めてこの差と、累積部３０の現在の累積値とを「累積」して、その結果を累積部３０に再び格納する。累積部３０はいわゆるアキュムレータから構成される。この累積部３０の出力信号がＳＴＦＴ出力を与える。【０１２０】
図５に示す形態のままで、前端付近と後端付近の荷重補償も行うことが可能である。また、矩形以外の窓関数では、非零の窓関数の微分の個数だけ遅延部からの乗算と加算の段数が増えることに注意されたい。演算量の低減は、窓関数が不要（境界以外で一様に−１）であるために実現されていることに注意すべきである。【０１２１】
４．４  分析区間の初期位相の補償
さて、ＳＴＦＴによって求められるのは、
【数１０８】
2008136443.tif 000110
である。これに対して、移動フーリエ変換で求められるのは、
【数１０９】
2008136443.tif 000111
であり、ｔ＋τを新たにｔと置き換えることによって、
【数１１０】
2008136443.tif 000112
のように互いに変換可能であることが理解されよう。この際の補償位相は分析の離散周波数ごとに異なることに注意するべきである。積分範囲における１／２標本点シフト（上述したようなサンプリング周期の１／２の期間のシフト）のために位相補償の乗算は元々必要である。したがって、これに含ませることによって、この演算は演算量の増大には結びつかない。【０１２２】
第５．参考発明
これまで、近似ではなく厳密にパラメータを求めることができる正弦波パラメータ推定方法を説明したが、この方法は、以下のようなハードウェアを用いて実施することが好適である。【０１２３】
５．１  装置
上で述べた正弦波パラメータ推定方法を実施する為に、任意の荷重関数を導入した有限区間の積分形式で、前記信号の支配方程式を構成する手段と、前記方程式を解く手段と、を含むことを特徴とする正弦波パラメータ推定支援装置を用いることが好ましい。【０１２４】
このような手段は、例えば、デジタルシグナルプロセッサ等で実現することが好適である。例えば、ＦＦＴを実行するする半導体チップや、暗号化処理を実行する暗号チップ等が知られているので、このような半導体チップと同様に所定の処理を高速に実施する半導体チップを作成することが好ましい。【０１２５】
なお、種々のアプリケーションに対応する為に、種々のサンプリング周期や、様々な音源の数に対応する為には、ソフトウェアとそのソフトウェアを実行する汎用的なプロセッサとで実現することも好ましいが、速度的には遅くなってしまう可能性がある。【０１２６】
また、音源の数や、データの入出力部分や、サンプリング数の設定部分のみをソフトウェアで構成し、方程式を解く演算部分のみを半導体チップを用いて構成することも好適である。このように、一部をハードウェアで、他の一部をソフトウェアで構成すれば、種々のアプリケーションに対応できると共に高速なデータ処理（信号処理）を実現できる可能性がある。【０１２７】
５．２  具体的な構成
また、参考発明としては、これまで述べた正弦波パラメータ推定方法を支援する装置において、前記解析信号に前記荷重を乗算するＮ個の乗算手段と、前記Ｎ個の乗算手段毎にそれぞれ設けられ、前記乗算手段の各出力信号をＮ周期遅延させるＮ個の遅延手段と、前記Ｎ個の遅延手段毎にそれぞれ設けられ、遅延後の出力信号と、遅延前の信号との差を求め、その差を累積値に累積するＮ個の加算手段と、前記加算手段毎にそれぞれ設けられ、前記加算手段の出力信号を累積するＮ個の累積手段と、を含むことを特徴とする正弦波パラメータ推定方法支援装置である。ここで、前記Ｎは正の整数である。【０１２８】
このような構成は、図５に示した通りである。【０１２９】
第６  方程式の他の解法
上記第２章（ｃ）ＤＣ成分を含む連立法や、同（ｄ）ＤＣ成分を含む３周波数を用いた連立法等において、方程式の種々の解法を説明した。しかし、連立方程式の解法は、他にもあり、以下、他の解法を説明する。【０１３０】
ａ）  単一荷重周波数と窓関数を用いる方法
本節でも、荷重積分観測値の種類が増加する欠点を承知の上で、積分境界値を未知数から除去する方法を導くことにする。出発点となる荷重積分観測方程式
【数１１１】
2008136443.tif 000113
において、
【数１１２】
2008136443.tif 000114
かつ、この窓関数ｐ（ｔ）を原点対称に選ぶと、
【数１１３】
2008136443.tif 000115
である。ゆえに、荷重積分観測方程式は
【数１１４】
2008136443.tif 000116
となる。ただし、
【数１１５】
2008136443.tif 000117
である。特に
【数１１６】
2008136443.tif 000118
となるように、窓関数を選ぶと、荷重積分観測方程式は、
【数１１７】
2008136443.tif 000119
となり、
【数１１８】
2008136443.tif 000120
のように求められる。ｐ（ｔ）の具体例としてＨａｎｎの窓関数を利用することが考えられる。Ｈａｎｎの窓関数を利用し、
【数１１９】
2008136443.tif 000121
とすると（ここで、
【数１２０】
2008136443.tif 000122
である）、原点対象で
【数１２１】
2008136443.tif 000123
となる。また、
【数１２２】
2008136443.tif 000124
である。これらのグラフが図６に示されている。図６のグラフは、横軸が時間を表し、縦軸が荷重を表す。図６（ａ）は、窓関数を示し、図６（ｂ）は荷重関数（Ｔ＝８）を表す。図６（ａ）では、実線が
【数１２３】
2008136443.tif 000125
を表し、破線が
【数１２４】
2008136443.tif 000126
を表し、一点鎖線が
【数１２５】
2008136443.tif 000127
をそれぞれ表す。図６（ｂ）は、これらに、ｎ＝４のｃｏｓ正弦波荷重関数を乗じたものである。
【図１】
2008136443.tif 000128
【図２】
2008136443.tif 000129
【図３】
2008136443.tif 000130
【図４】
2008136443.tif 000131
【図５】
2008136443.tif 000132
【図６】
2008136443.tif 000133
【国際調査報告】
2008136443.tif 000134
2008136443.tif 000135
2008136443.tif 000136
2008136443.tif 000137
2008136443.tif 000138
2008136443.tif 000139
2008136443.tif 000140
2008136443.tif 000141