1

MATLAB で指数カーネル (k) を使用して、2 つのスパイク (Spike と呼びます) を含む時系列を畳み込みたいと思います。畳み込み応答を「calcium1」と呼びます。カーネルでデコンボリューションを使用して、元のスパイク (「reconSpike」) データを復元したいと考えています。次のコードを使用しています。

k1=zeros(1,5000);
k1(1:1000)=(1.1.^((1:1000)/100)-(1.1^0.01))/((1.1^10)-1.1^0.01);
k1(1001:5000)=exp(-((1001:5000)-1001)/1000);
k1(1)=k1(2);

spike = zeros(100000,1);
spike(1000)=1;
spike(1100)=1;

calcium1=conv(k1, spike);
reconSpike1=deconv(calcium1, k1);

問題は、reconSpike の最後に、元のデータにはなかった非常に大きく高振幅の波の塊が得られることです。誰もがそれを修正する理由と方法を知っていますか?

ありがとう！

4

3 に答える 3

1

デコンボリューションが単純にコンボリューションを元に戻すことができるとは期待しないでください。これは、デコンボリューションが不適切な問題であるためです。

int f(x) g(x-t) dx問題は、畳み込みが整数演算子であるという事実から生じます (連続の場合、積分または類似のものを書き留めます)。ここで、積分の計算の逆 (デコンボリューション) は、微分を適用することです。残念ながら、差動は入力のノイズを増幅します。したがって、積分にわずかな誤差しかない場合 (浮動小数点の不正確さで十分な場合もあります)、微分後にはまったく異なる結果になります。

この増幅を軽減する方法はいくつかありますが、これらはアプリケーションごとに試す必要があります。

于 2012-08-23T16:58:46.537 に答える
1

MATLAB のデコンボリューション アルゴリズムの問​​題か、浮動小数点の精度の問題 (あるいはその両方) が発生しています。デコンボリューション中に行われるすべての除算と減算による浮動小数点の精度だと思いますが、MathWorks に直接連絡して、彼らの考えを尋ねる価値があるかもしれません。

MATLAB ドキュメンテーションによると、 if[q,r] = deconv(v,u)もv = conv(u,q)+r成立する必要があります (つまり、 の出力はdeconv常にこれを満たす必要があります)。あなたの場合、これは激しく違反しています。スクリプトの最後に次を追加します。

[reconSpike1 rem]=deconv(calcium1, k1);
max(conv(k1, reconSpike1) + rem - calcium1)

0 ではない 6.75e227 が得られます ;-) 次に、 の長さspikeを 6000 に変更してみてください。小さな数 (~1e-15) が得られます。spike;の長さを徐々に増やします。誤差はますます大きくなります。スパイクにゼロ以外の要素を 1 つだけ入れた場合、この動作は発生しないことに注意してください。エラーは常にゼロです。それは理にかなっている; MATLAB が行う必要があるのは、すべてを同じ数で除算することだけです。

ランダムなベクトルを使用した簡単なデモを次に示します。

v = random('uniform', 1,2,100,1);
u = random('uniform', 1,2,100,1);
[q r] = deconv(v,u);
fprintf('maximum error for length(v) = 100 is %f\n', max(conv(u, q) + r - v))
v = random('uniform', 1,2,1000,1);
[q r] = deconv(v,u);
fprintf('maximum error for length(v) = 1000 is %f\n', max(conv(u, q) + r - v))

出力は次のとおりです。

maximum error for length(v) = 100 is 0.000000
maximum error for length(v) = 1000 is 14.910770

あなたが本当に何を達成しようとしているのかわからないので、これ以上のアドバイスをするのは難しい. ただし、パルスが積み重なっている問題があり、各パルスに関する情報を抽出したい場合、これは難しい問題になる可能性があることを指摘しておきます。このようなことに取り組んでいる人を何人か知っているので、参照が必要な場合はお知らせください。

于 2012-08-22T21:42:03.270 に答える
1

スパイク ベクトルを k1 ベクトルと同じ長さに保つとうまくいきます。すなわち:

    k1=zeros(1,5000);
    k1(1:1000)=(1.1.^((1:1000)/100)-(1.1^0.01))/((1.1^10)-1.1^0.01);
    k1(1001:5000)=exp(-((1001:5000)-1001)/1000);
    k1(1)=k1(2);

    spike = zeros(5000, 1);
    spike(1000)=1;
    spike(1100)=1;

    calcium1=conv(k1, spike);
    reconSpike1=deconv(calcium1, k1);

違うものにした理由はありますか？

于 2012-08-22T14:28:09.343 に答える