为什么我尝试在 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(在所有维度上)。因此,如果您正在使用大小为 M
xN
的图像 X
和大小为 P
xQ
的遮罩 Z
,则您应该用零填充两个数组,以便它们至少具有 个维度 M+P-1
xN+Q-1
。任何额外的零都不会造成伤害,因此如果可能的话匹配 'fft-friendly' 大小会很方便(例如使用 nextpow2
)。您只需取前 M+P-1
xN+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
顺序的错误
我想在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(在所有维度上)。因此,如果您正在使用大小为 M
xN
的图像 X
和大小为 P
xQ
的遮罩 Z
,则您应该用零填充两个数组,以便它们至少具有 个维度 M+P-1
xN+Q-1
。任何额外的零都不会造成伤害,因此如果可能的话匹配 'fft-friendly' 大小会很方便(例如使用 nextpow2
)。您只需取前 M+P-1
xN+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