为什么我尝试在 matlab 中验证卷积定理时得到一个很大的 MSE?

Why I got a big MSE when I try to verify convolution theorem in matlab?

我想在matlab中验证卷积定理

首先,我对二维高斯进行二维离散卷积 图像灰度图(x, y).

其次,我计算了傅里叶变换 相同的 2D 高斯和原始图像。然后执行标量乘法 这两个傅立叶变换,然后是结果的逆傅立叶变换。

最后,我将计算两个结果之间的均方误差。但是,我发现错误是800+。

这是我的代码:

[row, col] = size(graymap);
[row_2, col_2] = size(z);
result = zeros(row, col);
for i = 1: col
    for j = 1:row
    accumulation_value = 0;
    for k = -4:4
        for h = -4:4
            if ((i+k > 0 && i+k < col + 1) && (j+h > 0 && j+h < row + 1))
                value_image = double(graymap(i+k, j+h));
            else
                value_image = 0;
            end
            accumulation_value = accumulation_value + value_image * double(z(5 + k, 5 + h));
            weighted_sum = weighted_sum + z(5 + k, 5 + h);
        end
    end
    result(i,j) = (accumulation_value);
end
result_blur_1 = uint8(255*mat2gray(result));

M = size(graymap,1);
N = size(graymap,2);

resIFFT = ifft2(fft2(double(graymap), M, N) .* fft2(double(z), M, N));
result_blur_2 = uint8(255*mat2gray(resIFFT));
err = immse(result_blur_1, result_blur_2);

z是9*9的高斯核。我不翻转它,因为它是对称的。

我认为我的卷积实现是正确的,因为结果与 conv2(graymap, z, 'same').

因此,我认为第二部分有问题。事实上,我对填充的工作原理感到困惑。可能是大MSE的原因吧。

你的第二部分实现确实有问题。通过 fft 实现卷积时要记住的最重要的规则是,您实际上是在计算 circular 卷积,而不是 线性卷积。幸运的是,有一个条件可以让两者变得等价。这个条件是 两个数组应该被零填充以使其大小等于每个数组的大小之和减去 1(在所有维度上)。因此,如果您正在使用大小为 MxN 的图像 X 和大小为 PxQ 的遮罩 Z,则您应该用零填充两个数组,以便它们至少具有 个维度 M+P-1xN+Q-1。任何额外的零都不会造成伤害,因此如果可能的话匹配 'fft-friendly' 大小会很方便(例如使用 nextpow2)。您只需取前 M+P-1xN+Q-1 个值。

现在,如果您只想要卷积的完整结果,那将直接起作用。但是因为你想要卷积的中心部分(选项 'same'),你需要 select 正确的索引。第一个索引将是 ceil(([P Q] - 1)/2) + 1,然后你取与图像大小一样多的连续索引。

这是一个将所有内容放在一起的示例:

M = randperm(1024,1);
N = randperm(1024,1);

X = rand(M,N);

P = randperm(64,1);
Q = randperm(64,1);

Z = rand(P,Q);

% 'standard' convolution with option 'same'
C1 = conv2(X,Z,'same');


R = 2^nextpow2(M+P-1);
S = 2^nextpow2(N+Q-1);

% convolution with fft. Notice the zero-padding to R,S
C2 = real(ifft2(fft2(X,R,S) .* fft2(Z,R,S)));

n       = ceil(([P Q] - 1)/2);
ind{1}  = n(1) + (1:M);
ind{2}  = n(2) + (1:N);

C2 = C2(ind{:});

err = immse(C1,C2)

我得到 1e-26

顺序的错误