Matlab中4D图像的带通滤波器

Question

我在 Matlab 中为 4D 图像（4D 矩阵）实现了带通滤波器。前三个维度是空间维度，最后一个维度是时间维度。代码如下：

function bandpass_img = bandpass_filter(img)
% Does bandpass filtering on input image
%
% Input:
%   img: 4D image
%
% Output:
%   bandpass_img: Bandpass filtered image

TR = 1; % Repetition time
n_vols = size(img,3);
X = [];

% Create matrix (voxels x time points)
for z = 1:size(img,3)
    for y = 1:size(img,2)
        for x = 1:size(img,1)
            X = [X; squeeze(img(x,y,z,:))']; %#ok<AGROW>
        end
    end
end

Fs = 1/TR;
nyquist = 0.5*Fs;

% Pass bands
F = [0.01/nyquist, 0.1/nyquist];
type = 'bandpass';

% Filter order
n = floor(n_vols/3.5);

% Ensure filter order is odd for bandpass
if (mod(n,2) ~= 0), n=n+1; end
fltr = fir1(n, F, type);

% Looking at frequency response
% freqz(fltr)

% Store plot to file
% set(gcf, 'Color', 'w');
% export_fig('freq_response', '-png', '-r100');

% Apply to image
X = filter(fltr, 1, X);

% Reconstructing image
i = 1;
bandpass_img = zeros(size(img));
for z = 1:size(img,3)
    for y = 1:size(img,2)
        for x = 1:size(img,1)
            bandpass_img(x,y,z,:) = X(i,:)';
            i = i + 1;
        end
    end
end

end

我不确定实施是否正确。有人可以验证它还是有人发现失败？

编辑：感谢 SleuthEye，当我使用 bandpass_img = filter(fltr, 1, img, [], 4); 时，它现在可以正常工作了。但是还有一个小问题。我的图像尺寸为 80x35x12x350，即有 350 个时间点。我已经绘制了应用带通滤波器前后的平均时间序列。

带通滤波前：

带通滤波后：

为什么这个峰值出现在过滤图像的最开始？

编辑 2：现在在开始和结束处都有一个峰值。见：

我制作了第二个图，其中我用 * 标记了每个点。见：

所以第一个和最后两个时间点似乎更低。

看来我要去掉开头的2个时间点，最后还要去掉2个时间点，所以一共要去掉4个时间点。

你怎么看？

Answer 1

使用一维filter 函数过滤所有元素的串联可能会导致图像失真，因为平滑会使每一行的末尾混合到下一行的开头。除非您试图获得 4D 数据的迷幻再现，否则这不太可能达到您的预期。

现在，根据Matlab's filter documentation：

y = filter(b,a,x,zi,dim) 作用于维度dim。例如，如果x 是一个矩阵，则filter(b,a,x,zi,2) 返回每​​一行的过滤数据。

因此，要随着时间的推移平滑 3D 图像（您指出这是矩阵 img 的第四维），您应该使用

bandpass_img = filter(fltr, 1, img, [], 4);

这样使用时，信号将从 0 开始，因为这是滤波器的默认初始状态，并且滤波器需要一些样本来加速。如果您知道初始状态的值，您可以使用 zi 参数（第 4 个参数）指定它：

bandpass_img = filter(fltr, 1, img, zi, 4);

否则，典型阶 N 线性 FIR 滤波器的延迟为 N/2 个样本。为了摆脱最初的加速，您可以丢弃那些 N/2 初始样本：

bandpass_img = bandpass_img(:,:,:,(N/2+1):end);

同样，如果您想查看与最后一个 N/2 输入值相对应的输出，则必须使用 N/2 额外样本填充您的输入（零即可）：

img = padarray(img, [0 0 0 N/2], 0, 'post');