畳み込み

質問

ある回路の出力y(t)は，回路のインパルス応答がf(t)のとき，入力x(t)との畳み込みによって
y(t)=f(t)*x(t)=∫(-∞,∞)f(τ)(t-τ)dτ
です．（*は畳み込み演算）

さて，出力と応答特性から入力を求める場合は，どうするのでしょう？時系列データが等時間間隔でないため，fftが使えないので困っています．

回答No.6

noname#11476 · ベストアンサー

書き方にちょっと誤解があったかもしれませんので。
>やはりFFTが正統なのでしょうね
ガウス－ザイデルの方法は、FFTは使いません。
直接デコンボリューション（内部ではコンボリューションもしている）を行う手法です。

> #3にも書きましたとおり，補完してFFTでやってみました．
> 今のところ，期待した結果とはなっておりません．
単純にＦＦＴして、わり算して求めても、測定値が実測の場合はノイズの影響でうまく復元しないでしょう。

実測データに対してやらなければならない手順は、

1.データ補間
　No.5の方の言われる窓関数を使った方法はかなり有効です。
　詳しくはNumerical Recipes in C などの本を参照して下さい。

2.ノイズフィルタで可能な限りノイズの除去
　窓関数を使う方法をとれば、同時にフィルタの役目もしますので、これは必要なくなります。
　保管方法によってはこれがないとまともな結果が得られません。

3.デコンボリューション
　もっとも数学的にわかりやすいのはガウスザイデルの方法（ガウスの消去法で連立方程式を解く）ですね。
　これ以外にも沢山の方法が考案されています。

h(t)が推定とのことですが、それだけではそれらしい結果が得られない理由にはなりません。
うまくデコンボリューション出来ない理由の大半はノイズによる影響です。
このデコンボリューションというのはかなり難しく、色々なやり方が研究されています。
デコンボリューション(deconvolution)、ガウス、ウィナーフィルタなどでwebを検索すると色々出てくると思いますよ。

では、がんばって下さい。

回答No.5

nuubou

ｙの帯域に上限があることがはっきりしている場合には
直線補間より高度な
ｓｉｎ（ｘ）／ｘを使う方法があります
これだと無限個のサンプリング値を必要とするので
窓関数ｗ（ｘ）をかけて
ｗ（ｘ）・ｓｉｎ（ｘ）／ｘ
を補間に使うと良いでしょう
するとサンプリング値を補間点近辺に限定できます
サンプリング間隔が大きいと単なる直線補間は余りにも荒すぎます

回答No.4

noname#11476

>逆畳み込み積分は応答関数の複素共役と畳み込み積分をするということで良いのでしょうか？
どのようにイメージされているかが、今一つ理解できませんでしたが、

y(t) = f(t) * x(t)

の操作をしてy(t)を得ているという点はよろしいわけですね。
原理的には、No.1の方が書いたように、フーリエ変換で表すと、x(t)をy(t)とx(t)から逆算出来ることになります。
（これが最も簡単で早い解法です）

が、もし実測波形の場合は、
y(t) = f(t)*x(t) + n(t)
とノイズの項n(t)が入ってしまうために、最もx(t)に近い波形を推定しなければ求めることが出来ません。
ノイズn(t)も実測できれば良いのでしょうけど、事実上無理ですよね？
また、わずかなノイズでもx(t)に大きく影響してしまうためにうまく復元できません。

なお、どの手法を選ぶにしてもデータは等間隔でないと、扱いがとても面倒（私にはその扱い方は理解できませんでした）なので、普通は補間などの方法で、等間隔にマッピングし直して処理します。

>単に相関を取ったのではなぜだめなのでしょうか？
デコンボリューションは、応答関数f(t)によってゆがめられた信号から、元の波形を推定することになりますので、相関（コリレーションの意味ですか？）では解けないと思いますが。
それとも相関という意味を、最小２乗法でfittingするようなイメージで使われていますか？
ならば、ガウス-ザイデル法がそれに該当します。
y1(t) = f(t) * x1(t)
x1(t)を初期値として、y(t)-y1(t)の残差から、次のx2(t)を推定して、また残差をみて、、、と繰り返して、残差が小さくなるようなxn(t)を求めることで、元波形を推定します。

では。

お礼

