二维矩阵数据的 MATLAB 加窗 FFT（图像）

Question

数据：假设我有一个 2000 行 x 500 列的矩阵（图像）

我需要什么：计算 64 行乘以 10 列块的上述数据的 FFT。换句话说，我想计算在整个数据矩阵中运行的 64X10 窗口的 F​​FT。 FFT 结果用于计算标量值（例如峰值幅度频率），用于创建新的“FFT 值”图像。

现在，我需要最终的 FFT 图像与原始数据的大小相同 (2000 X 500)。

在 MATLAB 中完成此任务的最快方法是什么？我目前正在使用相对较慢的 for 循环。我还使用插值将最终图像的大小调整为原始数据大小。

Answer 1

正如@EitanT 指出的那样，您可以使用blockproc 对图像进行批量块处理J. 但是您应该将函数句柄定义为

fun = @(block_struct) fft2(block_struct.data);
B = blockproc(J, [64 10], fun);

对于[2000 x 500] 矩阵，这将为您提供复杂傅立叶值的[2000 x 500] 输出，在局部支持（FFT 的输入大小）为[64 x 10] 的子采样像素位置进行评估。现在，用单个值替换这些值，例如使用峰值对数幅度，您可以进一步指定

fun = @(block_struct) max(max(log(abs(fft2(block_struct.data)))));
B = blockproc(J, [64 10], fun);

然后输出是块补丁值的 [2000/64 x 500/10] 输出，您可以通过最近邻插值（或其他更平滑的版本）将其大小调整为所需的 [2000 x 500] 原始值大小

C = imresize(B, [2000 500], 'nearest');

如果有进一步的帮助，我可以提供一个真实的图像示例。

更新：要获得重叠块，您可以使用blockproc 的'Bordersize' 选项，方法是设置重叠[V H]，这样最终窗口大小[M + 2*V, N + 2*H] 仍将是[64 , 10] 大小。示例：

fun = @(block_struct) log(abs(fft2(block_struct.data)));
V = 16; H = 3; % overlap values
overlap = [V H]; 
M = 32; N = 4; % non-overlapping values
B1 = blockproc(J, [M N], fun, 'BorderSize', overlap); % final windows are 64 x 10

但是，这将适用于保留完整的傅立叶响应，而不是上面带有 max(max()) 的单值版本。

另请参阅此帖子以使用 blockproc:Dealing with “Really Big” Images: Block Processing 进行过滤。

Answer 2

如果您想对较大矩阵中的各个不同块应用相同的函数（在您的情况下为二维傅里叶变换），您可以使用 blkproc 函数来实现，该函数在较新的 MATLAB 版本中被替换blockproc.

但是，我推断您希望以“滑动窗口”方式将 apply fft2 应用于重叠块。为此，您可以将colfilt 与'sliding' 选项一起使用。请注意，我们在每个块上应用的函数是 fft：

block_size = [64, 10];
temp_size = 5 * block_size;
col_func = @(x)cellfun(@(y)max(max(abs(fft2(y)))), num2cell(x, 1), 'Un', 0);
B = colfilt(A, block_size, 10 * block_size, 'sliding', col_func);

这是如何工作的？ colfilt 通过将每个“滑动”块重新排列到新临时矩阵的单独列中来处理矩阵 A，然后将 col_func 应用于这个新矩阵。 col_func 依次将每一列恢复到原始块并对其应用fft2，返回每一列的最大振幅值。

需要注意的重要事项：

1. 由于这个提到的临时矩阵包括所有可能的“滑动”块，内存可能是一个限制。因此，为了在计算中使用更少的内存，colfilt 将原始矩阵A 分解为temp_size 的子矩阵，并分别对每个子矩阵进行计算。当然，得到的矩阵 B 还是一样的。

2. 结果矩阵B 中的每个元素都是根据相应的块邻域计算得出的。您的图像越大，您需要处理的块就越多，因此计算时间将呈几何级数增加。我相信您必须等待一段时间，直到 MATLAB 处理完您的 2000×500 矩阵上的所有个滑动窗口。