JP3986457B2 - Input signal estimation method and apparatus, input signal estimation program, and recording medium therefor - Google Patents

Description

但し、ａ（ｚ）＝ａ₁ ｚ^-1＋ａ₂ ｚ^-2＋…＋ａ_p ｚ^-p
または、
【数１】

Ｃ＝Ｅ｛ｕ _n ^T ｕ _n-1 ｝（Ｅ｛ｕ _n-1 ^T ｕ _n-1 ｝）^-1 （2c）
但し、ｅ（ｎ）：線形予測誤差， ^T：転置
【０００６】
入力端１２とセンサ１３₁ ，１３₂ との間の各伝達関数ｈ₁(ｚ），ｈ₂(ｚ）（ｚ：複素変数）及びフィルタ１５₁ ，１５₂ の各伝達関数ｗ₁(ｚ），ｗ₂(ｚ）をそれぞれ次の様に表す。

但し、ｈ_j(z)（ｊ＝１，２）は共通零点を持たない。
【０００７】
この時、減算器１７の出力信号ｙ（ｎ）は次の様に表される。
【数２】

フィルタ１５_j の次数Ｌと系１１の伝達関数の次数ｍとの関係をＬ≧ｍ−１と定め、減算器１７の出力信号ｙ（ｎ）の二乗平均値Ｅ｛｜ｙ（ｎ＋１）｜² ｝を最小化すると、式（２ｃ）より、フィルタ１５₁ ，１５₂ の伝達関数を表すベクトルｗは次の様に求められる。

但し、ｓ² ：小さな正数、Ｉ：単位行列
【０００８】
ここで（ＨＥ｛ｕ _n-1 ^T ｕ _n-1 ｝Ｈ ^T ＋ｓ² Ｉ）^-1：行列ＨＥ｛ｕ _n-1 ^T ｕ _n-1 ｝Ｈ ^T の擬似逆行列
ｓ² を十分小さな正数とする事より、式（５）の関係は次の様に整理される（非特許文献１、３３９頁）。
【数３】

但し、Ａ ⁺ ：行列ＡのMoore-Penrose一般逆行列，Ｕ _n-1 ^T Ｕ _n-1 ＝Ｅ
｛ｕ _n-1 ^T ｕ _n-1 ｝
式（５）を式（４）へ代入し式（２ｂ）の関係を用いれば、出力信号ｙ（ｎ）は次の様に求められる。

【０００９】
このようにｙ（ｎ）はｕ（ｎ）とならない、つまりＡＲ過程より生成される入力信号ｕ（ｎ）は特許文献１に示す技術では推定されない。
この発明の目的は、室内音声収音等の１入力多出力線形系において、系の伝達関数に関わる事前情報を用いず、各センサの出力信号のみから、ＡＲ過程より生成される入力信号でも高精度にブラインド（Blind）推定する方法、及び装置、そのプログラム及び記録媒体を提供することである。
【００１０】
【課題を解決するための手段】
この発明によれば、１入力多出力線形系の出力端に備える第１，第２，…，第Ｎのセンサの出力信号各々に１サンプルの遅延を付与し、上記Ｎ個の出力信号と上記Ｎ個の遅延信号を用いてＮ個の前段フィルタ伝達関数と１個の後段フィルタ伝達関数を計算し、上記Ｎ個の遅延信号をそれぞれ上記Ｎ個の前段フィルタ伝達関数をフィルタ処理し、これらフィルタ処理した信号を加算し、この加算信号を前記第１センサの出力信号から減算し、この減算信号を上記後段フィルタ伝達関数でフィルタ処理する。
【００１１】
フィルタ伝達関数の計算は例えば次のように求めることができる。Ｎ個の遅延信号の自己相関行列を求め、更にこの自己相関行列の擬似逆行列を求め、また上記Ｎ個の遅延信号と上記Ｎ個のセンサの各出力信号との相互相関行列を求め、上記自己相関行列の擬似逆行列と上記相互相関行列との積行列を求める。この積行列の第１行のＮ個の成分を上記Ｎ個の前段フィルタ伝達関数とする。また、上記積行列の特性方程式を求め、これより上記ＡＲ生成過程の分母多項式と同等な多項式±｛１−ａ（ｚ）｝を求め、この多項式の逆元±１／（１−ａ（ｚ））を求めて後段フィルタ伝達関数とする。この後段フィルタ２２で前記減算信号ｈ_1:0 ｅ（ｎ）が処理され、その処理結果としてｅ（ｎ）は式（２ａ）の関係から、次式（８）に示す様に入力信号ｕ（ｎ）に比例する信号ｈ_1:0 ｕ（ｎ）が得られる。
±１／（１−ａ（ｚ））ｈ_1:0 ｅ（ｎ）＝±ｈ_1:0 ｕ（ｎ） （８）
即ち、系への入力信号は高精度に推定される。
【００１２】
【発明の実施の形態】
以下、図面を用いてこの発明の実施の形態を説明する。この発明の入力信号推定装置の実施形態の機能構成を図１に示す。１入力２出力線形系１１の入力端１２に入力された、時刻ｎの入力信号ｕ（ｎ）が出力端の２個のセンサ１３₁ ，１３₂ で電気信号として検出され、これら各センサ出力信号は遅延器１４₁ ，１４₂ で１サンプル遅延され、これら遅延信号は前段フィルタ２１₁ ，２１₂ でフィルタ処理され、これらフィルタ出力信号は加算器１６で加算されこの加算信号が第１センサ１３₁ の出力信号より減算器１７で差し引かれ、この減算信号が後段フィルタ２２でフィルタ処理される。フィルタ計算部２３でセンサ１３₁ ，１３₂ の各出力信号と遅延器１４₁ ，１４₂ の各遅延信号に基づき、前段フィルタ伝達関数ｗ ₁ ，ｗ ₂ 、後段フィルタ伝達関数を計算し、これら伝達関数が前段フィルタ２１₁ ，２１₂ 、後段フィルタ２２にそれぞれセットされる。
【００１３】
フィルタ計算
フィルタ計算部２３の機能構成例を図２に示す。
遅延器１４₁ ，１４₂ の各出力信号より自己相関計算部２４で自己相関行列Ｒ _A が計算され、この擬似逆行列（Ｒ _A ＋ｓ² Ｉ）^-1が擬似逆行列計算部２５で計算される。遅延器１４₁ ，１４₂ の各出力信号とセンサ１３₁ ，１３₂ の各出力信号との相互相関行列Ｒ _c が相互相関計算部２６で計算される。
【００１４】
上記擬似逆行列Ｒ _A ^-と相互相関行列Ｒ _c との積行列Ｗ＝Ｒ _c(Ｒ _A ＋ｓ² Ｉ）^-1が乗算部２７で計算される。この積行列Ｗの第１行の２個の成分が前段フィルタ伝達関数ｗ ₁ ，ｗ ₂ として前段フィルタ２１₁ ，２１₂ にセットされる。更に積行列Ｗの特性多項式が多項式計算部２８で計算され、その特性多項式中の、ＡＲ生成過程の分母多項式と同等な多項式の逆元が逆元計算部２９で計算され、後段フィルタ伝達関数ｗ_F として後段フィルタ２２にセットされる。
【００１５】
次にフィルタ計算部２３での処理を詳細に説明する。系１１への入力信号ｕ（ｎ）、及びｕ（ｎ）を生成するＡＲ過程、系１１の伝達関数ｈ_j(ｚ）、前段フィルタ２１ _j(ｊ＝１，２）の伝達関数ｗ_j(ｚ）、及びｈ_j(ｚ）の次数ｍとｗ_j(ｚ）の次数Ｌとの関係は、それぞれ前記の式（２）〜（５）の場合と同じとする。遅延器１４₁ ，１４₂ の出力信号より得られる自己相関行列Ｒ_A は、次式（９）の計算により自己相関計算部２４の計算で求められる。
Ｒ_A ＝ＨＥ｛ｕ_n-1 ^T ｕ_n-1 ｝Ｈ^T （９）
擬似逆行列Ｒ_A ^-を逆行列計算部２５で求めるとＲ_Aは（Ｒ_A ＋ｓ²Ｉ）^-1となり、この擬似逆行列に、式（９）を代入すると、次式（１０）に示すように式（５）中の逆行列の部と同一になる。
Ｒ_A ^-＝(Ｒ_A＋ｓ²Ｉ)^-1＝(ＨＥ｛ｕ_n-1 ^Tｕ_n-1｝Ｈ^T＋ｓ² Ｉ)^-1 （10）
【００１６】
相互相関行列Ｒ _C は、遅延器１４₁ ，１４₂ 、及びセンサ１３₁ ，１３₂ の各出力より、次式（１１）の計算により相互相関計算部２６で求められる。
Ｒ _C ＝ＨＥ｛ｕ _n-1 ^T ｕ _n ｝Ｈ ^T （11）
よって、乗算部２７で計算される積行列Ｗは次式（１２）の様になる。

式（１２）に式（２ｃ）の関係を代入し、ｓ² を十分小さな正数とする事より、式（６）と同様に、式（１２）の関係は次の様に整理される。
【数４】

【００１７】
この積行列Ｗの第１行Ｗ _1rの各成分が先に前段フィルタ２１₁ ，２１₂ にそれぞれ伝達関数としてセットされるが、この積行列Ｗの第１行は従来技術の項で述べたサンプル遅延信号をフィルタ処理するフィルタ１５₁ ，１５₂ の伝達関数ベクトルｗ、つまり式（６）に等しい。

従って、遅延器１４₁ ，１４₂ の出力遅延信号を前段フィルタ２１₁ ，２１₂ でそれぞれフィルタ処理し、これら前段フィルタ２１₁ ，２１₂ の出力信号の加算信号を第１センサ１３₁ の出力信号より減算して得られる減算器１７の出力信号ｙ（ｎ）は系１１への入力信号ｕ（ｎ）の生成過程がＡＲ過程である場合、式（７）に示したように式（２ａ）で示されるその入力信号ｕ（ｎ）を生成するＡＲ過程１／（１−ａ（ｚ））を駆動する線形予測誤差ｅ（ｎ）に比例した信号ｈ_1:0 ｅ（ｎ）になる。
【００１８】
ところで積行列Ｗの零でない固有値λ_i｛Ｗ｝は次の関係を満足する（非特許文献１、１６６〜１６７頁）。

行列Ｃの特性多項式ｃ_c(ｚ)は、式（２ｂ）より、次の様に求まる。

式（１５）の関係があるから、積行列Ｗは行列Ｕに対し１行の要素が２（Ｌ＋１）−ｐだけ多いから、Ｗの特性多項式ｃ_w(ｚ)は次の様に表される。

【００１９】
上記特性多項式ｃ_w(ｚ)にｚ^-2(L+1)を乗ずる事より、式（2a）に示すＡＲ生成過程の分母多項式に同等な多項式±｛１−ａ（ｚ）｝を得ることができる。実際には多項式計算部２８で積行列Ｗの特性多項式ｃ_w(ｚ)を計算し、その中のＡＲ生成過程の分母多項式と同等な多項式±｛１−ａ（ｚ）｝を取り出し、この多項式の逆元±１／（１−ａ（ｚ））を逆元計算部２９で計算する。
この逆元±１／（１−ａ（ｚ））＝ｗ_F を後段フィルタに伝達関数ｗ_F としてセットする。減算器１７の出力信号ｈ_1:0ｅ(ｎ）は後段フィルタ２２でフィルタ処理される。この処理により式（８）が演算され、系１１の入力信号ｕ（ｎ）に比例する信号±ｈ_1:0ｕ(ｎ）が得られる。即ち、系１１への入力信号が高精度に推定される。
【００２０】
次にこの発明の方法を図１に示した１入力２出力線形系１１に適用した実施形態を図３を参照して簡単に説明する。この方法の各手順を、図４に簡略に示すコンピュータにより実行させることを想定する。
センサ１３₁ ，１３₂ の出力信号を図に示してないが通常はそれぞれ一定周期で標本化してディジタル信号系列としてインタフェース３１を通して記憶部３２に一旦取り込む（Ｓ１）。これらセンサ出力信号に対し、それぞれ１サンプルの遅延を付与し（Ｓ２）、これら遅延信号と、センサ出力信号とに基づき前段フィルタ伝達関数ｗ ₁ ，ｗ ₂ 及び後段フィルタ伝達関数ｗ_F を計算する（Ｓ３）。
【００２１】
つまり例えば図２に示した処理と同様に両遅延信号の自己相関行列Ｒ _A ＝ＨＥ｛ｕ _n-1 ^T ｕ _n-1｝Ｈ ^T を計算し（Ｓ３−１）、更にその擬似逆行列Ｒ _A ^-＝（Ｒ _A ＋ｓ² Ｉ）^-1を計算して一時記憶する（Ｓ３−２）。また両センサ出力信号と両遅延信号との相互相関行列Ｒ _C ＝ＨＥ｛ｕ _n-1 ^T ｕ _n ｝Ｈ ^T を計算し（Ｓ３−３）、相互相関行列Ｒ _Cと擬似逆行列Ｒ _A ^-との積行列Ｗを計算する（Ｓ３−４）。この積行列Ｗの第１行の各成分をそれぞれ前段フィルタ伝達関数ｗ ₁ ，ｗ ₂ として記憶する（Ｓ３−５）。
【００２２】
また積行列Ｗの特性多項式ｃ_w(ｚ)を計算し（Ｓ３−６）、その特性多項式中のＡＲ生成過程の分母多項式に同等な多項式の逆元±１／（１−ａ（ｚ））を計算し、これを後段フィルタ伝達関数ｗ_F として記憶する（Ｓ３−７）。
このようにして前段フィルタ伝達関数と後段フィルタ伝達関数を求めた後、各遅延信号を対応する前段フィルタ伝達関数ｗ ₁ ，ｗ ₂ を用いてフィルタ処理し（Ｓ４）、これら前段フィルタ処理された信号を加算し（Ｓ５）、この加算信号を第１センサ出力信号から差し引き（Ｓ６）、その差し引かれた信号に対し後段フィルタ伝達関数ｗ_F でフィルタ処理し（Ｓ７）、その結果を出力部３３を通じて出力する（Ｓ８）。
【００２３】
この図３に示す各過程をコンピュータに実行させるための入力信号推定プログラムをＣＤ−ＲＯＭ、磁気ディスクなどの記録媒体から、あるいは通信回線を介してプログラムメモリ３４にダウンロードして、中央演算処理部（ＣＰＵ）により、プログラムメモリ３４に格納されているプログラムを実行させればよい。
入力信号ｕ（ｎ）の生成方法が決れば、後段フィルタ伝達関数ｗ_F は一定であるから、最初に１回計算して、図１では後段フィルタ２２にセットし、図３では記憶しておけばよい。前段フィルタ伝達関数ｗ ₁ ，ｗ ₂ は信号源、センサの配置の変更、その系１１内の空間の変更があれば、その都度、計算する必要がある。従って、前段フィルタ伝達関数は適当な間隔で計算して、適応的に変化させるとよい。
【００２４】
上述では１入力２出力線形系を対象としているが、１入力Ｎ出力線形系（Ｎは２以上の整数）を対象とする際には、センサ数、遅延器数、及び遅延器の出力を処理するフィルタ数をそれぞれＮへ変更し、これらフィルタの伝達関数の次数Ｌと系の伝達関数の次数ｍとの関係をＬ＞（ｍ−１）／Ｎとすれば良く、上記フィルタ演算部２３の計算原理や、減算器１７の出力を処理する後段フィルタ２２の処理、及び効果は変わらない。また入力信号ｕ（ｎ）としてはＡＲ生成過程に基づき生成されたものに限らない。またこの発明は室内の音声収録に限らず、移動無線におけるマルチパスによる影響の除去にも適用でき、その場合はセンサとしてアンテナが用いられ、そのアンテナ出力信号はベースバンド信号に変換されて、前述した場合と同様に処理される。その他同様の問題に対しこの発明を適用することができる。
【００２５】
【発明の効果】
以上説明した様に、この発明によれば、ＡＲ過程より生成される信号が１入力多出力線形系へ入力される際にも、系の出力端に備える複数センサの各出力を遅延させ、これら遅延信号それぞれをフィルタ処理した後加算し、第１センサの出力より減算して得られる信号（上記ＡＲ過程を駆動する線形予測誤差）を、ＡＲ過程の分母多項式に同等な多項式の逆元に相当する伝達関数によりフィルタで処理する事より、系への入力信号を高精度に推定する事が出来る。
【図面の簡単な説明】
【図１】この発明の一実施形態に係る入力信号推定装置の機能構成を示すブロック図。
【図２】図１中のフィルタ計算部２３の具体的機能構成を示すブロック図。
【図３】この発明の方法の実施形態を示す流れ図。
【図４】コンピュータの機能構成を簡略に示すブロック図。
【図５】従来の入力信号推定装置の機能構成を示すブロック図。[0001]
BACKGROUND OF THE INVENTION
This invention is caused by, for example, sound collection technology for reducing reverberation distortion, which is the main cause of quality deterioration of recorded sound, and for preventing intelligibility deterioration, such as buildings and highways when receiving mobile radios, for example, when recording sound in a room It is used for reception technology that reduces multipath distortion and maintains good speech quality. In a single-input multi-output linear system, when the output signal is distorted due to the transfer function of the system, the transfer function-specific prior The present invention relates to a method and apparatus for blindly estimating an input signal with high accuracy without using information, an input signal estimation program, and a recording medium thereof.
[0002]
[Prior art]
As this type of prior art, the one disclosed in Patent Document 1 is known. This will be briefly described with reference to FIG. An input signal u (n) at time n is input to the input terminal 12 of the 1-input 2-output linear system 11, and this input signal is supplied to two sensors 13 ₁ provided at the output terminal of the 1-input 2-output linear system 11. 13 ₂ is detected as an electrical signal. Sensor 13 _1, 13 the output signals of the _two are respectively one sample delayed by the delay unit 14 _1, 14 _2, it is supplied to the filter 15 _1, 15 _2. The output signals filtered by the filters 15 ₁ and 15 ₂ are added to each other by the adder 16, and this added output is subtracted from the output signal of the sensor 13 _{1 by} the subtractor 17. The output signal y (n) of the subtractor 17 and the delayed output signals from the delay units 14 ₁ and 14 ₂ are supplied to the filter calculation unit 18, and the filter calculation unit 18 uses the mean square value of the output signal y (n). the calculated transfer function of the filter 15 _1, 15 ₂ so as to minimize, these transfer functions are set in the filter 15 _1, 15 _2.
[0003]
When the input signal u (n) at the input terminal 12 satisfies the following relationship, the output signal y (n) of the subtractor 17 matches the input signal u (n), and the input signal is estimated with high accuracy. That is, the following equation is established and the input signal u (n) that is not affected by the reverberant sound is obtained as the output signal y (n).
E {u (n) u (n−j)} = 1 (j = 0) (1a)
E {u (n) u (n−j)} = 0 (j ≠ 0) (1b)
However, E {}: Average [0004]
[Patent Document 1]
Japanese Patent No. 3341811 (issued on November 5, 2002)
[Non-Patent Document 1]
Kodama, Suda, Matrix theory for system control, Corona, 1978
[Problems to be solved by the invention]
This conventional technique can provide good results when the input signal is a random signal.
However, when the generation process of the input signal u (n) is an autoregressive (AR) process, the output signal y (n) of the subtractor 17 does not match the input signal u (n), and the input signal u (n n) is not estimated accurately. The reason for this is explained as follows. If the coefficients of the AR generation process of the input signal u (n) are a ₁ ,..., A _p , the input signal u (n) satisfies the following relationship.

However, a (z) = a ₁ z ⁻¹ + a ₂ z ⁻² +... + A _p z ^−p
Or
[Expression 1]

_{^{C = E {u n T u}} n-1} (E {u n-1 T u n-1}) -1 (2c)
Where e (n): linear prediction error, ^T : transpose
Transfer functions h ₁ (z) and h ₂ (z) (z: complex variable) between the input terminal 12 and the sensors 13 ₁ and 13 ₂ and the transfer functions w ₁ (z) of the filters 15 ₁ and 15 ₂ , W ₂ (z) are expressed as follows.

However, h _j (z) (j = 1, 2) does not have a common zero.
[0007]
At this time, the output signal y (n) of the subtractor 17 is expressed as follows.
[Expression 2]

The relationship between the order L of the filter 15 _j and the order m of the transfer function of the system 11 is determined as L ≧ m−1, and the root mean square value E {| y (n + 1) | ² of the output signal y (n) of the subtractor 17. } Is minimized, the vector w representing the transfer functions of the filters 15 ₁ and 15 ₂ is obtained as follows from the equation (2c).

Where s ^{2 is a} small positive number, I is a unit matrix
Here _{(H E {u n-1} T u n-1} H T + s 2 I) -1: a matrix ^{_{H E {u n-1 T}} u n-1} pseudoinverse s ² of H ^T sufficiently small By taking a positive number, the relationship of equation (5) is arranged as follows (Non-Patent Document 1, page 339).
[Equation 3]

However, A ^+: Moore-Penrose generalized inverse matrix of the matrix _{^{A, U n-1 T U}} n-1 = E
{ U _n-1 ^T u _n-1 }
By substituting equation (5) into equation (4) and using the relationship of equation (2b), the output signal y (n) can be obtained as follows.

[0009]
Thus, y (n) does not become u (n), that is, the input signal u (n) generated from the AR process is not estimated by the technique shown in Patent Document 1.
An object of the present invention is to use an input signal generated from an AR process only from an output signal of each sensor without using prior information related to a transfer function of the system in a one-input multi-output linear system such as room sound pickup. To provide a method and apparatus for estimating blindness with accuracy, a program thereof, and a recording medium.
[0010]
[Means for Solving the Problems]
According to the present invention, a delay of one sample is given to each of the output signals of the first, second,..., Nth sensors provided at the output end of the one-input multi-output linear system, and the N output signals and the above-mentioned N delay signals are used to calculate N preceding filter transfer functions and one succeeding filter transfer function, and the N delayed signals are filtered to the N preceding filter transfer functions, respectively. The processed signals are added, the added signal is subtracted from the output signal of the first sensor, and the subtracted signal is filtered by the post-filter transfer function.
[0011]
The filter transfer function can be calculated as follows, for example. Obtain an autocorrelation matrix of N delayed signals, further obtain a pseudo inverse matrix of the autocorrelation matrix, obtain a cross correlation matrix between the N delayed signals and the output signals of the N sensors, A product matrix of the pseudo inverse matrix of the autocorrelation matrix and the cross-correlation matrix is obtained. The N components in the first row of the product matrix are defined as the N preceding-stage filter transfer functions. Further, a characteristic equation of the product matrix is obtained, and from this, a polynomial ± {1-a (z)} equivalent to the denominator polynomial in the AR generation process is obtained, and an inverse element ± 1 / (1-a (z) of this polynomial is obtained. )) To obtain the subsequent filter transfer function. The subtracted signal h _{1: 0} e (n) is processed by the post-stage filter 22, and as a result of the processing, e (n) is obtained from the relationship of the equation (2a) as shown in the following equation (8). A signal h _{1: 0} u (n) proportional to n) is obtained.
± 1 / (1-a (z)) h _{1: 0} e (n) = ± h _{1: 0} u (n) (8)
That is, the input signal to the system is estimated with high accuracy.
[0012]
DETAILED DESCRIPTION OF THE INVENTION
Embodiments of the present invention will be described below with reference to the drawings. FIG. 1 shows a functional configuration of an embodiment of the input signal estimation apparatus of the present invention. An input signal u (n) at time n input to the input terminal 12 of the one-input two-output linear system 11 is detected as an electrical signal by the _two sensors 13 ₁ and 13 ₂ at the output terminal, and each of these sensor output signals Is delayed by one sample by the delay units 14 ₁ and 14 ₂ , these delayed signals are filtered by the pre-filters 21 ₁ and 21 ₂ , and these filter output signals are added by the adder 16, and this added signal is the first sensor 13 _1. Is subtracted by the subtracter 17, and the subtracted signal is filtered by the post-stage filter 22. Based on the output signals of the sensors 13 ₁ and 13 ₂ and the delayed signals of the delay devices 14 ₁ and 14 ₂ , the filter calculation unit 23 calculates the pre-stage filter transfer functions w ₁ and w ₂ and the post-stage filter transfer function, Functions are set in the pre-filters 21 ₁ and 21 ₂ and the post-filter 22 respectively.
[0013]
Filter calculation An example of the functional configuration of the filter calculation unit 23 is shown in FIG.
An autocorrelation matrix R _A is calculated by the autocorrelation calculation unit 24 from the output signals of the delay devices 14 ₁ and 14 ₂ , and this pseudo inverse matrix ( R _A + s ² I ) ⁻¹ is calculated by the pseudo inverse matrix calculation unit 25. The A cross-correlation matrix R _c between the output signals of the delay devices 14 ₁ and 14 _{2 and} the output signals of the sensors 13 ₁ and 13 ₂ is calculated by the cross-correlation calculating unit 26.
[0014]
_A multiplication unit 27 calculates a product matrix W = R _c ( R _A + s ² I ) ^{−1 of the} pseudo inverse matrix R _A ⁻ and the cross-correlation matrix R _c . Two components in the first row of the product matrix W are set in the pre-filters 21 ₁ and 21 ₂ as the pre-filter transfer functions w ₁ and w ₂ . Further, a characteristic polynomial of the product matrix W is calculated by the polynomial calculation unit 28, and an inverse element of a polynomial equivalent to the denominator polynomial in the AR generation process in the characteristic polynomial is calculated by the inverse element calculation unit 29, and the post-stage filter transfer function w _F is set in the subsequent filter 22.
[0015]
Next, the processing in the filter calculation unit 23 will be described in detail. The input signal u (n) to the system 11 and the AR process for generating u (n), the transfer function h _j (z) of the system 11, and the transfer function w _j (of the pre-stage filter 21 _j (j = 1, 2) The relationship between the order m of z) and h _j (z) and the order L of w _j (z) is the same as in the case of the above-described equations (2) to (5). The autocorrelation matrix R _A obtained from the output signals of the delay devices 14 ₁ and 14 ₂ is obtained by calculation of the autocorrelation calculation unit 24 by calculation of the following equation (9).
R _A = HE {u _n-1 ^T u _n-1 } H ^T (9)
When the pseudo inverse matrix R _A ^{− is obtained} by the inverse matrix calculation unit 25, R _A becomes (R _A + s ² I) ⁻¹ , and when the formula (9) is substituted into this pseudo inverse matrix, the following formula (10) is obtained. Thus, it becomes the same as the part of the inverse matrix in equation (5).
R _A ⁻ = (R _A + s ² I) ⁻¹ = (HE {un _n−1 ^T u _n−1 } H ^T + s ² I) ⁻¹ (10)
[0016]
The cross-correlation matrix R _C is obtained by the cross-correlation calculating unit 26 from the outputs of the delay devices 14 ₁ and 14 ₂ and the sensors 13 ₁ and 13 ₂ by the calculation of the following equation (11).
_{R C = H E {u n} -1 T u n} H T (11)
Therefore, the product matrix W calculated by the multiplication unit 27 is expressed by the following equation (12).

By substituting the relationship of equation (2c) into equation (12) and making s ² a sufficiently small positive number, the relationship of equation (12) can be arranged as follows, as in equation (6).
[Expression 4]

[0017]
Samples each component of the first row W _1r of the product matrix W but are set as respective transfer functions upstream filter 21 _1, 21 ₂ above, the first row of the product matrix W is discussed in the section of the prior art It is equal to the transfer function vector w _{1 of} the filters 15 ₁ and 15 ₂ for filtering the delayed signal, that is, the expression (6).

Therefore, delay unit 14 _1, 14 front filter 21 ₁ to output delay signals of the _2, 21 ₂ filters respectively, these front filter 21 _1, 21 a sum signal of the _second output signal the first sensor 13 the _first output signal When the generation process of the input signal u (n) to the system 11 is an AR process, the output signal y (n) of the subtractor 17 obtained by further subtraction is given by the expression (2a) as shown in the expression (7). The signal h _{1: 0} e (n) proportional to the linear prediction error e (n) that drives the AR process 1 / (1-a (z)) that generates the input signal u (n) indicated by
[0018]
Incidentally, the non-zero eigenvalues λ _i { W } of the product matrix W satisfy the following relationship (Non-Patent Document 1, pages 166 to 167).

The characteristic polynomial c _c (z) of the matrix C is obtained from the equation (2b) as follows.

From a relationship of equation (15), since the product matrix W elements of one row to the matrix U is 2 (L + 1) -p by more, W of the characteristic polynomial c _w (z) can be expressed as follows .

[0019]
By multiplying the characteristic polynomial c _w (z) by z ^{−2 (L + 1)} , a polynomial ± {1-a (z)} equivalent to the denominator polynomial in the AR generation process shown in Expression (2a) is obtained. Can do. Actually, the polynomial calculation unit 28 calculates the characteristic polynomial c _w (z) of the product matrix W, and takes out the polynomial ± {1-a (z)} equivalent to the denominator polynomial in the AR generation process. Inverse element ± 1 / (1-a (z)) is calculated by the inverse element calculation unit 29.
This inverse element ± 1 / (1-a (z)) = w _F is set as a transfer function w _{F in the} subsequent filter. The output signal h _{1: 0} e (n) of the subtractor 17 is filtered by the post filter 22. By this processing, equation (8) is calculated, and a signal ± h _{1: 0} u (n) proportional to the input signal u (n) of the system 11 is obtained. That is, the input signal to the system 11 is estimated with high accuracy.
[0020]
Next, an embodiment in which the method of the present invention is applied to the one-input two-output linear system 11 shown in FIG. 1 will be briefly described with reference to FIG. It is assumed that each procedure of this method is executed by a computer schematically shown in FIG.
Although the output signals of the sensors 13 ₁ and 13 ₂ are not shown in the figure, they are usually sampled at regular intervals and once taken into the storage unit 32 through the interface 31 as a digital signal sequence (S1). Each sensor output signal is given a delay of one sample (S2), and the pre-stage filter transfer functions w ₁ and w ₂ and the post-stage filter transfer function w _F are calculated based on these delay signals and the sensor output signal ( S3).
[0021]
That computes the autocorrelation matrix of the process as well as both the delay signal _{R A = H E {u n} -1 T u n-1} H T shown in FIG. 2, for example (S3-1), further the pseudo-inverse matrix R _A ⁻ = ( R _A + s ² I ) ⁻¹ is calculated and temporarily stored (S3-2). Also by calculating the cross-correlation matrix _{R C = H E {u n} -1 T u n} H T between the both sensor output signal and two delayed signals (S3-3), the cross-correlation matrix R _C and the pseudo inverse matrix R _A ^- calculating a product matrix W with (S3-4). Each component of the first row of the product matrix W is stored as the pre-stage filter transfer functions w ₁ and w ₂ (S3-5).
[0022]
Further, the characteristic polynomial c _w (z) of the product matrix W is calculated (S3-6), and the inverse element ± 1 / (1-a (z)) of the polynomial equivalent to the denominator polynomial of the AR generation process in the characteristic polynomial Is stored as a post-stage filter transfer function w _F (S3-7).
After obtaining the pre-stage filter transfer function and the post-stage filter transfer function in this way, each delayed signal is filtered using the corresponding pre-stage filter transfer functions w ₁ and w ₂ (S4), and these pre-filter processed signals are obtained. (S5), and subtracts the added signal from the first sensor output signal (S6), filters the subtracted signal with the post-stage filter transfer function w _F (S7), and outputs the result through the output unit 33. Output (S8).
[0023]
An input signal estimation program for causing the computer to execute each process shown in FIG. 3 is downloaded to a program memory 34 from a recording medium such as a CD-ROM or a magnetic disk or via a communication line, and a central processing unit ( The CPU may execute a program stored in the program memory 34.
If the generation method of the input signal u (n) is determined, the post-stage filter transfer function w _F is constant, so it is calculated once, set in the post-filter 22 in FIG. 1, and stored in FIG. Just keep it. The pre-stage filter transfer functions w ₁ and w ₂ need to be calculated each time there is a change in the signal source and sensor arrangement, and a change in the space in the system 11. Therefore, the pre-stage filter transfer function is preferably calculated at appropriate intervals and adaptively changed.
[0024]
In the above, a 1-input 2-output linear system is targeted. However, when a 1-input N-output linear system (N is an integer of 2 or more) is targeted, the number of sensors, the number of delay units, and the output of the delay unit are processed. The number of filters to be changed is changed to N, and the relationship between the order L of the transfer function of these filters and the order m of the transfer function of the system may be L > (m−1) / N. The calculation principle, the processing of the post-filter 22 that processes the output of the subtractor 17, and the effect remain unchanged. The input signal u (n) is not limited to that generated based on the AR generation process. The present invention is not limited to indoor audio recording, but can also be applied to the removal of the effects of multipath in mobile radio. In this case, an antenna is used as a sensor, and the antenna output signal is converted into a baseband signal. It is processed in the same way. The present invention can be applied to other similar problems.
[0025]
【The invention's effect】
As described above, according to the present invention, even when a signal generated from the AR process is input to a one-input multi-output linear system, the outputs of a plurality of sensors provided at the output end of the system are delayed, Each delayed signal is filtered and then added, and the signal obtained by subtracting from the output of the first sensor (the linear prediction error that drives the AR process) is equivalent to the inverse element of a polynomial equivalent to the denominator polynomial of the AR process. The input signal to the system can be estimated with high accuracy by processing with a filter according to the transfer function.
[Brief description of the drawings]
FIG. 1 is a block diagram showing a functional configuration of an input signal estimation apparatus according to an embodiment of the present invention.
FIG. 2 is a block diagram showing a specific functional configuration of a filter calculation unit 23 in FIG.
FIG. 3 is a flowchart illustrating an embodiment of the method of the present invention.
FIG. 4 is a block diagram schematically showing a functional configuration of a computer.
FIG. 5 is a block diagram showing a functional configuration of a conventional input signal estimation apparatus.

Claims (4)

A process of delaying the output signals from the first, second,..., Nth sensors provided at the integer N number of output terminals of two or more integers of the one-input multi-output linear system by one sample,
A filter calculation process for calculating and storing N pre-stage filter transfer functions and one post-stage filter transfer function based on the N output signals and the delayed N signals;
A pre-filtering process for filtering the delayed first, second,..., N signals with the first, second,.
Adding each filtered N signal;
Subtracting the added signal from the output signal from the first sensor;
The subtracted signal possess a subsequent filtering step of filtering by the subsequent filter transfer function,
The filter calculation process is as follows:
Calculating the autocorrelation matrix of the delayed first, second,..., Nth signals;
Calculating a pseudo inverse of the autocorrelation matrix;
Calculating the cross-correlation matrix between the output signals of the first, second,..., Nth sensors and the delayed first, second,.
Calculating a product matrix of the pseudo inverse matrix of the autocorrelation matrix and the cross-correlation matrix and obtaining N components of the first row as the first, second,.
Calculating the characteristic polynomial of the product matrix,
Multiplying the characteristic polynomial by z ^{−N (L + 1)} (L is the order of the transfer function of the preceding filter) to obtain the transfer function in the latter filter process;
An input signal estimation method comprising: 1st, 2nd, ..., Nth sensors provided at the output end of the 1-input multi-output linear system (N is an integer of 2 or more),
1st, 2nd,..., Nth delay devices that give a delay of 1 sample to each of the output signals of these sensors;
A first, second,..., Nth pre-filter that processes each of the output signals of these delay devices;
An adder for adding the output signals of these filters;
A subtractor for subtracting the output signal of the adder from the output signal of the first sensor;
A post filter for processing the output signal of the subtractor,
Each delayed output signal of the output signal and the respective delay units of the sensors is input, the first, second, ..., and the N filter and the transfer function of the post-stage filter calculated, respectively, which filter transfer function A filter calculator for setting each to the corresponding one of the filters ,
The filter coefficient calculator is
An autocorrelation calculation unit for calculating an autocorrelation matrix of the delayed output signals of the first, second,.
An inverse matrix calculation unit for calculating a pseudo inverse matrix of the autocorrelation matrix;
A cross-correlation calculator that calculates a cross-correlation matrix between the output signals of the first, second,..., Nth sensors and the delayed output signals of the first, second,.
A product matrix of the pseudo inverse matrix and the cross-correlation matrix is calculated, and N components in the first row of the product matrix are used as filter transfer functions to the first, second,. A multiplier to set;
A polynomial calculator for calculating the characteristic polynomial of the product matrix;
And an inverse element calculation unit that multiplies the characteristic polynomial by z ^{−N (L + 1)} (L is the order of the transfer function of the preceding filter) and sets the resultant polynomial as a filter transfer function to the subsequent filter. Signal estimation device. An input signal estimation program for causing a computer to execute each step of the input signal estimation method according to claim 1 . A computer-readable recording medium in which the input signal estimation program according to claim 3 is recorded.

Images