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

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

现在，根据Matlab's filter documentation：

y = filter(b,a,x,zi,dim) acts along dimension dim. For example, if x is a matrix, then filter(b,a,x,zi,2) returns the filtered data for each row.

因此，要随着时间的推移平滑 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');

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

Bandpass Filter for 4D image in Matlab

matlab

filtering

signal-processing

lowpass-filter

highpass-filter