1

我想在 MATLAB 中将包含两个尖峰（称为 Spike）的时间序列与指数内核（k）进行卷积。将卷积响应称为“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 询问他们的想法。

根据 MA​​TLAB 文档，如果[q,r] = deconv(v,u), thenv = conv(u,q)+r也必须成立（即 的输出deconv应始终满足这一点）。在您的情况下，这是严重违反的。将以下内容放在脚本的末尾：

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

我得到 6.75e227，它不为零；-) 接下来尝试将长度更改spike为 6000；你会得到一个小数字（~1e-15）。逐渐增加长度spike；误差会越来越大。请注意，如果您只将一个非零元素放入尖峰，则不会发生此行为：错误始终为零。这说得通; 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 向量相同的长度，它对我有用。IE：

    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 回答