またまた回答ありがとうございました．経過報告です．

本屋にて信号処理関係を漁りましたが，直接デコンボルーションをする方法については全く見つかりませんでした．やはりFFTが正統なのでしょうね．

#3にも書きましたとおり，補完してFFTでやってみました．
今のところ，期待した結果とはなっておりません．

元データが2万5000点を，補完して2万点にしてしまったので，この時点でおかしいかも知れません．もう少しやってみます．

又，f(t)またはf(ω)が正確に判っているわけでは無く，判っているのは時定数だけなのです．mHz程度の応答であるため，手持ちの測定器では調べることもままならずというところです．

このセンサ自体にヒステリシスも含まれているようなので，計測方法自体も見直すことにします．

回答No.3

nuubou

出力信号を直線補間をして等間隔にあるデータを求めｆｆｔを行えばいいと思います
もし精度を上げたかったら高次補間をすればいいと思います
ｘが未知なのだからｆは予測値で我慢するしかありません

お礼

補完処理->FFTしてみました．
f(t)は一次遅延で時定数のみが判っており，そのまま計算したのでは期待した結果は得られませんでした．

きっと周波数特性が一致していないのだと思います．他にもヒステリシスが含まれていそうで，これらを元信号からできる限り除去しないと，信号処理だけでは厳しいようでした．

助言ありがとうございました．

回答No.2

noname#11476

その操作は、一般に「逆畳み込み積分」または「デコンボリューション」と呼ばれます。
やり方は実に様々あるので、お探しになってみて下さい。一番簡単なのはガウス-ザイデルの方法です。
(連立一次方程式を解く形となります）

データが等時間間隔でない場合は、扱いが大変面倒なので、一度等時間間隔の空間へ補間などを用いてマッピングし直してから計算すると良いでしょう。

観測波形にノイズなどがあると、強調されてしまいますので、事前にフィルタリングが必要になったりします。
（単純なFFTによる復元では、このノイズのために実用にはなりません。）

では。

お礼

回答ありがとうございます．
まだよく判っていないのですが，逆畳み込み積分は応答関数の複素共役と畳み込み積分をするということで良いのでしょうか？

単に相関を取ったのではなぜだめなのでしょうか？

回答No.1

nuubou

Ｙ（ω）＝∫ｄｔ・ｅｘｐ（－ｉ・ω・ｔ）・ｙ（ｔ）
Ｘ（ω）＝∫ｄｔ・ｅｘｐ（－ｉ・ω・ｔ）・ｘ（ｔ）
Ｆ（ω）＝∫ｄｔ・ｅｘｐ（－ｉ・ω・ｔ）・ｆ（ｔ）
とする
Ｙ（ω）＝Ｆ（ω）・Ｘ（ω）
であるから
Ｘ（ω）＝Ｙ（ω）／Ｆ（ω）
である
これより
ｘ（ｔ）＝１／（２・π）・∫ｄω・ｅｘｐ（ｉ・ω・ｔ）・Ｘ（ω）
＝１／（２・π）・∫ｄω・ｅｘｐ（ｉ・ω・ｔ）・Ｙ（ω）／Ｆ（ω）
もしｙ（ｔ）が
ｙ（ｔ）＝Σ（ｎ）・ｘ［ｎ］・δ（ｔ－ｔ［ｎ］）
と表されるならば
Ｙ（ω）＝Σ（ｎ）・ｘ［ｎ］・ｅｘｐ（－ｉ・ω・ｔ［ｎ］）
である
ただしｉは虚数単位である（ｊにしたほうがいいと思いますが）

ｘ（ｔ），ｙ（ｔ），ｆ（ｔ）のうち何が時系列データになるのですか？
ａｄｃ、ｄａｃがどこかには行っているのでしょうか？
もっと詳しい構成を教えてください

お礼

回答ありがとうございます．
得られているのはtとy(t)だけで，f(t)というのは部品カタログからの予測値です．
y(t)の測定が等時間間隔で無い（ある条件から外れたときはデータを取り込まないようになっているため），fftは使えないのです．
tが判っているので，フーリエ積分をコツコツやればいいのかと思ったのですが，やはり等時間間隔で無いでない為，積分値に誤差が生じてしまいそうです．
そのためフーリエ積分を用いないでやればいいのかと思ったのでした．