在 MATLAB 中使用 3 种不同方法的连续卷积给出不同的结果

Question

我试图理解为什么当我尝试不同的方法时，MATLAB 中两个函数的 convolution 的结果是不同的。例如，假设我的函数是 sin(x) 和 cos(x)。第一种方法是使用 MATLAB 中的 conv() 命令。第二种方法是直接使用卷积的定义来计算它。第三种方法是使用卷积定理并使用 fft() 和 ifft() 计算它。代码附在此处。我使用 sin(x) 和 exp(-x) 作为示例。我想对任意两个任意函数执行相同的操作。

x = linspace(0,10,100);
dx = x(2)-x(1);
t = -100:0.01:100;

y = conv(sin(x),exp(-x),'same')*dx;
figure(1)
plot(x,y);

for i=1:length(x)
    zz = sin(t).*exp(-x(i)-t);
    zz = trapz(t,zz);
    z(i) = zz;
end
figure(2)
plot(x,z)

ww = fft(sin(x)).*fft(exp(-x));
w = ifft(ww).*dx;
figure(3)
plot(x,w)

如您所见，两个输入的大小相等。这三种方法的结果非常不同。我应该怎么做才能让 3 种方法给出相同的答案？如果有人在这里解释问题，我将不胜感激。

Answer 1

让我们看看OP中实现的三个方法：

x = linspace(0,10,100);
dx = x(2)-x(1);
t = -100:0.01:100;

y = conv(sin(x),exp(-x),'same')*dx;

此处您正在卷积的信号定义在 [0,10] 范围内，并假定在此范围之外为零。

for i=1:length(x)
    zz = sin(t).*exp(-x(i)-t);
    zz = trapz(t,zz);
    z(i) = zz;
end

此处信号定义在正弦波的 [-100,100] 范围内，指数的范围更大（实际范围发生变化，因此不会对范围外发生的情况做出任何假设） : [-110,100].

请注意，指数函数呈指数增长...在所用范围的下端，其值为 2.6881e+43，该值使其他一切相形见绌并决定了结果。

还有一个错误：exp(-x(i)-t) 应该读作 exp(-(x(i)-t))。

ww = fft(sin(x)).*fft(exp(-x));
w = ifft(ww).*dx;

这里的信号与第一种情况一样在 [0,10] 范围内定义，但假设在此范围之外重复。这相当于将 [sin(x),sin(x),sin(x)] 与 [exp(-x),exp(-x),exp(-x)] 进行卷积，而第一种情况是将 [zeros(size(x)),sin(x),zeros(size(x))] 与 [zeros(size(x)),exp(-x),zeros(size(x))] 进行卷积。

您可以像这样使第一个和最后一个案例相等：

a = [zeros(size(x)),sin(x),zeros(size(x))];
b = [zeros(size(x)),exp(-x),zeros(size(x))];
y = conv(a,b,'same');
w = ifft(fft(a).*fft(ifftshift(b)));
plot(y)
hold on
plot(w)

此处使用 ifftshift 将信号 b 的原点移动到 FFT 期望的第一个 bin。否则，结果将偏移信号长度的一半。如果你以同样的方式移动a，w必须以相反的方式移动，所以不移动它更简单。

这些都是离散时间卷积。采样理论说，对于两个带限信号，先卷积后采样与先采样后卷积是一样的，采样信号的插值可以return我们的连续信号。但只有当我们可以对函数进行采样直到无穷大时，这才是正确的，而我们做不到。另请注意，指数函数不受频带限制。所以这里把离散卷积和连续卷积联系起来有两个问题：采样域和函数的选择。

因此，让我们取两个信号，使这些比较有意义。对于第一个 a，我们可以再次采用正弦波。对于第二个，b，我们将采用高斯导数，它几乎是带限的，并且在时间上也是有限的。

x = linspace(-10, 10, 101); % odd number of samples so we sample the origin
dx = x(2) - x(1);

a = @(x) sin(x);
s = 3*dx;
b = @(x) -exp(-x.^2/(2*s^2)) / (sqrt(2*pi)*s^3) .* x;

y = conv(a(x), b(x), 'same') * dx;

t = linspace(-15, 15, 400); % some larger domain, and higher sampling rate, though this is not necessary!
z = zeros(size(x));
for i = 1:length(t)
    z(i) = trapz(t, a(t) .* b(t(i)-t));
end

w = ifft( fft(a(x)) .* fft(ifftshift(b(x))) ) * dx;

plot(x, y, 'DisplayName','conv()')
hold on
plot(t, z, 'DisplayName','"continuous"')
plot(x, w, 'DisplayName','FFT')
legend('-DynamicLegend');

请注意三种形式的卷积的输出在图中间的区域中几乎相同。在采样域的边缘附近，它们开始不同（该区域的大小由 b 的大小给出，这就是我们为其选择紧凑函数的原因）。

[另请注意，z实际上只是具有更高采样率的离散卷积，真的，考虑到 trapz 与 sum 相同，除了两个样本在末端减半。]

在 MATLAB 中使用 3 种不同方法的连续卷积给出不同的结果

Continuous Convolution Using 3 Different Methods in MATLAB Gives Different Results

matlab

math

signal-processing

fft

convolution