(57)【特許請求の範囲】
【請求項１】
フィルタリング部において参照信号をフィルタ係数を適用してフィルタリングすることで推定信号を求め、差分器で上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新部において上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小になるように上記フィルタ係数を更新する伝達系推定装置において、
上記フィルタ係数更新部が、第一周波数領域変換部において上記参照信号を周波数領域の周波数領域参照信号に変換し、参照信号バッファにおいて現時刻の上記周波数領域参照信号を蓄積することで参照信号バッファの内容を更新し、乗算器で上記参照信号バッファの内容とフィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量の内積を周波数毎に取りその結果として周波数領域推定誤差信号を生成し、
第二時間領域変換部において上記周波数領域推定誤差信号を時間領域の選択前推定誤差信号に変換し、第二要素選択部において上記選択前推定信号の一部を選択し推定誤差信号とし、差分器で上記推定誤差信号と上記誤差信号の差である修正量更新用誤差を求め、
第二周波数領域変換部において上記修正量更新用誤差を周波数領域信号の周波数領域修正量更新用誤差に変換し、更新量推定部において上記周波数領域修正量更新用誤差と上記参照信号バッファの内容から、周波数領域修正量更新量を求め、
第一平滑部において上記周波数領域非拘束フィルタ係数修正量に上記周波数領域修正量更新量を重み付け加算し、
第一拘束部分決定部において第一拘束分割番号を決定し、
第二拘束番号決定部で上記第一拘束分割番号以外の任意の番号の中から第二拘束分割番号を決定し、
拘束部において、
上記第一平滑部における重み付け加算結果のうち上記第一拘束分割番号に相当しない拘束分割番号に相当する上記周波数領域非拘束フィルタ係数修正量のうち、
時間領域に変換した際に上記第二拘束分割番号の部分に相当する上記周波数領域非拘束フィルタ係数修正量の時間軸方向の前半部分又は後半部分の値が０になるように修正し、
上記第二拘束分割番号に相当しない上記周波数領域非拘束フィルタ係数修正量についてはそのまま上記フィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量に上書きし、
第一時間領域変換部において上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域信号の非拘束フィルタ係数修正量に変換し、第一要素選択部において上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、第二平滑部において上記フィルタ係数に上記フィルタ係数修正量を加算し、
上記第二平滑部における加算結果をFGフィルタ係数バッファの内容に上書きすると共に、
再初期値計算部において再初期値を求め、この再初期値を上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量に設定することを特徴とする伝達系推定装置。
【請求項２】
フィルタリング部において参照信号をフィルタ係数を適用してフィルタリングすることで推定信号を求め、差分器で上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新部において上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小になるように上記フィルタ係数を更新する伝達系推定装置において、
上記フィルタ係数更新部が、第一周波数領域変換部において上記参照信号を周波数領域の周波数領域参照信号に変換し、参照信号バッファにおいて現時刻の上記周波数領域参照信号を蓄積することで参照信号バッファの内容を更新し、乗算器で上記参照信号バッファの内容とフィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量の内積を周波数毎に取りその結果として周波数領域推定誤差信号を生成し、第二周波数領域変換部において上記誤差を周波数領域信号の周波数領域誤差信号に変換し、周波数領域差分器において上記周波数領域誤差信号と上記周波数領域推定誤差信号との差をとることで周波数領域修正量更新用近似誤差を求め、更新量計算部において上記周波数領域修正量更新用近似誤差と上記参照信号バッファの内容から周波数領域修正量更新量を求め、第一平滑部において上記周波数領域非拘束フィルタ係数修正量に上記周波数領域修正量更新量を重み付け加算し、第一拘束部分決定部において第一拘束分割番号を決定し、
第二拘束番号決定部で上記第一拘束分割番号以外の任意の番号の中から第二拘束分割番号を決定し、
拘束部において、
上記第一平滑部における重み付け加算結果のうち上記第一拘束分割番号に相当しない拘束分割番号に相当する上記周波数領域非拘束フィルタ係数修正量のうち、
時間領域に変換した際に、上記第二拘束分割番号の部分に相当する上記周波数領域非拘束フィルタ係数修正量の時間軸方向の前半部分又は後半部分の値が０になるように修正し、
上記第二拘束分割番号に相当しない上記周波数領域非拘束フィルタ係数修正量についてはそのまま上記フィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量に上書きし、
第一時間領域変換部において上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域信号の非拘束フィルタ係数修正量に変換し、第一要素選択部において上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、第二平滑部において上記フィルタ係数に上記フィルタ係数修正量を加算し、
上記第二平滑部における加算結果をFGフィルタ係数バッファの内容に上書きすると共に、
再初期値計算部において再初期値を求め、この再初期値を上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量に設定することを特徴とする伝達系推定装置。
【請求項３】
フィルタリング部において参照信号をフィルタ係数を適用してフィルタリングすることで推定信号を求め、差分器で上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新部において上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小になるように上記フィルタ係数を更新する伝達系推定装置において、
上記フィルタ係数更新部が、第一周波数領域変換部において上記参照信号を周波数領域の周波数領域参照信号に変換し、参照信号バッファにおいて現時刻の上記周波数領域参照信号を蓄積することで参照信号バッファの内容を更新し、乗算器で上記参照信号バッファの内容とフィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量の内積を周波数毎に取りその結果として周波数領域推定誤差信号を生成し、第二周波数領域変換部において上記誤差を周波数領域信号の周波数領域誤差信号に変換し、周波数領域差分器において上記周波数領域誤差信号と上記周波数領域推定誤差信号との差をとることで周波数領域修正量更新用近似誤差を求め、更新量計算部において上記周波数領域修正量更新用近似誤差と上記参照信号バッファの内容から周波数領域修正量更新量を求め、第一平滑部において上記周波数領域非拘束フィルタ係数修正量に上記周波数領域修正量更新量を重み付け加算し、第一拘束部分決定部において第一拘束分割番号を決定し、第二拘束番号決定部で上記第一拘束分割番号以外の任意の番号の中から第二拘束分割番号を決定し、
拘束部において、
上記第一平滑部における重み付け加算結果のうち上記第一拘束分割番号に相当しない拘束分割番号に相当する上記周波数領域非拘束フィルタ係数修正量のうち、
時間領域に変換した際に上記第二拘束分割番号の部分に相当する上記周波数領域非拘束フィルタ係数修正量の時間軸方向の前半部分又は後半部分の値が０になるように修正し、
上記第二拘束分割番号に相当しない上記周波数領域非拘束フィルタ係数修正量についてはそのまま上記フィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量に上書きし、
第一時間領域変換部において上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域信号の非拘束フィルタ係数修正量に変換し、第一要素選択部において上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、第二平滑部において上記フィルタ係数に上記フィルタ係数修正量を加算し、
上記第二平滑部における加算結果をFGフィルタ係数バッファの内容に上書きすると共に、
再初期値計算部において再初期値を求め、この再初期値を上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量に設定することを特徴とする伝達系推定装置。
【請求項４】
請求項２又は３の何れかに記載の伝達系推定装置において、上記フィルタ係数更新部の上記乗算器の出力を上記周波数領域差分器に入力することに替え、上記フィルタ係数更新部の乗算器の出力を近似拘束部に入力し、
上記周波数領域推定誤差信号を時間領域に変換した際に、時間軸方向の前半部分又は後半部分の上記周波数領域推定誤差信号の値が０になるようにする拘束処理を近似する行列と、上記周波数領域推定誤差信号との積により、近似拘束周波数領域推定誤差信号を求め、
上記近似拘束部の出力である上記近似拘束周波数領域推定誤差信号を上記周波数領域差分器に入力することを特徴とする伝達系推定装置。
【請求項５】
条件C1で参照信号と所望信号がNサンプル得られるまで待ち、参照信号と所望信号がNサンプル得られると、
ステップS1において参照信号を周波数領域に変換することで周波数領域参照信号を求め、
ステップS2において上記周波数領域参照信号を参照信号バッファに蓄積し、
ステップS3においてＦＧフィルタ係数バッファにあるフィルタ係数あるいは周波数領域フィルタ係数と参照信号あるいは参照信号バッファにある周波数領域参照信号から推定信号を求め、
ステップS4において上記推定信号と所望信号の差をとることで誤差信号を求め、
ステップS10においてフィルタ係数修正量バッファにある周波数領域非拘束フィルタ係数修正量と参照信号バッファの内容との内積を周波数毎にとることにより周波数領域推定誤差信号を求め、
ステップS11において上記周波数領域推定誤差信号を時間領域に変換した信号の前半を０にすることで推定誤差信号を求め、
ステップS12において上記推定誤差信号と誤差信号の差により修正用誤差信号を求め、
ステップS13において、上記修正用誤差信号を周波数領域に変換し周波数領域修正用誤差を求め、
ステップS14において上記周波数領域修正用誤差と参照信号バッファの内容から周波数領域修正量更新量を求め、
ステップS15において上記周波数領域修正量更新量をフィルタ係数修正量バッファの内容に重み付け加算し周波数領域非拘束フィルタ係数修正量を得て上記周波数領域非拘束フィルタ係数修正量を上記フィルタ係数修正量バッファに保存し、
ステップS16において第二拘束分割番号ｑ’を決定し、上記周波数領域非拘束フィルタ係数修正量のうち時間領域に変換した際に上記第二拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量の時間軸方向の前半部分又は後半部分の値が０になるようにした後、上記フィルタ係数修正量をバッファに保存し、
条件C1に戻り、また、上記条件C1とステップS1、S2、S3、S4、S10、S11、S12、S13、S14、S15、S16と同期あるいは非同期に条件C2で更新間隔が経過するのを待ち、
更新間隔が経過したなら、ステップS17において第一拘束分割番号を決定し上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域の非拘束フィルタ係数修正量に変換し、上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、
ステップS18において上記フィルタ係数修正量を上記FGフィルタ係数バッファの内容に加算するか、あるいは、上記フィルタ係数修正量を周波数領域に変換し周波数領域フィルタ係数修正量を得、上記周波数領域フィルタ係数修正量を上記FGフィルタ係数バッファの内容に加算し、
ステップS19において更新したフィルタ係数あるいは周波数領域フィルタ係数の上記第一拘束分割番号に対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量を再初期値にし、
条件C2に戻ることを特徴とする伝達系推定方法。
【請求項６】
条件C1で参照信号と所望信号がNサンプル得られるまで待ち、参照信号と所望信号がNサンプル得られると、
ステップS1において参照信号を周波数領域に変換することで周波数領域参照信号を求め、ステップS2において上記周波数領域参照信号を参照信号バッファに蓄積し、
ステップS3において時間領域で参照信号とフィルタ係数をフィルタリングするか、あるいは、それに等価な処理を周波数領域で行うことで推定信号を求め、
ステップS4において上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新ステップにおいて上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小となるように上記フィルタ係数を更新し、条件C1に戻る伝達系推定方法において、フィルタ係数更新ステップが、
ステップS10において上記参照信号バッファの内容と周波数領域非拘束フィルタ係数修正量の内積を周波数毎にとることで周波数領域推定誤差信号を得、
ステップS5において上記誤差信号を周波数領域に変換することで周波数領域誤差信号を求め、
ステップS30において上記周波数領域誤差信号から上記周波数領域推定誤差信号を減じることで周波数領域修正量更新用近似誤差を得、
ステップS31において上記周波数領域修正量更新用近似誤差と上記参照信号バッファの内容から周波数領域修正量更新量を得、
ステップS15において上記周波数領域修正量更新量をフィルタ係数修正量バッファの内容に重み付け加算し上記周波数領域非拘束フィルタ係数修正量を更新し、上記周波数領域非拘束フィルタ係数修正量を上記フィルタ係数修正量バッファに保存し、
ステップS16において第二拘束分割番号ｑ’を決定し、上記周波数領域非拘束フィルタ係数修正量のうち、時間領域に変換した際に上記第二拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量の時間軸方向の前半部分又は後半部分の値が０になるようにした後、上記フィルタ係数修正量バッファに保存し、条件C1に戻り、また、上記条件C1とステップS1、S2、S3、S4、S10、S5、S30、S31、S15と同期、あるいは非同期で、条件C2において更新間隔が経過するのを待ち、更新間隔が経過したなら、
ステップS17において第一拘束分割番号を決定し上記第一拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域の非拘束フィルタ係数修正量に変換し、上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、
ステップS18において上記フィルタ係数修正量をFGフィルタ係数バッファの内容に加算するか、あるいは、上記フィルタ係数修正量を周波数領域に変換し周波数領域フィルタ係数修正量を得、上記周波数領域フィルタ係数修正量を上記FGフィルタ係数バッファの内容に加算し、ステップS19において更新したフィルタ係数あるいは周波数領域フィルタ係数の上記第一拘束分割番号に対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量を微少値にし、条件C2に戻る、ことを特徴とする伝達系推定方法。
【請求項７】
請求項６記載の伝達系推定方法において、ステップS5の直後で、ステップS40において上記周波数領域推定誤差信号を近似拘束処理し、近似拘束周波数領域誤差信号を得、ステップS30をステップS41に替えステップS41において上記周波数領域誤差信号から上記近似拘束周波数領域推定誤差信号を減じることで周波数領域修正量更新用近似誤差を得る、ことを特徴とする伝達系推定方法。
【請求項８】
フィルタリング部において参照信号をフィルタ係数を適用してフィルタリングすることで推定信号を求め、差分器で上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新部において上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小になるように上記フィルタ係数を更新する伝達系推定装置において、
上記フィルタ係数更新部が、第一周波数領域変換部において上記参照信号を周波数領域の周波数領域参照信号に変換し、参照信号バッファにおいて現時刻の上記周波数領域参照信号を蓄積することで参照信号バッファの内容を更新し、乗算器で上記参照信号バッファの内容とフィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量の内積を周波数毎に取りその結果として周波数領域推定誤差信号を生成し、
第二時間領域変換部において上記周波数領域推定誤差信号を時間領域の選択前推定誤差信号に変換し、
第二要素選択部において上記選択前推定誤差信号の一部を選択し推定誤差信号とし、
差分器で上記推定誤差信号と上記誤差信号の差である修正量更新用誤差を求め、
第二周波数領域変換部において上記修正量更新用誤差を周波数領域信号の周波数領域修正量更新用誤差に変換し、
更新量推定部において上記周波数領域修正量更新用誤差と上記参照信号バッファの内容から、周波数領域修正量更新量を求め、
第一平滑部において上記周波数領域非拘束フィルタ係数修正量に上記周波数領域修正量更新量を重み付け加算し、拘束部分決定部において拘束分割番号を決定し、
上記第一平滑部における重み付け加算結果のうち上記拘束分割番号に相当しない上記周波数領域非拘束フィルタ係数修正量を上記フィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量に上書きし、第一時間領域変換部において上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域信号の非拘束フィルタ係数修正量に変換し、第一要素選択部において上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、
第二平滑部において上記フィルタ係数に上記フィルタ係数修正量を加算し、上記第二平滑部における加算結果を上記ＦＧフィルタ係数バッファの内容に上書きすると共に、再初期値計算部において再初期値を求め、この再初期値を上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量に設定することを特徴とする伝達系推定装置。
【請求項９】
フィルタリング部において参照信号をフィルタ係数を適用してフィルタリングすることで推定信号を求め、差分器で上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新部において上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小になるように上記フィルタ係数を更新する伝達系推定装置において、
上記フィルタ係数更新部が、第一周波数領域変換部において上記参照信号を周波数領域の周波数領域参照信号に変換し、参照信号バッファにおいて現時刻の上記周波数領域参照信号を蓄積することで参照信号バッファの内容を更新し、乗算器で上記参照信号バッファの内容とフィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量の内積を周波数毎に取りその結果として周波数領域推定誤差信号を生成し、第二周波数領域変換部において上記誤差信号を周波数領域信号の周波数領域誤差信号に変換し、周波数領域差分器において上記周波数領域誤差信号と上記周波数領域推定誤差信号との差をとることで周波数領域修正量更新用近似誤差を求め、更新量計算部において上記周波数領域修正量更新用近似誤差と上記参照信号バッファの内容から周波数領域修正量更新量を求め、
第一平滑部において上記周波数領域非拘束フィルタ係数修正量に上記周波数領域修正量更新量を重み付け加算し、
拘束部分決定部において拘束分割番号を決定し、上記第一平滑部における重み付け加算結果のうち上記拘束分割番号に相当しない上記周波数領域非拘束フィルタ係数修正量の一部あるいは全部を上記フィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量に上書きし、
第一時間領域変換部において上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域信号の非拘束フィルタ係数修正量に変換し、第一要素選択部において上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、
第二平滑部において上記拘束分割番号に対応する上記フィルタ係数に上記フィルタ係数修正量を加算し、
上記第二平滑部における加算結果を上記FGフィルタ係数バッファの内容に上書きすると共に、再初期値計算部において再初期値を求め、この再初期値を上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量に設定することを特徴とする伝達系推定装置。
【請求項１０】
フィルタリング部において参照信号をフィルタ係数を適用してフィルタリングすることで推定信号を求め、差分器で上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新部において上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小になるように上記フィルタ係数を更新する伝達系推定装置において、
上記フィルタ係数更新部が、第一周波数領域変換部において上記参照信号を周波数領域の周波数領域参照信号に変換し、参照信号バッファにおいて現時刻の上記周波数領域参照信号を蓄積することで参照信号バッファの内容を更新し、乗算器で上記参照信号バッファの内容とフィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量の内積を周波数毎に取りその結果として周波数領域推定誤差信号を生成し、第二周波数領域変換部において上記誤差信号を周波数領域信号の周波数領域誤差信号に変換し、周波数領域差分器において上記周波数領域誤差信号と上記周波数領域推定誤差信号との差をとることで周波数領域修正量更新用近似誤差を求め、更新量計算部において上記周波数領域修正量更新用近似誤差と上記参照信号バッファの内容から周波数領域修正量更新量を求め、
第一平滑部において上記周波数領域非拘束フィルタ係数修正量に上記周波数領域修正量更新量を重み付け加算し、拘束部分決定部において拘束分割番号を決定し、上記第一平滑部における重み付け加算結果のうち上記拘束分割番号に相当しない上記周波数領域非拘束フィルタ係数修正量の一部あるいは全部を上記フィルタ係数修正量バッファの内容である周波数領域非拘束フィルタ係数修正量に上書きし、第一時間領域変換部において上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域信号の非拘束フィルタ係数修正量に変換し、第一要素選択部において上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、第三周波数領域変換部において上記フィルタ係数修正量を周波数領域の信号である周波数領域フィルタ係数修正量に変換し、第二平滑部において上記拘束分割番号に対応する上記フィルタ係数に上記フィルタ係数修正量を加算し、
上記第二平滑部における加算結果を上記ＦＧフィルタ係数バッファの内容に上書きすると共に、再初期値計算部において再初期値を求め、この再初期値を上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量に設定することを特徴とする伝達系推定装置。
【請求項１１】
請求項９又は１０の何れかに記載の伝達系推定装置において、上記フィルタ係数更新部の上記乗算器の出力を上記周波数領域差分器に入力することに替え、上記フィルタ係数更新部の乗算器の出力を近似拘束部に入力し、
上記周波数領域推定誤差信号を時間領域に変換した際に、時間軸方向の前半部分又は後半部分の上記周波数領域推定誤差信号の値が０になるようする拘束処理を近似する行列と、上記周波数領域推定誤差信号との積により、近似拘束周波数領域推定誤差信号を求め、
上記近似拘束部の出力である上記近似拘束周波数領域推定誤差信号を上記周波数領域差分器に入力することを特徴とする伝達系推定装置。
【請求項１２】
条件C1で参照信号と所望信号がNサンプル得られるまで待ち、参照信号と所望信号がNサンプル得られると、ステップS1において参照信号を周波数領域に変換することで周波数領域参照信号を求め、ステップS2において上記周波数領域参照信号を参照信号バッファに蓄積し、ステップS3においてＦＧフィルタ係数バッファにあるフィルタ係数あるいは周波数領域フィルタ係数と参照信号あるいは参照信号バッファにある周波数領域参照信号から推定信号を求め、
ステップS4において上記推定信号と所望信号の差をとることで誤差信号を求め、
ステップS10においてフィルタ係数修正量バッファにある周波数領域非拘束フィルタ係数修正量と参照信号バッファの内容との内積により周波数領域推定誤差信号を求め、
ステップS11において上記周波数領域推定誤差信号を時間領域に変換し、その前半を０にして推定誤差信号を求め、
ステップS12において上記推定誤差信号と誤差信号の差により修正用誤差信号を求め、
ステップS13において、上記修正用誤差信号を周波数領域に変換し周波数領域修正用誤差を求め、
ステップS14において上記周波数領域修正用誤差と参照信号バッファの内容から周波数領域修正量更新量を求め、
ステップS15において上記周波数領域修正量更新量をフィルタ係数修正量バッファの内容に重み付け加算し上記周波数領域非拘束フィルタ係数修正量を得て、上記周波数領域非拘束フィルタ係数修正量を上記フィルタ係数修正量バッファに保存し、
条件C1に戻り、また、上記条件C1とステップS1、S2、S3、S4、S10、S5、S30、S31、S15と同期、あるいは非同期で、条件C2において更新間隔が経過するのを待ち、更新間隔が経過したなら、
ステップS17において拘束分割番号を決定し、上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域の非拘束フィルタ係数修正量に変換し、上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、
ステップS18において上記フィルタ係数修正量をFGフィルタ係数バッファの内容に加算するか、あるいは、上記フィルタ係数修正量を周波数領域に変換し周波数領域フィルタ係数修正量を得、上記周波数領域フィルタ係数修正量を上記FGフィルタ係数バッファの内容に加算し、ステップS19において更新したフィルタ係数あるいは周波数領域フィルタ係数の拘束分割番号に対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量を再初期値にし、
条件C2に戻る、ことを特徴とする伝達系推定方法。
【請求項１３】
条件C1で参照信号と所望信号がNサンプル得られるまで待ち、参照信号と所望信号がNサンプル得られると、
ステップS1において参照信号を周波数領域に変換することで周波数領域参照信号を求め、
ステップS2において上記周波数領域参照信号を参照信号バッファに蓄積し、
ステップS3において時間領域で参照信号とフィルタ係数をフィルタリングするか、あるいは、それに等価な処理を周波数領域で行うことで推定信号を求め、
ステップS4において上記推定信号と所望信号の差をとることで誤差信号を求め、フィルタ係数更新ステップにおいて上記フィルタ係数が参照信号と所望信号の間の伝達系を模擬し上記誤差信号のパワーが最小となるように上記フィルタ係数を更新し、条件C1に戻る伝達系推定方法において、フィルタ係数更新ステップが、
ステップS10において上記参照信号バッファの内容と周波数領域非拘束フィルタ係数修
正量の内積を周波数毎にとることで周波数領域推定誤差信号を得、
ステップS5において上記誤差信号を周波数領域に変換することで周波数領域誤差信号を求め、
ステップS30において上記周波数領域誤差信号から上記周波数領域推定誤差信号を減じることで周波数領域修正量更新用近似誤差を得、
ステップS31において上記周波数領域修正量更新用近似誤差と上記参照信号バッファの内容から周波数領域修正量更新量を得、
ステップS15において上記周波数領域修正量更新量をフィルタ係数修正量バッファの内容に重み付け加算し上記周波数領域非拘束フィルタ係数修正量を更新し、上記周波数領域非拘束フィルタ係数修正量を上記フィルタ係数修正量バッファに保存し、条件C1に戻り、また、上記条件C1とステップS1、S2、S3、S4、S10、S5、S30、S31、S15と同期、あるいは非同期で、条件C2において更新間隔が経過するのを待ち、更新間隔が経過したなら、
ステップS17において拘束分割番号を決定し、上記拘束分割番号に対応する上記周波数領域非拘束フィルタ係数修正量を時間領域の非拘束フィルタ係数修正量に変換し、上記非拘束フィルタ係数修正量の一部を選択することでフィルタ係数修正量とし、
ステップS18において上記フィルタ係数修正量をFGフィルタ係数バッファの内容に加算するか、あるいは、上記フィルタ係数修正量を周波数領域に変換し周波数領域フィルタ係数修正量を得、上記周波数領域フィルタ係数修正量を上記FGフィルタ係数バッファの内容に加算し、
ステップS19において更新したフィルタ係数あるいは周波数領域フィルタ係数の拘束分割番号に対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量を微少値にし、条件C2に戻る、ことを特徴とする伝達系推定方法。
【請求項１４】
請求項１３記載の伝達系推定方法において、ステップＳ５の直後で、ステップＳ４０において上記周波数領域推定誤差信号を時間領域に変換した際に、時間軸方向の前半部分又は後半部分の上記周波数領域推定誤差信号が０になるようにする拘束処理を近似する行列と、上記周波数領域推定誤差信号との積により、近似拘束周波数領域推定誤差信号を得、ステップＳ３０をステップＳ４１に替えステップＳ４１において上記周波数領域誤差信号から上記近似拘束周波数領域推定誤差信号を減じることで周波数領域修正量更新用近似誤差を得る、ことを特徴とする伝達系推定方法。
【請求項１５】
コンピュータが解読可能なプログラム言語によって記述され、コンピュータに請求項１乃至４又は８乃至１１の何れかに記載の伝達系推定装置として機能させる伝達系推定プログラム。
【請求項１６】
コンピュータが読み取り可能な記録媒体によって構成され、この記録媒体に請求項１５記載の伝達系推定プログラムを記録したことを特徴とする記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、音響エコーキャンセラ、回線エコーキャンセラ、無線および有線の適応等価、劣化画像の元画像の復元、ソナー、レーダー、アンテナアレー受信、マイクロホンアレー収音、その他、参照信号と所望信号の間の伝達系の推定について、適応フィルタの推定誤差のパワーを最小にする伝達系推定装置、方法、プログラム、および記録媒体に関する。
【背景技術】
【０００２】
１．まえがき
推定する伝達系をディジタルフィルタでモデル化し、推定する伝達系への入力信号と伝達系からの出力信号を用いて伝達系を推定する基本的な技術として時間領域で動作する適応フィルタがある。ここでは、まず、時間領域の適応フィルタについて述べ、次に、時間領域の適応フィルタの演算量の低減、収束速度の改善のための技術である周波数領域で動作する適応フィルタについて述べる。背景技術の最後に従来の適応フィルタの演算量と収束特性（収束速度と定常状態での係数誤差）を比較する。
【０００３】
なお、記号の用法は以下のとおりとする。
・時間領域の信号を要素とするベクトルは小文字で表す。また、ベクトルの下付文字で    ベクトルの要素数を示すことがある。要素との区別が必要なときのみ、小文字の後に    ［］を付す。
・時間領域の信号は小文字で表す。ベクトルと区別し時間領域の信号であることを示す    ために小文字の前に＠を付すことがある。
・周波数領域の信号を要素とするベクトルは大文字で表す。また、ベクトルの下付文字    でベクトルの要素数を示すことがある。
・行列はその要素が時間領域の信号でも周波数領域の信号でも大文字で表す。また、行
列の下付文字で行列の次元を示すことがある。例えばXL,Nは行列XがL行N列の行
列であることを示す。
・上付きのT,＊,Hはそれぞれ、転置、複素共役、複素共役転置を表す。
【０００４】
２．時間領域適応フィルタ
まず、時間領域で動作する適応フィルタについて述べる。図１（時間領域適応フィルタ）に適応フィルタを用いた適応制御系のブロック図を示す。図１において、参照信号をx(k)、推定信号をy(k)、所望信号をd(k)、誤差信号をe(k)で表す。この明細書では、時間（時刻）表現は離散時間として整数kで時間を表すものとする。なお、推定する伝達系への入力信号が適応フィルタへの参照信号に、伝達系からの出力信号が適応フィルタの所望信号に相当する。また、適応フィルタはFIR形フィルタであって、時変係数｛h1(k),h2(k),…,hL(k)｝を持つものとする。ただし、Lはフィルタ係数の数である。また、フィルタ係数はベクトルとして、
h［］(k)＝［h1(k),h2(k),…,hL(k)］T
と表す。適応フィルタは、誤差信号e(k)と参照信号x(k)に基づいて、フィルタ係数h(k)を調整して、誤差信号e(k)のパワーを最小にする。この時の係数修正の手順は適応アルゴリズムとして知られている。
【０００５】
次に、適応アルゴリズムについて説明する。適応フィルタは機能的には、例えば図１（時間領域適応フィルタ）に示すように参照信号x(k)がフィルタリング部１１に供給され、フィルタリング部１１において推定信号y(k)が
y(k)＝xLH(k-m)h(k)      (FLT)
と求められる。ただし、
xL［］(k)＝［x(k−L＋1),x(k−L＋2),…,x(k)］T
である。推定信号y(k)と所望信号ｄ（ｋ）との差である誤差信号ｅ（ｋ）が差分器１３で
e(k)＝ｄ(k)-y(k)      (EER)
と求められ、誤差信号ｅ（ｋ）と参照信号ｘ（ｋ）とがフィルタ係数更新部１２に入力され、フィルタリング部１１のフィルタ係数ｈ（ｋ）の更新がなされる。適応アルゴリズムはこのような適応フィルタの内部の演算手順であって、時刻ｋにおいて、フィルタ係数ｈ（ｋ）をフィルタ係数更新部１２で修正して、誤差信号ｅ（ｋ）のパワーが小さくなるような係数ｈ（ｋ）を得る方法である。毎時刻、フィルタ係数ｈ（ｋ）を修正する場合には、その修正は、次式で表される。
h(k)＝h(k-1)+μΔh(k)
ここで、Δｈ（ｋ）は修正ベクトル、μは修正の大きさを制御する量でステップサイズと呼ばれる。一方、主に演算量の低減を目的としてブロック間隔N時刻ごとにフィルタ係数ｈ（ｋ）を修正する場合の修正式は次式で表される。
h(k)＝h(k-N)+μΔh(k)    (ADFUPD)
【０００６】
以降、より一般的なブロック間隔N毎にフィルタ係数ｈ（ｋ）を修正する形式で説明を進める。
修正ベクトルやステップサイズは誤差信号パワーの定義によって異なるが、ここでは、代表的な２つの方法を説明する。
誤差信号パワーＪ（ｋ）を現時刻と過去ｐ-１時刻の現時刻でフィルタ係数を更新した後の誤差の重み付け二乗和として以下のように定義する。
J(k)＝Σ_｛m＝0｝^｛p-1｝λm［d(k-m)-xT(k-m)h(k)］2（Ｊ１）
ここで、Σ_｛m＝0｝^｛p-1｝はｍ＝０からｍ＝ｐ-１に関する総和をとることを意味する。ｐは誤差評価対象時間長であり、この誤差評価対象時間長ｐにより適応アルゴリズムが二つに分類される。λは現時刻に近い誤差により大きい重みをつけるための忘却定数で、０＜λ≦１である。
【０００７】
式（Ｊ１）が最小になるように、修正ベクトルΔｈ（ｋ）を求める。
（ADFUPD）に（Ｊ１）を代入すると、
J(k)＝Σ_｛m＝0｝^｛p-1｝λm［d(k-m)-XT(k-m)｛h(k-N)+Δh(k)｝］2    （Ｊ３）
（Ｊ３）をベクトル行列形式で表すと、
J(k)＝‖Λp,p［dp(k)-XTL,p(k)｛h(k-N)+Δh(k)｝］‖2    （Ｊ４）
となる。ただし、Λp,pは対角行列でその要素をdiag(・)、平方根を√(・)で表すと、
Λp,p＝diag(1，√(λ)，√(λ2)，…，√(λp-1))．
また、
dp［］(k)＝［d(k−p＋1)，d(k−p＋2)，…，d(k)］T
XL,p(k)＝［xL［］(k−p＋1)，xL［］(k−p＋2)，…，xL［］(k)］T
さらに、時刻k-Nのフィルタ係数h(k-N)による誤差を並べた誤差ベクトルep(k|k-N)をep(k|k-N)＝dp(k)-XTL,p(k) h(k-N)
と定義すると、（Ｊ４）は
J(k)＝‖Λp,p［ep(k|k-N)-XTL,p(k) Δh(k)］‖2    （Ｊ５）
となる。
【０００８】
誤差評価対象時間長ｐとフィルタ長Lの大小関係によって修正ベクトルΔh(k)は以下のように二通りに分類される。
（１）  ｐ≦Ｌの場合
未知数の数Lが方程式の数ｐより多い場合には、誤差信号パワーJ(k)を０にする未知数の数が無数に存在するが、通常は、未知数である修正ベクトルΔh(k)を参照信号xL［］(k)，xL［］(k-1)，・・・xL［］(k-p+1)の線形結合として、つまり、修正ベクトルΔh(k)の二乗ノルムが最小となるように決定する。このとき、修正ベクトルΔh(k)は、（Ｊ５）の右辺の二乗ノルムの中を０とおいたΔh(k)に関する方程式の最小二乗ノルム解を求めることで
Δh(k)＝XL,N(k)［XTL,p(k)XL,p(k)］-1ep(k|k-N)
となる。ここではこの修正ベクトルを用いる適応アルゴリズムをブロックアフィン射影法と呼ぶことにする。推定信号y(k)の計算が無駄なく重なりなく行われるようにp＝Nとした、フィルタリングも含めたブロックアフィン射影法は、
yN(k)＝XTL,N(k)h(k-N)        (BAP-FLT)
eN(k|k-N)＝dN(k)-yN(k)      (BAP-ERR)
Δh(k)＝XL,N(k)［XTL,N(k)XL,N(k)］-1eN(k|k-N)        (BAP-GAIN)
h(k)＝h(k-N)+μΔh(k)        (BAP-UPD)
である。
【０００９】
（２）p>Lの場合
未知数の数Lが（独立な）方程式の数ｐより少ない場合には、誤差信号パワーJ(k)を一般には０にすることはできず、誤差信号パワーJ(k)の最小にするΔh(k)は（Ｊ５）の右辺の二乗ノルムの中を０とおいたΔh(k)に関する方程式の最小二乗誤差解として、
Δh(k)＝［XL,p(k)Λ2p,pXTL,p(k)］-1XL,p(k)Λ2p,pep(k|k-N)        （Ｊ６）
と求まる。（Ｊ６）の右辺を時刻kで得られるＮ個の誤差とそれ以外で分けると
Δh(k)＝［XL,p(k)Λ2p,pXTL,p(k)］-1XL,N(k)Λ2N,NeN(k|k-N)+
［XL,p(k)Λ2p,pXTL,p(k)］-1XL,p-N(k)λN-1Λ2p-N,p-Nep-N(k-N|k-N)  （Ｊ７）
となる。ここで、ステップサイズをμ＝１とし、ｐ＝ｋとし、毎N時刻毎にフィルタ係数を更新する場合を考えると、右辺第二項は０になることが以下のようにして分かる。まず、ep-N(k-N|k-N)は、
ep-N(k-N|k-N)＝dk-N(k-N)-XTL,k-N(k-N)｛h(k-2N)+Δh(k-N)｝    （Ｊ８）
となる。
ek-N(k-N|k-2N)＝dk-N(k-N)-XTL,k-N(k-N)h(k-2N),
Δh(k-N)＝［XL,k-N(k-N)Λ2p-N,p-NXTL,k-N(k-N)］-1XL,k-N(k-N)Λ2k-N,k-Nek-N(k-N|k-2N)
であることから、（Ｊ８）は
ep-N(k-N|k-N)＝｛I-XTL,k-N(k-N)［XL,k-N(k-N)Λ2p-N,p-NXTL,k-N(k-N)］-1XL,k-N(k-N)Λ2k-N,k-N｝ek-N(k-N|k-2N)    （Ｊ９）
となる。（Ｊ７）の右辺第二項の
XL,p-N(k)λN-1Λ2p-N,p-Nep-N(k-N|k-N)
に（Ｊ９）を代入すると０になり、右辺第二項が０となることが分かり、修正ベクトルは、
Δh(k)＝［XL,k(k)Λ2k,kXTL,k(k)］-1XL,N(k)Λ2N,NeN(k|k-N)
となる。ここでは、この修正ベクトルを用いる適応アルゴリズムをブロックRecursive Least Squares(RLS)法と呼ぶことにする。フィルタリングも含めたブロックRLS法は、
yN(k)＝XTL,N(k)h(k-N)        (BRLS-FLT)
eN(k|k-N)＝dN(k)-yN(k)      (BRLS-ERR)
Δh(k)＝［XL,k(k)Λ2k,kXTL,k(k)］-1XL,N(k)Λ2N,NeN(k|k-N)    (BRLS-GAIN)
h(k)＝h(k-N)+Δh(k)        (BRLS-UPD)
である。
【００１０】
以上、時間領域の代表的な適応アルゴリズムについて参照信号が単一のチャネルの場合で説明した。参照信号が複数のチャネルである場合には、各チャネルの参照信号ベクトルをチャネル数分だけ連結し一つのベクトルとしたものを新たに参照ベクトルとし、同様に、フィルタ係数もチャネル分連結し一つのベクトルとしたものを新たにフィルタ係数とすることで同様の議論が行えるため省略する。
【００１１】
３．周波数領域適応フィルタ
次に周波数領域の適応フィルタについて説明する。周波数領域で適応フィルタを行なう目的は演算量の低減と、各周波数帯域で適応フィルタを動作させることで参照信号の周波数特性に依存しない速い収束速度を得ることである。
演算量の低減については以下の従来の手法が用いられている。
巡回畳み込みによる直線畳み込み：Discrete Fourier Transform(DFT)の高速演算であるFast Fourier Transform(FFT)等を用いた直線畳み込み。時間領域の適応フィルタで  述べたフィルタリング(BAP-FLT)、(BRLS-FLT)と修正ベクトルの計算(BAP-GAIN)、(BRLS-GAIN)に現れる参照信号からなる行列XL,N(k-N)とベクトル乗算に適用可能である。
巡回畳み込みによる高速なフィルタリングの手法はDiscrete Fourie Transform(DFT)を用いた方法が有名だが、DFTの他にもDiscrete Sine/Cosine Transform(DST,DCT)や、Discrete Hartley Transform(DHT)、Harr Transform、Hadamard Transformなどの直交変換を用いることもできる。ここではDFTについてのみ述べるが上記DCT,DHTを用いる手法については非特許文献１に説明がある。
【００１２】
信号の分割：時刻ｋでの参照信号からなる行列XL,N(k)と１ブロック前の時刻k-NでのXL,N(k-N)は@x(k)、@x(k-1)、…、@x(k-N+1)が含まれるN行を除いて重複しているため、@x(k)、@x(k-1)、…、@x(k-N+1)を含む区間長を一単位として参照信号とフィルタ係数を分割して処理する方が、演算量が少ない、固定小数点演算ではまるめ誤差が少ない、また、誤差信号算出までの遅延時間が短いという点で有利である。これはMulti Delay Filter(MDF)と呼ばれる手法で、詳細は非特許文献２や非特許文献４にある。
また、フィルタリングと修正ベクトルの更新で異なる参照信号の分割をすることも、演算量は増加するが、可能である。本明細書では最も演算量の少ない参照信号をＮ個のデータが重なる２N個データに、フィルタ係数は重なりのないN個ごとに分割する方法で説明を行う。この説明のために、フィルタ係数hをL/N個に分割し、そのm番目の分割をhの第m分割と呼びh(m)と表す。周波数領域のフィルタ係数やフィルタ係数の修正量についても同様である。
【００１３】
共分散行列の周波数領域での近似更新と逆行列の計算：修正ベクトルの計算(BAP-GAIN)、(BRLS-GAIN)に現れる共分散行列XTL,N(k)XL,N(k)、XL,k(k)Λ2k,kXTL,k(k)に適用される。詳細は非特許文献３に譲るが、結果として、参照信号を周波数領域に変換した後、周波数毎に、XTL,N(k)XL,N(k)の場合はパワーの逆数を、XL,k(k)Λ2k,kXTL,k(k)の場合は共分散や相関行列の逆行列を求めることで近似される。
この近似により時間領域で、N行N列、あるいは、L行L列の逆行列を求める演算が、N個のスカラーの逆数、あるいは、L/N行L/N列の行列の逆行列を求める演算になる。N行N列の行列の逆行列を求めるには、N3に比例する演算が必要であるから、時間領域でN3やL3の演算量が、周波数領域ではそれぞれN(N/N)3＝N、N(L/N)3＝L3/N2になるため、大きな演算量の低減が可能となる。
【００１４】
以上の演算量低減手法を時間領域の二つのアルゴリズムに適用し得られる周波数領域のアルゴリズムを以下に示す。記号定義の都合上、まず、修正ベクトルを求める部分を示す。なお、Fは、DFTの処理に相当する行列、F-1は逆DFT処理に相当する行列を表し、ベクトルxのFFT、ベクトルXのIFFTという意味で、それぞれFx、F-1Xと表記する。
（１）周波数領域ブロックアフィン射影法
X2N,2N(k)＝diag(Fx2N(k))        (FDBAP-FX)
U2L,2N(k)＝［X2N,2N(k)，X2N,2N(k-N)，…，X2N,2N(k-L+N)］T  (FDBAP-FXUF)
yN(k)＝WN,2N0 1F-1UT2L,2N(k)H2L(k)                (FD-FLT)
ただし、WN,2N0 1は、周波数領域変換にFFTを用いる場合は、0N,NをN行N列の０行列、INをN行N列の単位行列として、WN,2N0 1＝［0N,N ，IN］であり、その右から乗じられるベクトルの後半分を抜き出す効果がある。
【００１５】
eN(k|k-N)＝dN(k)-yN(k)        (BAP-ERR)
E2N(k)＝F［0TN，eTN(k|k-N)］T        (FDBAP-FE)
R(k)＝UH2L,2N(k)U2L,2N(k)          (FDBAP-R)
あるいは、R(k)＝(1-N/L)R(k-N)+X2N,2N(k)X2N,2N＊(k)      (FDBAP-R2)
ΔH'2L(k)＝U2L,2N(k)｛R-1(k)E2N(k)｝        (FDBAP-GAIN)
Δh'(m)(k)＝F-1ΔH'(m)2N(k)、m＝0,1,…,L/N-1        (FDBAP-DH')
Δh(m)(k)＝WN,2N1 0Δh'(m)2N(k)、m＝0,1,…,L/N-1        (FDBAP-DH)
ただし、WN,2N1 0は、周波数領域変換にFFTを用いる場合は、0N,NをN行N列の０行列、INをN行N列の単位行列として、WN,2N1 0＝［IN ，0N,N］であり、その右から乗じられるベクトルの前半分を抜き出す効果がある。
【００１６】
ΔH(m)2N(k)＝F［ΔhT(m)(k)，0TN］T,m＝0,1,…,L/N-1        (FDBAP-FDH)
H2L(k)＝H2L(k-N)+μΔH2L(k),                      (FDBAP-UPD)
この手法では時間領域のNormalized Least-Mean-Square(NLMS)法を周波数毎に適用したものに等価であり、N＝Lである非特許文献６ではFast LMS、それをMDF構成にした非特許文献２ではMDFと呼ばれている。本明細書では、Frequency-domain Block Affine Projection (FDBAP)と称することにする。
【００１７】
（２）周波数領域ブロックRLS法
X2N,2N(k)＝diag(Fx2N(k))                      (FDBRLS-FX)
U2L,2N(k)＝［X2N,2N(k)，X2N,2N(k-N)，…，X2N,2N(k-L+N)］T 
(FDBRLS-FXBUF1)
U(i)L/N(k)＝［X2N,2N(k)のi行i列,…,X2N,2N(k-L+N)のi行i列］T
i＝1,2,…,L/N    (FDBRLS-FXBUF2)
yN(k)＝WN,2N0 1［F-1UT2L,2N(k)H2L(k)］                (FD-FLT)
eN(k|k-N)＝dN(k)-yN(k)      (BRLS-ERR)
E2N(k)＝F［0TN，eTN(k|k-N)Λ2N,N］T        (FDBRLS-FE)
R(i)(k)＝λNR(i)(k-N)+U(i)L/N(k)U(i)HL/N(k)    (FDBRLS-R)
K'(i)(k)＝R(i)-1(k)U(i)L/N(k),i＝0,1,…,2N-1    (FDBRLS-GAIN1)
K(m)(k)＝diag(K'(0)(k)のm番目の要素,K'(1)(k)のm番目の要素,…,K'(2N-1)(k)のm番目の要素),m＝0,1,…,L/N-1    (FDBRLS-GAIN2)
ΔH'(m)2N(k)＝K(m)(k)E2N(k),m＝0,1,…,L/N-1        (FDBRLS-GAIN3)
Δh'(m)(k)＝F-1ΔH'(m)2N(k)、m＝0,1,…,L/N-1        (FDBRLS-DH')
Δh(m)(k)＝WN,2N1 0Δh'(m)(k)、m＝0,1,…,L/N-1        (FDBRLS-DH)
ΔH(m)2N(k)＝F［ΔhT(m)(k)，0TN］T,m＝0,1,…,L/N-1        (FDBRLS-FDH)
H2L(k)＝H2L(k-N)+ΔH2L(k),                      (FDBRLS-UPD)
この手法は時間領域のRecursive Least Square(RLS)法を周波数毎に適用したものと等価であり、N＝Lでは非特許文献３でEMDFと呼ばれている。本明細書では、Frequency-Domain Block RLS(FDBRLS)と称する。なお、同文献中ではK(m)はFast Transversal Filter (FTF)非特許文献５を用いた高速算法を用いているため表記が上記とは異なる。
【００１８】
なお、上記FDBAPとFDRLSにおいて、参照信号xが実数の場合には周波数領域参照信号X、周波数領域フィルタ係数Hなどで、その要素の約半分が残りの半分の複素共役となる関係がある。この関係を利用すると演算を省略し演算量、使用メモリを減らすことが可能であるが、その実際的手段については、自明なため説明を省略する。
以上代表的な時間領域の二つの手法に適応した周波数領域の手法の修正ベクトルの計算法を示したが、前記演算量低減のための従来手法を用いることで、修正ベクトルの計算が周波数毎に時間領域の従来手法を適用することで行え、同様のことが誤差信号パワーJ(k)の最小化から導出した類似の手法についても言える。従って、以降の説明では主に上記FDBAPを参照しながら行うが、それに特化したものではない。
【００１９】
FDBAPとFDBRLSを含む従来の周波数領域の適応フィルタを用いた伝達系推定装置と方法について説明する。まず、装置から説明する。
図２（周波数領域従来法１）は従来の伝達系推定装置を説明する図であって、第一周波数領域変換部において(FDBAP-FX)等により参照信号を周波数領域に変換することで周波数領域参照信号を求め、上記周波数領域参照信号を参照信号バッファ１２０２に(FDBAP-FXBUF)等により蓄積し、フィルタリング部１１において(FD-FLT)等によりフィルタ係数バッファ１２０９にある周波数領域フィルタ係数と参照信号バッファ１２０２にある周波数領域参照信号から(FD-FLT)等により推定信号を求める。
【００２０】
差分器１３において（BAP-ERR）等により上記推定信号と所望信号の差をとることで誤差信号を求め、第二周波数領域変換部１２０３において（FDBAP-FE）等により上記誤差信号を周波数領域に変換することで周波数領域誤差信号を求め、更新量計算部１２０４において（FDBAP-GAIN）等により上記参照信号バッファ１２０２の内容と上記周波数領域誤差信号から周波数領域非拘束フィルタ係数修正量を求める。
第一時間領域変換部において（FDBAP-DH'）等により上記周波数領域非拘束フィルタ係数修正量を時間領域に変換することで非拘束フィルタ係数修正量を求め、第一要素選択部１２０６において（FDBAP-DH）等により上記非拘束フィルタ係数修正量の要素を取り出すことでフィルタ係数修正量を求め、第三周波数領域変換部１２０７において（FDBAP-FDH）等により上記フィルタ係数修正量を周波数領域に変換することで周波数領域フィルタ係数修正量を求める。
【００２１】
第一平滑部１２０８において（FDBAP-UDP）等により上記周波数領域フィルタ係数修正量と過去の周波数領域フィルタ係数を重み付け加算することで現時刻の周波数領域フィルタ係数を求め、上記周波数領域フィルタ係数でフィルタ係数バッファ１２０９の内容を上書きする。
次に、従来の伝達系推定方法を説明する。図１１（従来法１の方法）は従来の伝達系推定方法を説明するフロー図であって、条件C1で参照信号と所望信号がNサンプル得られるまで待ち、参照信号と所望信号がNサンプル得られると、
ステップS1において（FDBAP-FX）等により参照信号を周波数領域に変換することで周波数領域参照信号を求める。
ステップS2において上記周波数領域参照信号を参照信号バッファに（FDBAP-FXBUF）等により蓄積する。
ステップS3において（BAP-FLT）や（FD-FLT）等によりフィルタ係数バッファにあるフィルタ係数あるは周波数領域フィルタ係数と参照信号あるいは参照信号バッファにある周波数領域参照信号から推定信号を求める。
ステップS4において（BAP-ERR）等により上記推定信号と所望信号の差をとることで誤差信号を求める。
【００２２】
ステップS5において（FDBAP-FE）等により上記誤差信号を周波数領域に変換することで周波数領域誤差信号を求める。
ステップS6において（FDBAP-GAIN）等により上記参照信号バッファの内容と上記周波数領域誤差信号から周波数領域非拘束フィルタ係数修正量を求める。
ステップS7において（FDBAP-DH'）等により上記周波数領域非拘束フィルタ係数修正量を時間領域に変換することで非拘束フィルタ係数修正量を求め（FDBAP-DH）等により上記非拘束フィルタ係数修正量の要素を取り出すことでフィルタ係数修正量を求める。
ステップS8において（BAP-UPD）や（BRLS-UPD）により上記フィルタ係数修正量をフィルタ係数バッファの内容に重み付け加算するか、あるいは、（FDBAP-FDH）等により上記フィルタ係数修正量を周波数領域に変換することで周波数領域フィルタ係数修正量を求め、（FDBAP-UPD）等により上記周波数領域フィルタ係数修正量とフィルタ係数バッファの内容に重み付け加算し、上記重み付け加算結果を上記フィルタ係数バッファの内容に上書きする。
【００２３】
４．低演算量化
主に演算量を低減することを目的としてFDBAP、FDBRLSを修正したものが知られている。
周波数領域非拘束フィルタ係数修正量ΔH'を時間領域に変換し一部を選択し再び周波数領域に変換する処理（FDBAP-DH'）、（FDBAP-DH）を省略して、
ΔH'(m)2N(k)＝H'(m)2N(k-N)+μΔH'(m)2N(k),m＝0,1,…,L/N-1    (UMDF-UPD)
により周波数領域非拘束フィルタ係数H'を更新する手法で、非特許文献２でUnconstrained MDF(UMDF)と呼ばれている。UMDFにおけるフィルタリングはFDBAPの（FD-FLT）と同様に、
yN(k)＝WN,2N0 1F-1UT2L,2N(k)H'2L(k)                (UMDF-FLT)
により行う。周波数領域誤差E'の算出は、FDBAPなら（BAP-ERR）、（FDBAP-FE）、FDBRLSなら（BRLS-ERR）、（FDBRLS-FE）により行う。周波数領域非拘束フィルタ係数修正量ΔH'の算出はFDBAPの場合であれば（FDBAP-R）、（FDBAP-GAIN）により、FDBRLSであれば（FDBRLS-R）、（FDBRLS-GAIN1）、（FDBRLS-GAIN３）により行う。
【００２４】
図３（周波数領域従来法２）はこのUMDFによる装置を説明する図であって、図２（周波数領域従来法１）における第一時間領域変換部、第一要素選択部、第三周波数領域変換部で行う処理を略し、平滑部において（UMDF-UPD）等により上記非拘束周波数領域フィルタ係数を求め、上記周波数領域非拘束フィルタ係数をフィルタ係数バッファに蓄積する。
また、図１２（従来法２の方法）はこの手法による手法を説明する図であって、図１１（従来法１の方法）におけるステップS7とS8を略し、ステップS6の次に（UMDF-UPD）を行うステップS9を設けたものである。
【００２５】
UMDFの（UMDF-UPD）に替えて、周波数領域非拘束フィルタ係数の一部あるいは全部を更新することにし、更新する第ｑ分割については、周波数領域非拘束フィルタ係数修正量を時間領域に変換し一部を選択し再び周波数領域に変換する処理（FDBAP-DH'）、（FDBAP-DH）、（FDBAP-FDH）を
H'(q)2N(k)＝F［{WN,2N1 0F-1［H'(q)2N(k-N)+μΔH'(q)2N(k)］}T，0NT］T,  
(PFDLMS-UPD1)
により行い、その他の分割については、UMDFと同様に、
H'(m)2N(k)＝H'(m)2N(k-N)+μΔH'(m)2N(k),m≠q,        (PFDLMS-UPD2)
により行う手法で、非特許文献７でPartitioned Frequency-Domain LMS(PFDLMS)と呼ばれている。更新する周波数領域フィルタ係数Hの分割の番号は、q＝(k/N)mod(L/N)と均一な機会で更新されるように決める、あるいは、事前にH(m)のどの分割のノルムが大きいか分かれば、ノルムの大きい分割ほどより頻繁にするなど不均一な機会で更新されるように決める。
【００２６】
図４（周波数領域従来法３）はPFDLMSによる装置を説明する図であって、図３（周波数領域従来法２）における第一平滑部１２０８からフィルタ係数バッファ１２０９に向かうデータの流れの一部を第一時間領域変換部１２０５、第一要素選択部１２０６，第三周波数領域変換部１２０７において（PFDLMS-UPD1）右辺のF[{WN,2N10F-1}T，0N,2NT]Tに相当する処理を行い、その他はUMDFと同様に何も処理せずに、フィルタ係数バッファ１２０９の内容を上書きする。拘束部分決定部１２１０は、第１平滑部１２８０から第１時間領域変換部１２０５に向かう周波数領域非拘束フィルタ係数Ｈ’の分割の番号ｑを決定する。分割番号ｑの決定はPFDLMSの説明で述べたようにq＝(k/N)mod(L/N)のように均一な機会で分割が選ばれるようにするか、あるいは、事前にどの分割のノルムが大きいか分かれば、ノルムの大きい分割ほどより頻繁に選択されるようにする。
図１３（従来法３の方法）はPFDLMSによる方法を説明する図であって、図１２（従来法２の方法）に加えて、（PFDLMS-UPD1）右辺に相当する処理を
H'(q)2N(k)＝F［｛WN,2N1 0F-1H'(q)2N(k)｝T，0NT］T        (PFDLMS-UPD3)
により行うステップSaを条件C1とステップS1、S2、S3、S4、S5とS6とは非同期に行うようにしたものである。
【００２７】
５．フィルタ係数二重化との組み合せによる安定化
従来技術を比較する中でも後述するが、UMDFとPFDLMSは、（FDBAP-DH'）と（FDBAP-FDH）の処理を行わないか、行っても更新時にしか行われないため、周波数領域フィルタ係数がそれを時間領域に戻したときに半分が０となるという時間領域の畳み込みを周波数領域で等価に行うための拘束条件が全く満たされないか、更新時でしか満たされない。このためにUMDFとPFDLMSでは推定信号と誤差信号のブロック境界で波形が不連続になりクリック性の雑音が生じるという問題がある。
この問題を解決するために、UMDFやPFDLMS中の周波数領域非拘束フィルタ係数H'の他に新たに、H'に対して（FDBAP-DH'）、（FDBAP-DH）とを行ったフィルタ係数か、あるいは、（FDBAP-DH'）、（FDBAP-DH）に加え（FDBAP-FDH）を行った周波数領域フィルタ係数Hを保存するFGフィルタ係数バッファを設け、推定信号ｙを求めるためのフィルタリングを上記フィルタ係数ｈか上記周波数領域フィルタ係数Hを用いて、（BAP-FLT）か（FD-FLT）により行うというフィルタ係数の二重化をすることが考えられる。本明細書ではフィルタ係数の二重化をしたPFDLMSをdouble buffered PFDLMS（DBPFDLMS）と呼ぶ事にする。以下にDBPFDLMSを詳述する。
【００２８】
DBPFDLMSでは、参照信号を周波数領域に変換し参照信号バッファに保存する処理はFDBAPやFDBRLSと同じく、（FDBAP-FX）と（FDBAP-FXBUF）により行う。また、半固定のフィルタ係数である周波数領域フィルタ係数Hを用いた推定信号yN(k)とeN(k|k-N)の算出は、FDBAPと同様に（FD-FLT）、（BAP-ERR）により行う。フィルタ係数ｈを用いて時間領域でフィルタリングする場合は（BAP-FLT）により行う。毎ブロック、全係数を更新する周波数領域非拘束フィルタ係数H'を用いたフィルタリングはUMDFの（UMDF-FLT）と同様に、
Y'^(k)＝UT2L,2N(k)H'2L(k)                (DB-BGFLT1)
y'^(k)＝F-1Y'^(k)                      (DB-BGFLT2)
y^(k)＝WN,2N0 1y'^(k)                  (DB-BGFLT3)
によりBG側推定信号y^(k)を求める。BG側誤差信号e^(k)は、（BAP-ERR）と同様に、
e^(k)＝d(k)-y^(k)                    (DB-ERR)
により求める。周波数領域の誤差である周波数領域BG側誤差E^(k)は、UMDFと同様に、FDBAPの場合は、
E^(k)＝F［0TN，e^T(k)］T                (DB-FDBAP-FE)
により、FDBRLSの場合は、
E^(k)＝F［0TN，e^T(k)Λ2］T                (DB-FDBRLS-FE)
により変換する。周波数領域非拘束フィルタ係数修正量ΔH'の算出はFDBAPの場合であれば（FDBAP-R）、（FDBAP-GAIN）により、FDBRLSであれば（FDBRLS-R）、（FDBRLS-GAIN1）、（FDBRLS-GAIN3）により行う。周波数領域非拘束フィルタ係数H'および周波数領域フィルタ係数Hの更新は、更新する第ｑ分割については、（PFDLMS-UPD1）により、その他の分割については（PFDLMS-UPD3）により行う。更新した第ｑ分割については、
H(q)(k)＝H'(q)(k)                        (DB-TRANS)
により周波数領域非拘束フィルタ係数からFGフィルタ係数バッファ１２０９Ｆの内容である周波数領域フィルタ係数に複製する。時間領域のフィルタ係数ｈを使用する場合には、（PFDLMS-UPD1）の途中結果として、
h(q)N(k)＝｛WN,2N10F-1［H'(q)2N(k-N)+μΔH'(q)2N(k)］｝  (DB-TD-TRANS)
により得られるフィルタ係数をFGフィルタ係数バッファ１２０９Ｆの内容に上書きする。
【００２９】
図５（周波数領域従来法４）はDBPFDLMSを用いた伝達系推定装置を説明する図である。PFDLMSを用いた装置を説明する図４（周波数領域従来法３）において間引き拘束修正量加算部１２１の出力をフィルタリング部１１に供給することに替え、図５では前記周波数領域非拘束フィルタ係数H'をBGフィルタリング部１１Bに（BGはbackgroundを意味する）供給し、（PFDLMS-UPD1）と（DB-TRANS）により更新された周波数領域フィルタ係数H(q)あるいは（DB-TD-TRANS）により得たフィルタ係数ｈ(q)をFGフィルタ係数バッファ１２０９Ｆ（FGはForegroundの意）の内容に上書きする。
BGフィルタリング部１１Ｂにおいては、乗算器１２１２で上記間引き拘束修正量加算部１２１から供給される周波数領域非拘束フィルタ係数H'と参照信号バッファ１２０２の内容Uを（DB-BGFLT1）により乗じることで、周波数領域BG側推定信号Y'^(k)を生成し、第二時間領域変換部１２１３において上記周波数領域BG推定信号Y'^(k)を時間領域の非拘束BG側推定信号y'^(k)に（DB-BGFLT2）により変換し、第二要素選択部１２１４において時間領域の非拘束BG側推定信号y'^(k)の一部を（DB-BGFLT3）により選択しBG側推定信号y^(k)とし、上記BG側推定信号y^(k)をBGフィルタリング部１１Ｂの出力とする。
【００３０】
また、図４（周波数領域従来法３）において第二周波数領域変換部１２０３への入力を誤差信号eに替えて、図５では上記BG側推定信号y^(k)と所望信号d(k)の差を差分器１２１５で（DB-ERR）により求めたBG側誤差信号e^(k)を求め、第二周波数変換部１２０３において、FDBAPの場合は、（DB-FDBAP-FE）により、FDBRLSの場合は、（DB-FDBRLS-FE）により周波数領域信号の周波数領域BG側誤差E^(k)に変換する。
図１４（従来法４の方法）は、DBPFDLMSによる方法を説明する図であって、
PFDLMSによる従来法を説明する図１３（従来法３の方法）のステップS4の次に（DB-BGFLT1）を行うS3B1と（DB-BGFLT2）と（DB-BGFLT3）を行うS3B2と（DB-ERR）を行うステップS4Bを追加し、ステップS5を、（DB-FDBAP-FE）あるいは（DB-FDBRLS-FE）を行うステップS5Bで置き換え、さらに、ステップSaの後に（PFDLMS-UPD3）の結果の周波数領域フィルタ係数H'か、（PFDLMS-UPD3）の途中結果をFGフィルタ係数バッファに上書きするステップSaFを追加する。
【００３１】
６．従来法の比較
以上示した従来法の中で、時間領域の手法として最も一般的なブロックアフィン射影法でN＝p＝1とした通常Normalized Least-Mean-Square(NLMS)と呼ばれるもの、周波数領域の手法として周波数毎の演算が簡単なFDBAP、UMDF、PFDLMSとDBPFDLMSについて、演算量と収束特性（フィルタ係数誤差の収束速度、定常状態誤差レベル）を比較する。
○演算量
周波数領域の手法の演算の多寡は主に周波数領域変換、時間領域変換で用いられるFFT、Inverse FFT（IFFTと略）と周波数領域での複素積の回数で決まる。ブロック間隔Nあたりのこれらの回数はFDBAPの説明と各手法の説明より表１のとおりになる。なお、PFDLMSでの（PFDLMS-UPD1）を行う間隔はLとしている。
【００３２】
表１
FFT,IFFT            複素積          実数積
NLMS          ０                    ０              ２NL
FDBAP        ３＋２（L/N）        ４L+４N          ０
UMDF          ３                  ４L+４N          ０
PFDLMS        ３＋２              ４L+４N          ０
DBPFDLMS    ３＋２＋１          ４L+２L+４N        ０
表１でFFT、IFFTの回数のうち、NLMS以外の全てに共通な３回は、第一周波数領域変換部１２０１、第二周波数領域変換部１２０３とフィルタリング部１１での演算によるものであり、FDBAPの＋２（L/N）回とPFDLMS、DBPFDLMSの＋２回は、第一時間領域変換部１２０５と第三周波数領域変換部１２０７での演算によるものであり、DBPFDLMSの＋１回は第二時間領域変換部によるものである。
【００３３】
また、４Ｌ回の複素積は、フィルタリング部１１での２Ｌ回と更新量計算部１２０４での２Ｌ回の複素積によるものであり、DBPFDLMSの＋２Ｌ回は乗算器１２１２での複素積によるものである。４Ｎ回の複素積は（FDBAP-R2）の右辺第二項と（FDBAP-GAIN）の２N回とR-1(k)E2N(k)の２N回の複素積による。FFT、IFFTに要する複素積の数を、2(2N)log2(2N)として従来の各種法を比較する。
まず、時間領域と周波数領域の手法を比較すると、時間領域の手法NLMSは演算量がNとLの積に比例するのに対し、周波数領域の手法は、最大であるFDBAPでもLとlog2Nの積に比例する程度であり、周波数領域の手法の方が時間領域の手法に比べて(log2N)/N倍と少ない。
【００３４】
次に周波数領域の手法間を比較すると、UMDF、PFDLMSとDBPFDLMSは、適応フィルタ長Lに依存せず、FFT、IFFT回数が３回から６回であるのに対し、FDBAPでは、適応フィルタ長が長くなるに従いFFT、IFFT回数が増加するので、FDBAPよりUMDF、PFDLMS、DBPFDLMSの方が演算量が少ない、また、複素積の数はDBPFDLSMが他の周波数領域の手法に比べて２L回多いので、適応フィルタ長Lが長くなるにつれてUMDFが最も演算量が少なく、次いでPFDLMS、DBPFDLMS、FDBAPの順で演算量が少ない。
【００３５】
○収束特性
比較のための試験条件は以下のとおりである。
伝達系のインパルス応答長：  ２０００
（時間領域）適応フィルタ長L：    ２５６０
ブロック間隔  N  ：１２８
FFT点数：  ２５６
ステップサイズμ：  NLMSは０．２５、それ以外は０．５、
参照信号：  （１）白色ガウス雑音、  （２）男性音声
その他：伝達系の変化を模擬するため、２００００点で伝達系のインパルス応答を、符号を反転することで変化させた。
【００３６】
試験の結果を図２１、２２、２３（従来手法の係数誤差収束特性の比較１、２、３）に示す。図２１と図２２の横軸はサンプル時間を、縦軸は、時間領域の適応フィルタ係数と伝達系のインパルス応答の差の二乗和を伝達系のインパルス応答の二乗和で除した値であり、小さければフィルタ係数の誤差が少なくよい推定ができているといえる。
図２１に示す参照信号が白色ガウス雑音の場合、破線で示したUMDF以外はほぼ同じ収束速度、定常状態誤差レベルであり、UMDFだけが収束速度が遅く、また、定常状態誤差レベルも大きい。なお、UMDFの定常状態誤差レベルが大きい理由は、他の周波数領域の手法にはある（FDBAP-DH'）、（FDBAP-DH）と（FDBAP-FDH）が無いため、（FD-FLT）による周波数領域で行うフィルタリングと時間領域でのフィルタリングの結果が一致しないためである。
【００３７】
図２２に示す参照信号が男性音声の場合、NLMSの収束が他より遅く、NLMS以外では、PFDLMS、FDBAP、UMDFの順で２，３dB程度ずつ定常誤差が小さくなっている。なお、DBPFDLMSとPFDLMSの収束曲線はほぼ重なっている。また、UMDFとPFDLMSではブロックの境界で波形が不連続になることに起因するクリック雑音が発生する。図２３（従来手法の係数誤差収束特性の比較３）に各手法の第２１７ブロックから第２１８ブロックの境界付近の誤差信号を示す。DBPFDLMSでは不連続は見られないが、LMDFとPFDLMSではブロック境界で波形が不連続になっていることが分かる。
【非特許文献１】"An embedding approach to frequency-domain and subband adaptive filtering,"Ricardo Merched, Ali H. Sayed, IEEE Transaction on Signal Processing, vol. 48, No.9, pp.2607-2619 (2000.9)
【非特許文献２】"Multidelay block frequency domain adaptive filter,"Jia Sien Soo and Khee K. Pang, IEEE Transaction on Acoustics, Speech, and Signal Processing, vol. 38, No. 2, pp. 373-376(1990.2)
【非特許文献３】"Frequency-domain adaptive filtering revisited, generalization to the multi-channel case, and application to acoustic echo cancellation," Jacob Benesty, and Dennis R. Morgan, IEEE ICASSP-2000, pp.789-792
【非特許文献４】"The generalized multidelay adaptive filter: structure and convergence analysis," E.Moulines, O. Ait Amrane, and Y. Grenier, IEEE Transaction of Signal Processing, vol. 43, pp. 14-28, (1995.1)
【非特許文献５】"Fast recursive-least-squares transversal filters for adaptive filtering,"J. M. Cioffi and T. Kailath, IEEE Transaction on Acoustics, Speech, and Signal Processing, vol. 32, pp. 304-337 (1984.4)
【非特許文献６】"Fast implementation of LMS adaprive filters," E. R. Ferrara, IEEE Transaction on Acoustics, Speech, and Signal Processing, vol. 28, pp. 474-475 (1980.8)
【非特許文献７】"Connecting Partitioned Frequency-Domain Filters in Parallel or in Cascade, section I, II, III," M. Joho and G. S. Moschytz, vol 47. pp. 685-698 (2000.8)
【発明の開示】
【発明が解決しようとする課題】
【００３８】
DBPFDLMSでは、所望信号ｄ中の参照信号ｘに起因しない外乱成分が増加する、エコーキャンセラでは同時通話と呼ばれるような状況において、周波数領域非拘束適応フィルタ係数Ｈ’が乱れ、乱れた非拘束周波数領域適応フィルタ係数Ｈ’がＦＧフィルタ係数バッファ１２０９Ｆの内容である周波数領域フィルタ係数Ｈに複製されると、乱れた非拘束周波数領域フィルタ係数Ｈ’が乱れた影響を大きく受けて誤差信号ｅが大きくなる問題がある。また、DBPFDLMSは他の従来法に比べ、収束が速く定常消去量が小さく、かつ、雑音が誤差信号に生じないが、演算量は最小のものより、ブロック長Nの２倍である２N点のIFFTが１回と複素積がフィルタ長Lの２倍である２L回多い。この演算量の増加分だけ、伝達係推定装置または方法の規模が大きくなり高価になるか、同一の規模ならば性能が適応フィルタ長を削減したための性能劣化が生じるという問題がある。以下に示すこの発明による解決手段は、外乱による収束特性劣化の緩和と増加した１回のIFFTを削減することを目的としている。
【課題を解決するための手段】
【００３９】
本発明の理解を容易にするために、ここではDBPFDLMSと数学的には等価で同じ収束性を持つが、周波数領域非拘束フィルタ係数H'を毎ブロック更新することに替え、周波数領域非拘束フィルタ係数修正量ΔH'を更新するようにした修正量蓄積型DBPFDLMSについて説明する。修正量蓄積型DBPFDLMSに変形する理由は、通常フィルタ係数よりフィルタ係数修正量のノルムが小さいため、発明が解決しようとする課題の項で述べた外乱や、IFFTを除去する直接の影響がよりノルムの小さいフィルタ係数修正量にのみ及ぶようにすることで収束特性の劣化を抑えられることと、後述するように周波数領域非拘束フィルタ係数修正量ΔH'は繰り返しリセットされるので発明が解決しようとする課題の項で述べた外乱や、IFFTを除去した影響の蓄積を避けられることである。
【００４０】
修正量蓄積型DBPFDLMSでは、参照信号を周波数領域に変換し参照信号バッファに保存する処理はFDBAPやFDBRLSと同じく、（FDBAP-FX）と（FDBAP-FXBUF）により行う。また、半固定のフィルタ係数である周波数領域フィルタ係数を用いた推定信号とyN(k)と誤差eN(k|k-N)の算出は、FDBAPと同様に（FD-FLT）、（BAP-ERR）により行う。毎ブロック、全係数を更新する周波数領域非拘束フィルタ係数修正量ΔH'を用いたフィルタリングは
ΔY'^(k)＝UT2L,2N(k)ΔH'(k)                (DB2-BGFLT1)
Δy'^(k)＝F-1ΔY'^(k)                      (DB2-BGFLT2)
Δy^(k)＝WN,2N0 1Δy'^(k)                  (DB2-BGFLT3)
により行い、推定誤差信号Δy^(k)を求める。次に、誤差eN(k|k-N)と推定誤差信号Δy^(k)の差分である修正量更新用誤差Δe^(k)を
Δe^(k)＝eN(k|k-N)-Δy^(k)                (DB2-ERR)
により求める。周波数領域の修正量更新誤差ΔE^(k)は、UMDFと同様に、FDBAPの場合は、
ΔE^2N(k)＝F［0TN，Δe^T(k)］T                (DB2-FDBAP-FE)
により、FDBRLSの場合は、
ΔE^2N(k)＝F［0TN，Δe^T(k)Λ2］T                (DB2-FDBRLS-FE)
により求める。
【００４１】
周波数領域非拘束フィルタ係数修正量の修正量である周波数領域修正量更新量ΔΔH'はFDBAPの場合であれば（FDBAP-R）と（FDBAP-GAIN）で、ΔH'とΔΔH'、EとΔE^を入れ替えた
ΔΔH'2L(k)＝U2L,2N(k)｛R-1(k)ΔE^2N(k)｝          (DB2-FDBAP-GAIN)
により行う。FDBRSであれば（FDBRLS-R）、（FDBRLS-GAIN1）と（FDBRLS-GAIN3）においてE^とΔE^、ΔH'とΔΔH'を置き換えた
ΔΔH'(m)2N(k)＝K(m)(k)ΔE^2N(k)、m＝0,1,…,L/N-1    (DB2-FDBRLS-GAIN3)
により行う。周波数領域非拘束フィルタ係数修正量ΔH'の更新は、FDBAPを用いている場合は、（FDBAP-UPD）においてH2LをΔH'2Lと、ΔH2LをμΔΔH'2Lと入れ替えた
ΔH'2L(k)＝ΔH'2L(k-N)+μΔΔH'2L(k),                  (DB2-FDBAP-UPD)
により行う。FDBRLSを用いている場合には、（FDBRLS-UPD）においてH2LとΔH'2Lと、ΔH2LをΔΔH'2Lと入れ替えた
ΔH'2L(k)＝ΔH'2L(k-N)+ΔΔH'2L(k),                  (DB2-FDBRLS-UPD)
により行う。上記フィルタリング部１１でフィルタリングを時間領域で行う場合は、フィルタ係数ｈの拘束分割番号ｑに対応する分割h(q)の更新は、(FDBAP-DH')、(FDBAP-DH)でm＝qとした
Δh'(q)(k)＝F-1ΔH'(q)(k)                              (DB2-DH')
Δh(q)(k)＝WN,2N1 0Δh'(q)(k)                            (DB2-DH)
により、フィルタ係数修正量Δhに変換し
h(q)(k)＝h(q)(k(q)')+ηΔh(q)(k)                    (DB2-TD-TRANS)
により行う。ここでk(q)'は前回フィルタ係数h(q)(k)を更新した時刻であり、ηは更新調整量であり、定数としても時刻kと拘束分割番号に依存してもよく、通常は０から２の範囲の値をとる。上記フィルタリング部１１ではフィルタリングを周波数領域で行う場合の周波数領域フィルタ係数Hの上記拘束分割番号qに対応する分割H(q)の更新は、（DB2-DH）の後、（FDBAP-FDH）でm＝qとした
ΔH(q)(k)＝F［ΔhT(q)(k)， 0TN］T                        (DB2-FDH)
によりフィルタ係数修正量Δhから周波数領域フィルタ係数修正量ΔH(q)に変換し、
H(q)2N(k)＝H(q)2N(k)+ηΔH(q)2N(k(q)')                    (DB2-TRANS)
により更新する。さらに、
ΔH'0(q)＝α1H(q)2N(k)+α2ΔH(q)2N(k)+α3ΔH'(q)2N(k(q)')+β［］  
(DB2-HCLR)
なる再初期値ΔH'0(q)を周波数領域フィルタ係数修正量ΔH'0(q)に設定する。ここでα1、α2及びα3は０に近い微少値、β［］はその要素が０に近いベクトルである。
【発明の効果】
【００４２】
本発明法は従来法に比べて大きな収束特性の劣化がなく、演算量をIFFT一回分削減することが可能であり、その結果、従来法を用いた伝達系推定装置あるいは方法に比べて規模が小さくなり廉価になるか、適応フィルタ長が長くして規模を同じにするなら性能の向上が得られるという効果がある。
【発明を実施するための最良の形態】
【００４３】
本発明による伝達系推定装置は全てをハードウェアによって構成することもできるが、それよりもむしろ、コンピュータに本発明による伝達系推定プログラムをインストールし、コンピュータに、伝達系推定装置として機能させる実施形態が最良である。
図６に示す第一発明、図８に示す第二発明、図９に示す第三発明はそれぞれ、図７に示す拘束部１２５０を備えた間引き拘束修正量更新量加算部を用いた修正量蓄積型DBPFDLMSで構成される伝達系推定装置の実施例を示す。
同じく図６に示す第四発明、図８に示す第五発明、図９に示す第六発明はそれぞれ、図１０に示す拘束部を除去した間引き拘束修正量更新量が加算部を用いた修正量蓄積型DBPFDLMSで構成される伝達系推定装置の実施例を示す。
【実施例１】
【００４４】
図６及び図７を用いてこの発明の第一発明の実施例を説明する。図６は上述した修正量蓄積型DBPFDLMSによる伝達系推定装置の実施例、図７はこの伝達系推定装置を構成する間引き拘束修正量更新量加算部の実施例を示す。伝達系推定装置１はフィルタリング部において時間領域で参照信号xとフィルタ係数hを（BAP-FLT）によりフィルタリングするか、あるいは、それに等価な処理を周波数領域で（FD-FLT）により行うことで推定信号yを求め、差分器１３では上記推定信号yと所望信号dの差をとることで誤差信号eを（BAP-ERR）等により求め、フィルタ係数更新部１２においては上記フィルタ係数Hが参照信号xと所望信号dの間の伝達系を模擬し上記誤差信号eのパワーが最小になるように上記フィルタ係数Hを更新する。
【００４５】
上記フィルタ係数更新部１２において、第一周波数領域変換部１２０１において参照信号x(k)を（FDBAP-FX）により、周波数領域の周波数領域参照信号X2N(k)に変換し、参照信号バッファ１２０２において現時刻の上記周波数領域参照信号Xを蓄積することで参照信号バッファ１２０２の内容Uを（FDBAP-FXBUF）により更新し、乗算器１２１２で上記参照信号バッファの内容Uとフィルタ係数修正量バッファ１２２０（図７）の内容である周波数領域非拘束フィルタ係数修正量ΔH'(k)の内積を（DB2-BGFLT1）により求め、その結果として周波数領域推定誤差信号ΔY'^(k)を生成する。
【００４６】
第二時間領域変換部１２１３において上記周波数領域推定誤差信号ΔY'^(k)を時間領域の選択前推定誤差信号Δy'^(k)に（DB2-BGFLT2）により変換し、第二要素選択部１２１４において上記選択前推定信号Δy'^(k)の一部を（DB-BGFLT3）により選択し推定誤差信号Δy^(k)とし、差分器１２１５で上記推定誤差信号Δy^(k)と誤差信号e(k)の差である修正量更新用誤差Δe^(k)を（DB2-ERR）により求め、第二周波数領域変換部１２０３において上記修正量更新用誤差Δe^(k)をFDBAPの場合は、（DB2-FDBAP-FE）により、FDBRLSの場合は、（DB2-FDBRLS-FE）により周波数領域信号の周波数領域修正量更新量用誤差ΔE^(k)に変換する。
【００４７】
更新量計算部１２０４において上記周波数領域修正用更新用誤差ΔE^(k)と上記参照信号バッファ１２０２の内容U2L,2N(k)から、利用する適応フィルタ手法に従って周波数領域修正量更新量ΔΔH'(k)を求め、間引き拘束修正量更新量加算部１２２において、第一平滑部１２０８（図７）において上記周波数領域非拘束フィルタ係数修正量ΔH'(k)に上記周波数領域修正量更新量ΔΔH'(k)を利用する適応フィルタ手法に従って（DB2-FDBAP-UPD）あるいは（DB2-FDBRLS-UPD）により加算し、第一拘束部分決定部１２１０（図７）で第一拘束分割番号qを決定する。
【００４８】
上記第一平滑部１２０８（図７）の出力のうち上記第一拘束分割番号q以外のうち、第二拘束部分決定部１２１０' で任意の番号の中から決定された第二拘束分割番号ｑ’については、拘束部１２５０において（FDBAP-DH'）、（FDBAP-DH）と（FDBAP-FDH）により時間領域に変換した際に時間軸方向の前半部分又は後半部分の係数が０になるようにした後、第二拘束分割番号ｑ’以外の部分についてはそのまま、上記フィルタ係数修正量バッファ１２２０内の上記周波数領域非拘束フィルタ係数修正量ΔH'(k)に上書きする。第一時間領域変換部１２０５において上記第一平滑部１２０８の出力のうち上記拘束部分qに対応するΔH'(q)(k)を（DB2-DH'）により非拘束フィルタ係数修正量Δh'(q)に変換し、第一要素選択部１２０６において上記非拘束フィルタ係数修正量Δh'(q)を（DB2-DH）によりフィルタ係数修正量Δh(q)に変換し第二平滑部１２１１において上記FGフィルタ係数バッファ１２０９Ｆの内容に（DB2-TD-TRANS）により加算しFGフィルタ係数バッファ１２０９Ｆに保存するか、あるいは、第三周波数変換領域変換部１２０７において上記フィルタ係数修正量Δh(q)を（DB2-FDH）により周波数領域フィルタ修正量ΔH(q)に変換し、上記周波数領域非拘束フィルタ係数修正量ΔH'(k)と周波数領域フィルタ修正量ΔH(q)を間引き拘束修正量更新量加算部１２２の出力とし、第二平滑部１２１１において上記周波数領域フィルタ修正量ΔH(q)を周波数領域フィルタ係数H(q)に（DB2-TRANS）により加算しFGフィルタ係数バッファ１２０９Ｆに保存する。さらに再初期値計算部１２３０において（DB2-HCLR）により微少値の再初期値を求め、その値で上記フィルタ係数修正量バッファ１２２０内の上記周波数領域非拘束フィルタ修正量ΔH'(q)を上書きする。
【００４９】
図１５は修正量蓄積型DBPFDLMS法による第一発明の方法を説明するフローチャートである。条件C1で参照信号xと所望信号dがNサンプル得られるまで待ち、参照信号xと所望信号dがNサンプル得られると、
ステップS1において（FDBAP-FX）等により参照信号xを周波数領域に変換することで周波数領域参照信号Xを求める。
ステップS2において上記周波数領域参照信号Xを参照信号バッファの内容Uに（FDBAP-FXBUF）等により蓄積する。
【００５０】
ステップS3において（BAP-FLT）や（FD-FLT）等によりフィルタ係数バッファにあるフィルタ係数hあるいは周波数領域フィルタ係数Hと参照信号xあるいは参照信号バッファにある周波数領域参照信号Uから推定信号yを求める。
ステップS4において（BAP-ERR）等により上記推定信号yと所望信号dの差をとることで誤差信号eを求める。
ステップS10においてフィルタ係数修正量バッファにある周波数領域非拘束フィルタ係数修正量ΔH'(k)と参照信号バッファの内容Uから（DB2-BGFLT1）により周波数領域推定誤差信号ΔY'^を求める。
【００５１】
ステップS11において上記周波数領域推定誤差信号ΔY'^を（DB2-BGFLT2）と（DB2-BGFLT3）により推定誤差信号Δy^を求める。
ステップS12において上記推定誤差信号Δy^と誤差信号eの差として（DB2-ERR）により修正用誤差信号Δe^を求める。
ステップS13において上記修正用誤差信号Δe^を（DB2-FE）により周波数領域に変換し周波数領域修正用誤差ΔE^を求める。
ステップS14において上記周波数領域修正用誤差ΔE^と参照信号バッファの内容Uから（DB2-FDBAP-GAIN）等により周波数領域修正量更新量ΔΔH'を求める。
【００５２】
ステップS15において上記周波数領域修正量更新量ΔΔH'をフィルタ係数修正量バッファの内容ΔH'に（DB2-FDBAP-UPD）等により重み付け加算し周波数領域非拘束フィルタ係数修正量ΔH'を得て上記周波数領域非拘束フィルタ係数修正量ΔH'を上記フィルタ係数修正量バッファに保存し、ステップS16において第二拘束分割番号q'を決定し、上記周波数領域非拘束フィルタ係数修正量ΔH'のうち第二拘束分割番号q'に対応する部分ΔH'(q')を(FDBAP-DH')、(FDBAP-DH)と(FDBAP-FDH)により時間領域に変換した際に時間軸方向の前半部分又は後半部分の係数が０に成るようにした後、上記フィルタ係数修正量バッファに保存し、条件C1に戻るループと、上記条件C1とステップS1、S2、S3、S4、S10、S11、S12、S13、S14、S15、S16と同期あるいは非同期に条件C2で更新間隔が経過するのを待ち、更新間隔が経過したなら、
ステップS17において拘束分割番号qを決定し上記周波数領域非拘束フィルタ係数修正量ΔH'のうち上記拘束分割番号qに対応する部分ΔH'(q)を（BD2-DH'）により時間領域の非拘束フィルタ係数修正量Δh'(q)に変換し、上記非拘束フィルタ係数修正量Δh'(q)の一部を（DB2-DH）により選択することでフィルタ係数修正量Δh(q)とする。
【００５３】
ステップS18において上記フィルタ係数修正量Δh(q)を上記FGフィルタ係数バッファの内容に（DB2-TD-TRANS）により加算するか、あるいは、上記フィルタ係数修正量を（DB2-FDH）により周波数領域に変換し周波数領域フィルタ係数修正量ΔH(q)を得、上記周波数領域フィルタ係数修正量ΔH(q)を上記FGフィルタ係数バッファの内容に（DB2-TRANS）により加算する。
ステップS19において更新したフィルタ係数hあるいは周波数領域フィルタ係数Hの拘束分割番号qに対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量ΔH'(q)(k)を（DB2-HCLR）により得られる再初期値にし、条件C2に戻るループからなる。
ここで第二発明及び第三発明の概要を予め説明する。
【００５４】
第二の発明
修正量蓄積形DBPFDLMS法において周波数領域推定誤差信号ΔY'^から推定誤差信号Δy^に変換し、時間領域で（DB2-ERR）により修正量更新誤差Δe^を求めることに替えて、この第二発明では周波数領域で
ΔE'^(k)＝E(k)-ΔY'^(k)                (CLAIM1-FE)
により、（DB2-FE）左辺の周波数領域修正量更新誤差ΔE^に相当する周波数領域修正量更新用近似誤差ΔE'^を求めることで、IFFTを１回削減するようにしたのが第二の発明である。
【００５５】
なお、周波数領域修正量更新用近似誤差を用いた周波数領域修正量更新量
ΔΔH'は、FDBAPの場合は、（FDBAR-R）を行った後に、
ΔΔH'(k)＝U2L,2N(k)｛R-1(k)E^2N(k)｝        (CLAIM-FDBAP-GAIN)
によりFDBRLSの場合は、Λ＝Iの場合のみ計算可能で、（FDBRLS-R）、（FDBRLS-GAIN1）と（FDBRLS-GAIN2）を行った後に、
ΔΔH'(m)2N(k)＝K(m)(k)ΔE^2N(k),m＝0,1,…,L/N-1        (FDBRLS-dGAIN3)
により求められる。
本第二の発明は、修正量蓄積型DBPFDLMSに比べIFFT一回分の演算が少ないが、その一方で、後述するように、（DB2-BGFLT2）、（DB2-BGFLT3）、（DB2-ERR）、（DB2-FE）の一連の演算を大幅に省略したために収束速度が遅くなるという欠点がある。
【００５６】
第三の発明
第三の発明は、第二の発明における収束速度の劣化の欠点を若干の演算量の増加と引き換えに緩和するものである。第二の発明で省略した修正量蓄積型DBPFDLMSの周波数領域修正量更新用近似誤差ΔE'を求めるための処理（DB2-BGFLT2）、（DB2-BGFLT3）、（DB2-ERR）、（DB2-FDBAP-FE）、（DB2-FDBRLS-FE）は周波数領域修正量更新用近似誤差ΔE^が時間領域で要素の半分が０になるという拘束を課すものであることから拘束処理と呼ぶことにする。拘束処理はFDBAPの場合、（DB2-BGFLT2）、（DB2-BGFLT3）と（DB2-FDBAP-FE）より、
G＝F［0NT，(WN,2N0 1F-1)T］T
であり、FDBRLSの場合、（DB2-BGFLT2）、（DB2-BGFLT3）と（DB2-FDBRLS-FE）より、
G＝F［0NT，(Λ2WN,2N0 1F-1)T］T
である。この行列Gは、非特許文献３にあるように対角要素が優勢な巡回行列である。巡回行列の１列目の先頭と底に近い要素だけを窓掛けして抜き出し、他の要素を０としたベクトルから構成した巡回行列によってGを近似すれば精度をあまり落とさずに演算量を低減できる。さらに、Gの１列目の中央に向かってなめらかに０になるようにする必要がある。これは、近似したGを＃Gと表すと、１列目の中央に向かってなめらかに０になるように窓掛けする方が、F-1＃GFとWN,2N0 1との誤差の最大値が小さくなるためである。例えば、Hanning窓により、2N要素がある１列目の先頭D+1（D＜N/2）要素と底のD要素を取り出すには、m番目（0≦m≦2N-1）の要素に
S｛1+cos(πm/D)｝，0≦m≦D
S｛1+cos(π(2N-m)/D)｝，2N-D≦m≦2N-1
0,                    D＜m＜2N-D
とする。ここではSは窓振幅で、窓掛けによる信号のレベルの低下を補償するためのものである。上記Hanning窓の場合、窓振幅Sは通常２であるが、２以下の正数でもかまわない。また値が複素数の要素は実数部が虚数部より数倍以上大きいので、虚数部を無視し０としてもよい。Dの値はNに依存するが、典型的な例としては、N＝128のときにD＝5である。
【００５７】
行列＃Gを用いる拘束処理を近似拘束処理と呼ぶこととする。第三の発明は、第二の発明において、周波数領域推定誤差ΔY'^を近似拘束処理し、近似拘束周波数領域推定誤差信号＃ΔY'^を求め、次に、周波数領域修正量更新用近似誤差ΔE^を求めるものである。すなわち、（CLAIM1-FE）に替えて、
＃ΔY'^(k)＝＃GΔY'^(k)                (CLAIM2-CNSTRAN)
ΔE'^(k)＝ΔE(k)-＃ΔY'^(k)              (CLAIM2-FE)
を用いる。
【実施例２】
【００５８】
第二の発明のよる伝達系推定装置の実施例を、図７と図８（発明２による装置）を用いて説明する。一般に、伝達系推定装置は、フィルタリング部１１と差分器１３とフィルタ係数更新部１２を有する。フィルタリング部１１において時間領域で参照信号xとフィルタ係数hを（BAP-FLT）によりフィルタリングするか、あるいは、それに等価な処理を周波数領域で（FD-FLT）により行うことで推定信号yを求め、差分器１３では上記推定信号yと所望信号dの差をとることで誤差信号eをe＝d-yにより求め、フィルタ係数更新部１２においては上記フィルタ係数Hが参照信号xと所望信号dの間の伝達系を模擬し上記誤差信号eのパワーが最小となるように上記フィルタ係数hあるいは周波数領域フィルタ係数Hを更新する。
【００５９】
第二の発明の伝達系推定装置１の実施例では、上記フィルタ係数更新部１２が、以下の処理を有することを特徴とする。つまり、
（１）第一周波数領域変換部１２０１において上記参照信号xを（FDBAP-FX）により周波数領域の周波数領域参照信号Xに変換する。ただし、周波数領域への変換は、従来の技術で述べたように、DFT、DCT、DSTやDHTなどの直交変換を用いることもできる。
（２）参照信号バッファ１２０２において現時刻の上記周波数領域参照信号Xを（FDBAP-FXBUF）により蓄積することで参照信号バッファ１２０２の内容Uを更新する。
（３）乗算器１２１２で上記参照信号バッファ１２０２の内容Uとフィルタ係数修正量バッファ１２２０（図７）の内容である周波数領域非拘束フィルタ係数修正量ΔH'の内積を（DB2-BGFLT1）により周波数毎に取りその結果として周波数領域推定誤差信号ΔY'^を生成する。
（４）第二周波数領域変換部１２０３において（FDBAP-FE）により上記誤差信号eを周波数領域信号の周波数領域誤差Eに変換する。
（５）周波数領域差分器１２３０において上記周波数領域誤差信号Eと上記周波数領域推定誤差信号ΔY'^との差を（CLAIM1-FE）によりとることで周波数領域修正量更新用近似誤差ΔE'^を求める。
【００６０】
（６）更新量計算部１２０４において上記周波数領域修正量更新用近似誤差ΔE'^と上記参照信号バッファの内容Uから周波数領域修正量更新量ΔΔH'を求める。周波数領域修正量更新量ΔΔH'は、適応アルゴリズムによって異なり、FDBAPの場合は、（FDBAR-R）を行った後に（CLAIM-FDBAP-GAIN）により、FDBRLSの場合は、Λ＝Iの場合のみ計算可能で、（FDBRLS-R）、（FDBRLS-GAIN1）と（FDBRLS-GAIN2）を行った後に（FDBRLS-dGAIN3）により求める。
（７）間引き拘束修正量更新量加算部１２２において上記非拘束フィルタ係数修正量ΔH'に上記周波数領域修正量更新量ΔΔH'を以下のように加えることで更新するとともに周波数領域フィルタ係数修正量ΔHを得る。第一平滑部１２０８（図７）において上記周波数領域非拘束フィルタ係数修正量ΔH'に上記周波数領域修正量更新量ΔΔH'を(DB2-FDBAP-UPD)あるいは(DB2-FDBRLS-UPD)により重み付け加算する。第一拘束部分決定部１２１０において第一拘束分割番号qを決定する。qの決定方法は、q＝(k/N)mod(L/N)と均等な機会で更新されるように決める、あるいは、事前にH(m)のどの分割ノルムが大きいかが分かれば、ノルムの大きい分割ほどより頻繁にするなど不均一な機会で更新されるように決める。上記第一平滑部１２０８における重み付け加算結果である上記周波数領域非拘束フィルタ係数修正量ΔH'の一部あるいは全部分のうち、第二拘束部分決定部１２１０'で任意の番号から選ばれて決定された第二拘束分割番号q'については、拘束部１２５０において(FDBAP-DH')、(FDBAP-DH)と(FDBAP-FDH)により時間領域に変換した後に時間軸方向の前半部分又は後半部分の係数が０になるようにした後、第二拘束分割番号q'以外のものについてはそのまま上記フィルタ係数修正量バッファ１２２０の内容である周波数領域非拘束フィルタ係数修正量に上書きし、第一時間領域変換部１２０５において上記第一拘束分割番号qに対応する上記周波数領域非拘束フィルタ係数修正量ΔH'(q)を（DB2-DH'）により時間領域信号の非拘束フィルタ係数修正量Δh'(q)に変換し、第一要素選択部１２０６において上記非拘束フィルタ係数修正量のΔh'(q)の一部を（DB2-DH）により選択することでフィルタ係数修正量Δh(q)  とする。
【００６１】
（８）フィルタリング部１１ではフィルタリングを周波数領域で行っている場合は第三周波数変換部１２０７（図７）において上記フィルタ係数修正量Δh(q)を周波数領域フィルタ係数修正量ΔH(q)に（DB2-FDH）により変換する。
（９）第二平滑部１２１１において上記拘束分割番号qに対応する上記フィルタ係数h(q)に上記フィルタ係数修正量Δh(q)を（DB2-TD-TRANS）により加算するか、あるいは、第二平滑部１２１１において上記拘束分割番号qに対応する上記周波数領域フィルタ係数H(q)に上記周波数領域フィルタ係数修正量ΔH(q)を（DB2-TRANS）により加算するかし、上記第二平滑部１２１１における加算結果を上記FGフィルタ係数バッファ１２０９Ｆの内容に上書きする。
【００６２】
（１０）再初期値計算部１２３０において（DB2-HCLR）により再初期値を求め、その値で上記フィルタ係数修正量バッファ１２２０内の上記周波数領域非拘束フィルタ修正量ΔH'(q)を上書きする。
次に、第二の発明による伝達系推定方法の実施例を、図１６を（第二発明の方法）を用いて説明する。一般に伝達系推定方法は、データが蓄積するまで待つ条件とフィルタリングを行うステップと誤差信号eを求めるステップとフィルタ係数を更新するステップを有する。データが蓄積するまで待つ条件C1は参照信号と所望信号がNサンプル得られるまで待ち、参照信号xと所望信号dがNサンプル得られると、ステップS1において参照信号xを（FDBAP-FX）により周波数領域に変換することで周波数領域参照信号Xを求める。
ステップS2において上記周波数参照信号xを参照信号バッファに（FDBAP-FXBUF）により蓄積する。
フィルタリングを行うステップS3では時間領域で参照信号xとフィルタ係数hを（BAP-FLT）を用いてフィルタリングするか、あるいは、それに等価な処理を周波数領域で（FD-FLT）により行うことで推定信号yを求める。
誤差信号を求めるステップS4では上記推定信号yと所望信号dの差をとることで誤差信号をe＝d-yにより求める。フィルタ係数更新ステップでは上記フィルタ係数hが参照信号xと所望信号dの間の伝達系を模擬し上記誤差信号eのパワーが最小になるように上記フィルタ係数hあるいは周波数領域フィルタ係数Hを更新し、条件C1に戻る。
【００６３】
第二の発明の伝達系推定方法の実施例では、フィルタ係数更新ステップが以下のステップを有することを特徴とする。つまり、
（１）ステップS10において上記参照信号バッファの内容Uと周波数領域非拘束フィルタ係数修正量ΔH'の内積を（DB2-BGFLT1）により周波数毎にとることで周波数領域推定誤差信号ΔY'^を得る。
（２）ステップS5において（FDBAP-FE）により上記誤差信号eを周波数領域に変換することで周波数領域誤差信号Eを求める。
【００６４】
（３）ステップS30において（CLAIM1-FE）により上記周波数領域誤差信号Eから上記周波数領域推定誤差信号ΔY'^を減じることで周波数領域修正量更新用近似誤差ΔE'^を得る。
（４）ステップS31において上記周波数領域修正量更新用近似誤差ΔE'^と上記参照信号バッファの内容Uから周波数領域修正量更新量ΔΔH'を得る。周波数領域修正量更新量ΔΔH'は、適応アルゴリズムによってことなり、FDBAPの場合は、（FDBAR-R）を行った後に（CLAIM-FDBAP-GAIN）により、FDBRLSの場合は、Λ＝Iの場合のみ計算可能で、（FDBRLS-R）、（FDBRLS-GAIN1）と（FDBRLS-GAIN2）を行った後に（FDBRLS-dGAIN3）により求める。
【００６５】
（５）ステップS15において上記周波数領域修正量更新量ΔΔH'をフィルタ係数修正量バッファの内容ΔH'に（DB2-FDBAP-UPD）あるいは（DB2-FDBRLS-UPD）により重み付け加算し上記周波数領域非拘束フィルタ係数修正量ΔH'を更新し、上記周波数領域非拘束フィルタ係数修正量ΔH'を上記フィルタ係数修正量バッファに保存し、ステップS16において第二拘束分割番号q'を決定し、上記周波数領域非拘束フィルタ係数修正量ΔH'のうち第二拘束分割番号に対応する部分ΔH'を(FDBAP-DH')、(FDBAP-DH)と(FDBAP-FDH)により時間領域に変換した際に時間軸方向の前半部分又は後半部分の係数が０になるようにした後、上記フィルタ係数修正量バッファに保存し、
（６）上記条件C1とステップS1、S2、S3、S4、S10、S5、S30、S31、S15、S16と同期、あるいは、非同期で、条件C2において更新間隔が経過するのを待ち、更新間隔が経過したなら、以下のステップS17、S18、S19を行う。
【００６６】
（７）ステップS17において拘束分割番号qを決定し上記周波数領域非拘束フィルタ係数修正量ΔH'のうち上記拘束分割番号qに対応する部分ΔH'(q)を（DB2-DH'）により時間領域の上記拘束分割番号qに対応する非拘束フィルタ係数修正量Δh'(q)に変換し、上記拘束分割番号qに対応する上記非拘束フィルタ係数補正量Δh'(q)の一部を（DB2-DH）により選択することで上記拘束分割番号qに対応するフィルタ係数修正量Δh(q)とする。拘束分割番号qの決定方法は、q＝(k/N)mod(l/N)と均等な機会で更新されるように決める、あるいは、事前にH(m)のどの分割ノルムが大きいか分かれば、ノルムの大きい分割ほどより頻繁にするなど不均一な機会で更新されるように決める。
【００６７】
（８）ステップS18において上記拘束分割番号qに対応する上記フィルタ係数修正量Δh(q)を上記FGフィルタ係数バッファ１２２０の内容に（DB2-TD-TRANS）により加算するか、あるいは、上記拘束分割番号qに対応する上記フィルタ係数修正量Δh(q)を（DB2-FDH）により周波数領域に変換し上記拘束分割番号qに対応する周波数領域フィルタ係数修正量ΔH(q)を得、上記拘束分割番号qに対応する上記周波数領域フィルタ係数修正量ΔH(q)を上記FGフィルタ係数バッファの内容に（DB2-TRANS）により加算する。
（９）ステップS19において更新したフィルタ係数h(q)あるいは周波数領域フィルタ係数H(q)の拘束分割番号qに対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量ΔH'を、（DB2-HCLR）により微少値にし、条件C2に戻る。
【実施例３】
【００６８】
次に、第三の発明による伝達系推定装置の実施例を図７と図９（発明３による装置）を用いて説明する。この実施例では第２発明による伝達系推定装置の実施例において以下の差異を有することを特徴とする。つまり、
フィルタ係数更新部１２の中の上記乗算器１２１２の出力を上記周波数領域差分器１２１６に入力することに替え、上記演算器１２１２の出力ΔY'^を近似拘束部１２４０に入力し、上記近似拘束部１２４０において（CLAIM2-CNSTRN）により近似拘束周波数領域推定誤差信号＃ΔY'^を求め、上記近似拘束周波数領域推定誤差信号＃ΔY'^を上記周波数領域差分器１２１６に入力することを特徴とする。
【００６９】
次に、第三発明による伝達系推定方法の実施例を図１７（発明３の方法）により説明する。この実施例では第二の発明による伝達系推定方法の実施例において以下の差異を有することを特徴とする。つまり、ステップS5の直後で、ステップS40を置きステップS40において上記周波数領域推定誤差信号ΔY'^を（CLAIM2-CNSTRN）により近似拘束処理し近似拘束周波数領域推定誤差信号＃ΔY'^を得、ステップS30をステップS41に替えステップS41において上記周波数領域推定誤差信号ΔEから上記近似拘束周波数領域推定誤差信号＃ΔY'^を（CLAIM2-FE）により減じることで周波数領域修正量更新用近似誤差を得ることを特徴とする。
【００７０】
以上は図７に示した間引き拘束修正量更新量加算部１２２を用いて第一発明、第二発明、第三発明に対応する伝達系推定装置１の実施例である。
以下では伝達系推定装置１を構成する間引き拘束修正更新量が加算部として図１０に示す間引き拘束修正量更新量加算部１２２を用いた第四発明、第五発明、第六発明の実施例を説明する。
【実施例４】
【００７１】
この実施例は図６に示す伝達系推定装置１に図１０に示す間引き拘束修正量更新量加算部１２２を適用した実施例を示す。図６に示すフィルタリング部１１において時間領域で参照信号xとフィルタ係数hを（BAP-FLT）によりフィルタリングするか、あるいは、それに等価な処理を周波数領域で（FD-FLT）により行うことで推定信号yを求め、差分器１３では上記推定信号yと所望信号dの差をとることで誤差信号eをe＝d-yにより求め、フィルタ係数更新部１２においては上記フィルタ係数Hが参照信号xと所望信号dの間の伝達系を模擬し上記誤差信号eのパワーが最小となるように上記フィルタ係数hあるいは周波数領域フィルタ係数Hを更新する。
【００７２】
上記フィルタ係数更新部１２において、第一周波数領域変換部１２０１において参照信号x(k)を（FDBAP-FX）により、周波数領域の周波数領域参照信号X2N(k)に変換し、参照信号バッファ１２０２において現時刻の上記周波数領域参照信号Xを蓄積することで参照信号バッファ１２０２の内容Uを（FDBAP-FXBUF）により更新し、乗算器１２１２で上記参照信号バッファの内容Uとフィルタ係数修正量バッファ１２２０（図１０）の内容である周波数領域非拘束フィルタ係数修正量ΔH'(k)の積を（DB2-BGFLT1）により求め、その結果として周波数領域推定誤差信号ΔY'^(k)を生成する。
【００７３】
第二時間領域変換部１２１３において上記周波数領域推定誤差信号ΔY'^(k)を時間領域の選択前推定誤差信号Δy'^(k)に（DB2-BGFLT2）により変換し、第二要素選択部１２１４において上記選択前推定信号Δy'^(k)の一部を（DB-BGFLT3）により選択し推定誤差信号Δy^(k)とし、差分器１２１５で上記推定誤差信号Δy^(k)と誤差信号e(k)の差である修正量更新用誤差Δe^(k)を（DB2-ERR）により求め、第二周波数領域変換部１２０３において上記修正量更新用誤差Δe^(k)をFDBAPの場合は、（DB2-FDBAP-FE）により、FDBRLSの場合は、（DB2-FDBRLS-FE）により周波数領域信号の周波数領域修正量更新量用誤差ΔE^(k)に変換する。
【００７４】
更新量計算部１２０４において上記周波数領域修正用更新用誤差ΔE^(k)と上記参照信号バッファ１２０２の内容U2L,2N(k)から、利用する適応フィルタ手法に従って周波数領域修正量更新量ΔΔH'(k)を求め、間引き拘束修正量更新量加算部１２２において、第一平滑部１２０８（図１０）において上記周波数領域非拘束フィルタ係数修正量ΔH'(k)に上記周波数領域修正量更新量ΔΔH'(k)を利用する適応フィルタ手法に従って（DB2-FDBAP-UPD）あるいは（DB2-FDBRLS-UPD）により加算し、拘束部分決定部１２１０（図１０）で拘束分割番号qを決定する。
【００７５】
上記第一平滑部１２０８（図１０）の出力のうち上記拘束分割番号q以外は上記フィルタ係数修正量バッファ１２２０内の上記周波数領域非拘束フィルタ係数修正量ΔH'(k)に上書きする。第一時間領域変換部１２０５において上記第一平滑部１２０８の出力のうち上記拘束分割番号qに対応するΔH'(q)(k)を（DB2-DH'）により非拘束フィルタ係数修正量Δh'(q)に変換し、第一要素選択部１２０６において上記非拘束フィルタ係数修正量Δh'(q)を（DB2-DH）によりフィルタ係数修正量Δh(q)に変換し第二平滑部１２１１（図６）において上記FGフィルタ係数バッファ１２０９Ｆの内容に（DB2-TD-TRANS）により加算しFGフィルタ係数バッファ１２０９Ｆに保存するか、あるいは、第三周波数変換領域変換部１２０７において上記フィルタ係数修正量Δh(q)を（DB2-FDH）により周波数領域フィルタ修正量ΔH(q)に変換し、上記周波数領域非拘束フィルタ係数修正量ΔH'(k)と周波数領域フィルタ修正量ΔH(q)を間引き拘束修正量更新量加算部１２２の出力とし、第二平滑部１２１１において上記周波数領域フィルタ修正量ΔH(q)を周波数領域フィルタ係数H(q)に（DB2-TRANS）により加算しFGフィルタ係数バッファ１２０９Ｆに保存する。さらに再初期値計算部１２３０において（DB2-HCLR）により微少値の再初期値を求め、その値で上記フィルタ係数修正量バッファ１２２０内の上記周波数領域非拘束フィルタ修正量ΔH'(q)を上書きする。
【００７６】
図１８（発明四の方法）は、修正量蓄積型DBPFDLMS法による方法を説明する図であって、条件C1で参照信号と所望信号がNサンプル得られるまで待ち、参照信号と所望信号がNサンプル得られると、
ステップS1において（FDBAP-FX）等により参照信号を周波数領域に変換することで周波数領域参照信号xを求める。
ステップS2において上記周波数領域参照信号xを参照信号バッファに（FDBAP-FXBUF）等により蓄積する。
【００７７】
ステップS3において（BAP-FLT）や（FD-FLT）等によりフィルタ係数バッファにあるフィルタ係数hあるいは周波数領域フィルタ係数Hと参照信号xあるいは参照信号バッファにある周波数領域参照信号から推定信号yを求める。
ステップS4において（BAP-ERR）等により上記推定信号yと所望信号dの差をとることで誤差信号eを求める。
ステップS10においてフィルタ係数修正量バッファにある周波数領域非拘束フィルタ係数修正量ΔH'(k)と参照信号バッファの内容Uから（DB2-BGFLT1）により周波数領域推定誤差信号ΔY'^を求める。
【００７８】
ステップS11において上記周波数領域推定誤差信号ΔY'^を（DB2-BGFLT2）と（DB2-BGFLT3）により推定誤差信号Δy^を求める。
ステップS12において上記推定誤差信号Δy^と誤差信号eの差として（DB2-ERR）により修正用誤差信号Δe^を求める。
ステップS13において上記修正用誤差信号Δe^を（DB2-FE）により周波数領域に変換し周波数領域修正用誤差ΔE^を求める。
ステップS14において上記周波数領域修正用誤差ΔE^と参照信号バッファの内容Uから（DB2-FDBAP-GAIN）等により周波数領域修正量更新量ΔΔH'を求める。
【００７９】
ステップS15において上記周波数領域修正量更新量ΔΔH'をフィルタ係数修正量バッファの内容ΔH'に（DB2-FDBAP-UPD）等により重み付け加算し周波数領域非拘束フィルタ係数修正量ΔH'を得て上記周波数領域非拘束フィルタ係数修正量ΔH'を上記フィルタ係数修正量バッファに保存し、条件C1に戻るループと、上記条件C1とステップS1、S2、S3、S4、S10、S11、S12、S13、S14、S15と同期あるいは非同期に条件C2で更新間隔が経過するのを待ち、更新間隔が経過したなら、
ステップS17において拘束分割番号qを決定し上記周波数領域非拘束フィルタ係数修正量ΔH'のうち上記拘束分割番号qに対応する部分ΔH'(q)を（BD2-DH'）により時間領域の非拘束フィルタ係数修正量Δh'(q)に変換し、上記非拘束フィルタ係数修正量Δh'(q)の一部を（DB2-DH）により選択することでフィルタ係数修正量Δh(q)とする。
【００８０】
ステップS18において上記フィルタ係数修正量Δh(q)を上記FGフィルタ係数バッファの内容に（DB2-TD-TRANS）により加算するか、あるいは、上記フィルタ係数修正量を（DB2-FDH）により周波数領域に変換し周波数領域フィルタ係数修正量ΔH(q)を得、上記周波数領域フィルタ係数修正量ΔH(q)を上記FGフィルタ係数バッファの内容に（DB2-TRANS）により加算する。
ステップS19において更新したフィルタ係数hあるいは周波数領域フィルタ係数Hの非拘束分割番号qに対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量ΔH'(q)(k)を（DB2-HCLR）により最小値にし、条件C2に戻るループからなる。
第五発明と第六発明の概要は上述した第二発明と第三発明の概要と同様であるから、ここではその説明を省略する。
【実施例５】
【００８１】
第五の発明のよる伝達系推定装置の実施例を、図８と図１０（発明五による装置）を用いて説明する。一般に、伝達系推定装置は、フィルタリング部１１と差分器１３とフィルタ係数更新部１２を有する。フィルタリング部１１において時間領域で参照信号xとフィルタ係数hを（BAP-FLT）によりフィルタリングするか、あるいは、それに等価な処理を周波数領域で（FD-FLT）により行うことで推定信号yを求め、差分器１３では上記推定信号yと所望信号dの差をとることで誤差信号eをe＝d-yにより求め、フィルタ係数更新部１２においては上記フィルタ係数Hが参照信号xと所望信号dの間の伝達系を模擬し上記誤差信号eのパワーが最小となるように上記フィルタ係数hあるいは周波数領域フィルタ係数Hを更新する。
【００８２】
第五の発明の伝達系推定装置１の実施例では、上記フィルタ係数更新部１２が、以下の処理を有することを特徴とする。つまり、
（１）第一周波数領域変換部１２０１において上記参照信号xを（FDBAP-FX）により周波数領域の周波数領域参照信号Xに変換する。ただし、周波数領域への変換は、従来の技術で述べたように、DFT、DCT、DSTやDHTなどの直交変換を用いることもできる。
（２）参照信号バッファ１２０２において現時刻の上記周波数領域参照信号Xを（FDBAP-FXBUF）により蓄積することで参照信号バッファ１２０２の内容Uを更新する。
【００８３】
（３）乗算器１２１２で上記参照信号バッファ１２０２の内容Uとフィルタ係数修正量バッファ１２２０（図１０）の内容である周波数領域非拘束フィルタ係数修正量ΔH'の内積を（DB2-BGFLT1）により周波数毎に取りその結果として周波数領域推定誤差信号ΔY'^を生成する。
（４）第二周波数領域変換部１２０３において（FDBAP-FE）により上記誤差信号eを周波数領域信号の周波数領域誤差Eに変換する。
（５）周波数領域差分器１２３０において上記周波数領域誤差信号Eと上記周波数領域推定誤差信号ΔY'^との差を（CLAIM1-FE）によりとることで周波数領域修正量更新用近似誤差ΔE'^を求める。
【００８４】
（６）更新量計算部１２０４において上記周波数領域修正量更新用近似誤差ΔE'^と上記参照信号バッファの内容Uから周波数領域修正量更新量ΔΔH'を求める。周波数領域修正量更新量ΔΔH'は、適応アルゴリズムによって異なり、FDBAPの場合は、（FDBAR-R）を行った後に（CLAIM-FDBAP-GAIN）により、FDBRLSの場合は、Λ＝Iの場合のみ計算可能で、（FDBRLS-R）、（FDBRLS-GAIN1）と（FDBRLS-GAIN2）を行った後に（FDBRLS-dGAIN3）により求める。
（７）間引き拘束修正量更新量加算部１２２において上記非拘束フィルタ係数修正量ΔH'に上記周波数修正量更新量ΔΔH'を以下のように加えることで更新するとともに周波数領域フィルタ係数修正量ΔHを得る。第一平滑部１２０８（図１０）において上記周波数領域非拘束フィルタ係数修正量ΔH'に上記周波数領域修正量更新量ΔΔH'を(DB2-FDBAP-UPD)あるいは(DB2-FDBRLS-UPD)により重み付け加算する。拘束部分決定部１２１０において拘束分割番号qを決定する。qの決定方法は、q＝(k/N)mod(L/N)と均等な機会で更新されるように決める、あるいは、事前にH(m)のどの分割ノルムが大きいかが分かれば、ノルムの大きい分割ほどより頻繁にするなど不均一な機会で更新されるように決める。上記第一平滑部１２０８における重み付け加算結果のうち拘束分割番号qに相当しない上記周波数領域非拘束フィルタ係数修正量ΔH'の一部あるいは全部分を上記フィルタ係数修正量バッファ１２２０の内容である周波数領域非拘束フィルタ係数修正量に上書きし、第一時間領域変換部１２０５において上記拘束分割番号qに対応する上記周波数領域非拘束フィルタ係数修正量ΔH'(q)を（DB2-DH'）により時間領域信号の非拘束フィルタ係数修正量Δh'(q)に変換し、第一要素選択部１２０６において上記非拘束フィルタ係数修正量のΔh'(q)の一部を（DB2-DH）により選択することでフィルタ係数修正量Δh(q)  とする。
【００８５】
（８）フィルタリング部１１ではフィルタリングを周波数領域で行っている場合は第三周波数変換部１２０７（図１０）において上記フィルタ係数修正量Δh(q)を周波数領域フィルタ係数修正量ΔH(q)に（DB2-FDH）により変換する。
（９）第二平滑部１２１１において上記拘束分割番号qに対応する上記フィルタ係数h(q)に上記フィルタ係数修正量Δh(q)を（DB2-TD-TRANS）により加算するか、あるいは、第二平滑部１２１１において上記拘束分割番号qに対応する上記周波数領域フィルタ係数H(q)に上記周波数領域フィルタ係数修正量ΔH(q)を（DB2-TRANS）により加算するかし、上記第二平滑部１２１１における加算結果を上記FGフィルタ係数バッファ１２０９Ｆの内容に上書きする。
【００８６】
（１０）再初期値計算部１２３０において（DB2-HCLR）により再初期値を求め、その値で上記フィルタ係数修正量バッファ１２２０内の上記周波数領域非拘束フィルタ修正量ΔH'(q)を上書きする。
次に、第五の発明による伝達系推定方法の実施例を、図１９を（第五発明の方法）を用いて説明する。一般に伝達系推定方法は、データが蓄積するまで待つ条件とフィルタリングを行うステップと誤差信号eを求めるステップとフィルタ係数を更新するステップを有する。データが蓄積するまで待つ条件C1は参照信号と所望信号がNサンプル得られるまで待ち、参照信号xと所望信号dがNサンプル得られると、ステップS1において参照信号xを（FDBAP-FX）により周波数領域に変換することで周波数領域参照信号Xを求める。
【００８７】
ステップS2において上記周波数参照信号xを参照信号バッファに（FDBAP-FXBUF）により蓄積する。
フィルタリングを行うステップS3では時間領域で参照信号xとフィルタ係数hを（BAP-FLT）を用いてフィルタリングするか、あるいは、それに等価な処理を周波数領域で（FD-FLT）により行うことで推定信号yを求める。
誤差信号を求めるステップS4では上記推定信号yと所望信号dの差をとることで誤差信号をe＝d-yにより求める。フィルタ係数更新ステップでは上記フィルタ係数hが参照信号xと所望信号dの間の伝達系を模擬し上記誤差信号eのパワーが最小になるように上記フィルタ係数hあるいは周波数領域フィルタ係数Hを更新し、条件C1に戻る。
【００８８】
第五発明の伝達系推定方法の実施例では、フィルタ係数更新ステップが以下のステップを有することを特徴とする。つまり、
（１）ステップS10において上記参照信号バッファの内容Uと周波数領域非拘束フィルタ係数修正量ΔH'の内積を（DB2-BGFLT1）により周波数毎にとることで周波数領域推定誤差信号ΔY'^を得る。
（２）ステップS5において（FDBAP-FE）により上記誤差信号eを周波数領域に変換することで周波数領域誤差信号Eを求める。
（３）ステップS30において（CLAIM1-FE）により上記周波数領域誤差信号Eから上記周波数領域推定誤差信号ΔY'^を減じることで周波数領域修正量更新用近似誤差ΔE'^を得る。
（４）ステップS31において上記周波数領域修正量更新用近似誤差ΔE'^と上記参照信号バッファの内容Uから周波数領域修正量更新量ΔΔH'を得る。周波数領域修正量更新量ΔΔH'は、適応アルゴリズムによってことなり、FDBAPの場合は、（FDBAR-R）を行った後に（CLAIM-FDBAP-GAIN）により、FDBRLSの場合は、Λ＝Iの場合のみ計算可能で、（FDBRLS-R）、（FDBRLS-GAIN1）と（FDBRLS-GAIN2）を行った後に（FDBRLS-dGAIN3）により求める。
【００８９】
（５）ステップS15において上記周波数領域修正量更新量ΔΔH'をフィルタ係数修正量バッファの内容ΔH'に（DB2-FDBAP-UPD）あるいは（DB2-FDBRLS-UPD）により重み付け加算し上記周波数領域非拘束フィルタ係数修正量ΔH'を更新し、上記周波数領域非拘束フィルタ係数修正量ΔH'を上記フィルタ係数修正量バッファに保存し、条件C1に戻る。
（６）上記条件C1とステップS1、S2、S3、S4、S10、S5、S30、S31、S15と同期、あるいは、非同期で、条件C2において更新間隔が経過するのを待ち、更新間隔が経過したなら、以下のステップS17、S18、S19を行う。
（７）ステップS17において拘束分割番号qを決定し上記周波数領域非拘束フィルタ係数修正量ΔH'のうち上記拘束分割番号qに対応する部分ΔH'(q)を（DB2-DH'）により時間領域の上記拘束分割番号qに対応する非拘束フィルタ係数修正量Δh'(q)に変換し、上記拘束分割番号qに対応する上記非拘束フィルタ係数補正量Δh'(q)の一部を（DB2-DH）により選択することで上記拘束分割番号qに対応するフィルタ係数修正量Δh(q)とする。拘束分割番号qの決定方法は、q＝(k/N)mod(l/N)と均等な機会で更新されるように決める、あるいは、事前にH(m)のどの分割ノルムが大きいか分かれば、ノルムの大きい分割ほどより頻繁にするなど不均一な機会で更新されるように決める。
【００９０】
（８）ステップS18において上記拘束分割番号qに対応する上記フィルタ係数修正量Δh(q)を上記FGフィルタ係数バッファ１２２０の内容に（DB2-TD-TRANS）により加算するか、あるいは、上記拘束分割番号qに対応する上記フィルタ係数修正量Δh(q)を（DB2-FDH）により周波数領域に変換し上記拘束分割番号qに対応する周波数領域フィルタ係数修正量ΔH(q)を得、上記拘束分割番号qに対応する上記周波数領域フィルタ係数修正量ΔH(q)を上記FGフィルタ係数バッファの内容に（DB2-TRANS）により加算する。
（９）ステップS19において更新したフィルタ係数h(q)あるいは周波数領域フィルタ係数H(q)の拘束分割番号qに対応するフィルタ係数修正量バッファ中の周波数領域非拘束フィルタ係数修正量ΔH'を、（DB2-HCLR）により微少値にし、条件C2に戻る。
【実施例６】
【００９１】
次に、第六発明による伝達系推定装置の実施例を図９と図１０（第六発明による装置）を用いて説明する。この実施例では第五発明による伝達系推定装置の実施例において以下の差異を有することを特徴とする。つまり、
フィルタ係数更新部１２の中の上記乗算器１２１２の出力を上記周波数領域差分器１２１６に入力することに替え、上記演算器１２１２の出力ΔY'^を近似拘束部１２４０に入力し、上記近似拘束部１２４０において（CLAIM2-CNSTRN）により近似拘束周波数領域推定誤差信号＃ΔY'^を求め、上記近似拘束周波数領域推定誤差信号＃ΔY'^を上記周波数領域差分器１２１６に入力することを特徴とする。
【００９２】
次に、第六発明による伝達系推定方法の実施例を図２０（第六発明の方法）により説明する。この実施例では第五発明による伝達系推定方法の実施例において以下の差異を有することを特徴とする。つまり、ステップS5の直後で、ステップS40を置きステップS40において上記周波数領域推定誤差信号ΔY'^を（CLAIM2-CNSTRN）により近似拘束処理し近似拘束周波数領域推定誤差信号＃ΔY'^を得、ステップS30をステップS41に替えステップS41において上記周波数領域推定誤差信号ΔEから上記近似拘束周波数領域推定誤差信号＃ΔY'^を（CLAIM2-FE）により減じることで周波数領域修正量更新用近似誤差を得ることを特徴とする。
【００９３】
以下に、本発明を含む応用装置をいくつか説明するが、方法、プログラム及び記録媒体についても同様である。
（音響エコーキャンセラ）
図２６（音響EC）は本発明の第一の応用装置例を説明する図であって、
スピーカSPとマイクロホンMを用いたハンズフリー通信中にマイクロホンMで拾われるスピーカSPから放射された受話信号が相手側に戻るのを防止するための音響エコーキャンセラにおいて、受話信号を参照信号x、マイクロホンMで収音した音を所望信号dとして本発明による伝達系推定装置１を使用することを特徴とするエコーキャンセラ。
従来法より演算量の小さい本発明による伝達系推定装置１を用いることで、装置が小型化、廉価になる、適応フィルタ長を伸ばすことでエコー消去量が向上し、通話品質が向上するという利点がある。
【００９４】
（回線エコーキャンセラ）
図２７（回線EC）は本発明の第二の応用装置例を説明する図であって、通信回線の２線４線変換や回線の分岐によるエコーを消去するための回線エコーキャンセラにおいて、最寄りの２線４線変換器Hへ向かう信号を参照信号x、最寄りの２線４線変換器Hからの信号を所望信号dとして本発明による伝達系推定装置１を使用することを特徴とする回線エコーキャンセラ。
従来法より演算量の小さい本発明による伝達系推定装置を用いることで、装置が小型化、廉価になる、適応フィルタ長を伸ばすことでエコー消去量が向上し通信品質あるいはデータ転送速度が向上するという利点がある。
【００９５】
（適応等価器、教師信号有り１）
図２８（適応等価器、教師信号有り１）は本発明の第三の応用装置例を説明する図であって、データ伝送の受信側で通信路の伝達特性を補正するための適応等価器において、あらかじめ決められている学習用データあるいは判定器４の出力を参照信号x、受信信号を所望信号dとして本発明による伝達系推定装置１を使用することを特徴とする適応等価器。本発明により通信路の伝達特性を推定伝達特性フィルタとして求め、逆フィルタ計算部２において上記推定伝達特性フィルタの伝達特性逆フィルタを求め、等価フィルタリング部３において上記伝達特性逆フィルタで受信信号をフィルタリングすることで等価する。
従来法より演算量の小さい本装置を用いることで、装置が小型化、廉価になる、適応フィルタ長を伸ばすことで伝達経路の推定精度が向上し通話品質あるいはデータ転送速度が向上するという利点がある。
【００９６】
（適応等価器、教師信号有り２）
図２９（適応等価器、教師信号有り２）は本発明の第四の応用装置例を説明する図であって、データ伝送の受信側で通信路の伝達特性を補正するための適応等価器において、受信信号を参照信号x、あらかじめ決められている学習用データあるいは判定器４の出力を所望信号dとして本発明による伝達系推定装置１を使用することを特徴とする適応等価器。本発明の伝達系推定装置１により通信路の伝達特性の逆特性を推定伝達特性逆フィルタとして求め、等価フィルタリング部３において上記推定伝達特性逆フィルタで受信信号をフィルタリングすることで等価する。
従来法より演算量の小さい本発明装置１を用いることで、装置が小型化、廉価になる、適応フィルタ長を伸ばすことで伝達経路の推定精度が向上し通話品質あるいはデータ転送速度が向上するという利点がある。
【００９７】
（適応等価器、教師信号無し）
図３０（適応等価器、教師信号無し）は本発明の第五の応用装置例を説明する図であって、データ伝送の受信側で通信路の伝達特性を補正するための適応等価器において、受信信号と受信信号を遅延したものを参照信号x、受信信号あるいは判定器４の出力を所望信号dとして本発明を使用することを特徴とする適応等価器。本発明により通信路の伝達特性の逆特性を推定伝達特性逆フィルタとして求め、等価フィルタリング部３において上記推定伝達特性逆フィルタで受信信号をフィルタリングすることで等価する。
従来法より演算量の小さい本発明装置１を用いることで、装置が小型化、廉価になる、適応フィルタ長を伸ばすことで伝達経路の推定精度が向上し通話品質あるいはデータ転送速度が向上するという利点がある。
【００９８】
（ソナー、レーダー、医療用イメージング機器、非破壊検査機器）
図３１（ソナー・レーダー）は本発明の第六の応用装置例を説明する図であって、音波、電波や物質の振動などの入力波を送信器５により媒質に放射しその応答を受信器６で受信し上記本発明による伝達系推定装置１で伝達系の特性を推定し、推定した伝達系の特性から媒質中の波の伝播速度が変化する点の集まりを推定するソナー、レーダー、医療用イメージング機器、非破壊検査機器において、上記入力波を参照信号x、上記受信器６で受信した信号を所望信号dとして本発明装置１を使用することを特徴とするソナー、レーダー、医療用イメージング機器、非破壊検査機器。本発明装置１により上記媒質の伝達系を推定する。
従来法より演算量の小さい本発明装置１を用いることで、装置が小型化、廉価になる、あるいはこれまでと同一の処理装置規模で送信器や受信器の数を増加させ伝達速度が変化する点をより多く推定できソナー、レーダー、医療用イメージング機器、非破壊検査装置としての性能が向上するという利点がある。
【００９９】
（受信器アレー装置）
図３２（受信器アレー装置）は本発明の第七の応用装置例を説明する図であって、音波あるいは電波を複数の受信器６で受信し、目的の音波あるいは電波と目的以外の電波あるいは電波のレベル比を改善するために、目的の音波あるいは電波を受信する指向特性を有する強調アレー処理部７と目的の音波あるいは電波以外を受信する指向特性を有する遮断アレー処理部８と強調アレー処理部７の出力から遮断アレー処理部８の出力に含まれる目的の音波あるいは電波を除去する雑音キャンセラ部９を有する受信器アレー装置において、上記遮断アレー処理部８で強調アレー処理部７の出力を参照信号x、受信器６で受信した信号を所望信号dとして本発明装置１を利用し、上記雑音キャンセラ部９において遮断アレー処理部８で上記遮断アレー処理部８の出力を参照信号x、強調アレー処理部７の出力を所望信号dとして本発明装置１を利用することを特徴とする受信器アレー装置。
従来法より演算量の小さい本発明装置１を用いることで、装置が小型化、廉価になる、あるいは、従来と同じ演算量で受信器数を増やすことで遮断アレー処理装置の出力に混入する所望信号が減り雑音キャンセラの出力で所望信号レベルが増加するという利点がある。
【０１００】
本発明と従来法の比較
従来法と第二、第三及び第五，第六の発明について、演算量と収束特性（フィルタ係数誤差の収束速度、定常状態誤差レベル）を比較する。
○演算量
従来での技術と同様に、ブロック間隔あたりのFFT、Inverse FFT（IFFT）と周波数領域での複素積の回数を比較する。PFDLMS、DBPDFLMSと本発明における周波数領域フィルタ係数の転送間隔をLとしたときの値は表２のとおりである。
【０１０１】
表２
FFT,IFFT        複素積
PFDLMS              ３＋２        ４Ｌ＋４Ｎ
DBFDLMS             ３＋２＋１    ６Ｌ＋４Ｎ
第一発明            ３＋２        ６Ｌ＋４Ｎ
第二発明            ３＋２        ６Ｌ＋（４＋（２Ｄ＋１））Ｎ
表２において、従来法と本発明を比較すると、本発明は従来法DBPFDLMSに比してFFTの回数が１回少ない。また、第三及び第六の発明で増加した（２Ｄ＋１）Ｎの複素積は近似拘束部１２４０での演算によるものであるが、N＜＜Ｌの場合には６Ｌに比べて十分小さい。
【０１０２】
○収束特性
収束速度の遅いUMDF以外は、どの従来手法も同じ収束特性であることから、従来手法としてはPFDLMSだけを取り上げ、本発明の手法と比較する。
比較のための試験条件は従来技術の従来法の比較と同じである。
試験の結果を図２４と図２５（従来手法と本発明の係数誤差収束特性の比較１、２）に示す。同図の横軸はサンプル時間を、縦軸は、時間領域の適応フィルタ係数と伝達係数のインパルス応答の差の二乗和を伝達系インパルス応答の二乗和で除した値であり、小さければフィルタ係数の誤差が少なくよい推定ができているといえる。
図２４に示す参照信号が白色ガウス雑音の場合、点線で示した第五番目及び第六番目の発明以外はほぼ同じ収束速度、定常状態誤差レベルであり、一番目の発明だけ収束速度が３割程度遅い。
図２１に示す参照信号が男性音声の場合も点線で示した一番目の発明だけ収束速度が若干遅い。
【０１０３】
上述した本発明の伝達系推定方法はコンピュータが解読可能なプログラム言語によって記述された伝達系推定プログラムをコンピュータにインストールし、このプログラムをコンピュータに備えられたCPUにより解読させて実行させることにより実現することができる。本発明による伝達系推定プログラムはコンピュータが読み取り可能な例えば磁気ディスクあるいはCD-ROM、半導体メモリ等の記録媒体に記録され、これらの記録媒体からコンピュータにインストールするか、又は通信回線を通じてインストールすることができる。
【産業上の利用可能性】
【０１０４】
本発明による伝達系推定装置、方法は図２６〜図３２に示したように音響エコーキャンセラの分野、適応等価器の分野、ソナー、レーダー、医療用イメージング機器の分野、非破壊検査機器の分野、受信器アレー等の分野に広く利用される。
【図面の簡単な説明】
【０１０５】
【図１】従来の時間領域適応フィルタの構成を説明するためのブロック図。
【図２】従来の周波数領域適応フィルタ手法であるFDBAP(UMDF)による装置を説明するためのブロック図。
【図３】従来の周波数領域適応フィルタ手法であるUMDFによる装置を説明するためのブロック図。
【図４】従来の周波数領域適応フィルタ手法であるPFDLMSによる装置を説明するためのブロック図。
【図５】従来の周波数領域適応フィルタ手法であるPFDLMSに従来法であるフィルタ係数二重化を適用したDBPFDLMSによる装置を説明するためのブロック図。
【図６】第一発明及び第四発明による装置の実施例を説明するためのブロック図。
【図７】第一発明、第二発明、第三発明に用いる間引き拘束修正量更新量加算部の構成を説明するためのブロック図。
【図８】第二発明及び第五発明による装置の実施例を説明するためのブロック図。
【図９】第三発明及び第六発明による装置の実施例を説明するためのブロック図。
【図１０】第四発明、第五発明、第六発明に用いる拘束修正量更新量加算部の構成を説明するためのブロック図。
【図１１】従来の周波数領域適応フィルタの手法であるFDBAP(UMDF)による方法を説明するためのフローチャート。
【図１２】従来の周波数領域適応フィルタの手法であるUMDFによる方法を説明するためのフローチャート。
【図１３】従来の周波数領域適応フィルタの手法であるPFDLMSによる方法を説明するためのフローチャート。
【図１４】従来の周波数領域適応フィルタの手法であるPFDLMSに従来法であるフィルタ係数二重化を適用したDBPFDLMSによる装置を説明するためのフローチャート。
【図１５】第一発明による方法の実施例を説明するためのフローチャート。
【図１６】第二発明による方法の実施例を説明するためのフローチャート。
【図１７】第三発明による方法の実施例を説明するためのフローチャート。
【図１８】第四発明による方法の実施例を説明するためのフローチャート。
【図１９】第五発明による方法の実施例を説明するためのフローチャート。
【図２０】第六発明による方法の実施例を説明するためのフローチャート。
【図２１】従来手法の白色ガウス信号に対する係数誤差収束特性を説明するためのグラフ。
【図２２】従来法の男性音声に対する係数誤差収束特性を説明するためのグラフ。
【図２３】従来手法のブロック境界での誤差信号の連続性を説明するためのグラフ。
【図２４】従来法と本発明法の白色ガウス信号に対する係数誤差収束特性を説明するためのグラフ。
【図２５】従来法と本発明法の男性音声に対する係数誤差収束特性を説明するためのグラフ。
【図２６】本発明を音響エコーキャンセラに応用した応用装置例を説明するためのブロック図。
【図２７】本発明を回線エコーキャンセラに応用した応用装置例を説明するためのブロック図。
【図２８】本発明を適応等価器、教師信号ありに応用した応用装置例を説明するためのブロック図。
【図２９】本発明を適応等価器、教師信号ありに応用した応用装置例を説明するためのブロック図。
【図３０】本発明を適応等価器、教師信号なしに応用した応用装置例を説明するためのブロック図。
【図３１】本発明をソナー、レーダー、医療用イメージング機器、非破壊検査機器に応用した応用装置例を説明するためのブロック図。
【図３２】本発明を受信器アレーに応用した応用装置例を説明するためのブロック図。
【符号の説明】
【０１０６】
１    伝達系推定装置                １２０７    第三周波数領域変換部
４    判定器                        １２０８    第一平滑部
１１    フィルタリング部              １２０９Ｆ  FGフィルタ係数バッファ
１１Ｂ  BGフィルタリング部            １２１０    第一拘束部分決定部
１２    フィルタ係数更新部            １２１０’  第二拘束部分決定部
１３    差分器                        １２１１    第二平滑部
１２０    更新量加算部                  １２１２    乗算器
１２１    間引き拘束修正量加算部        １２１３    第二時間領域変換部
１２２    間引き拘束修正量更新量加算部  １２１４    第二要素選択部
１２０１    第一周波数領域変換部          １２１５    差分器
１２０２    参照信号バッファ              １２１６    周波数領域差分器
１２０３    第二周波数領域変換部          １２２０    フィルタ係数修正量バッファ
１２０４    更新量計算部                  １２３０    再初期値計算部
１２０５    第一時間領域変換部            １２４０    近似拘束部
１２０６    第一要素選択部                １２５０    拘束部
【図１】
0004510795.tif 000002
【図２】
0004510795.tif 000003
【図３】
0004510795.tif 000004
【図４】
0004510795.tif 000005
【図５】
0004510795.tif 000006
【図６】
0004510795.tif 000007
【図７】
0004510795.tif 000008
【図８】
0004510795.tif 000009
【図９】
0004510795.tif 000010
【図１０】
0004510795.tif 000011
【図１１】
0004510795.tif 000012
【図１２】
0004510795.tif 000013
【図１３】
0004510795.tif 000014
【図１４】
0004510795.tif 000015
【図１５】
0004510795.tif 000016
【図１６】
0004510795.tif 000017
【図１７】
0004510795.tif 000018
【図１８】
0004510795.tif 000019
【図１９】
0004510795.tif 000020
【図２０】
0004510795.tif 000021
【図２１】
0004510795.tif 000022
【図２２】
0004510795.tif 000023
【図２３】
0004510795.tif 000024
【図２４】
0004510795.tif 000025
【図２５】
0004510795.tif 000026
【図２６】
0004510795.tif 000027
【図２７】
0004510795.tif 000028
【図２８】
0004510795.tif 000029
【図２９】
0004510795.tif 000030
【図３０】
0004510795.tif 000031
【図３１】
0004510795.tif 000032
【図３２】
0004510795.tif 000033