Bandpass Filter for 4D image in Matlab

Question

I have implemented in Matlab a bandpass filter for a 4D image (4D matrix). The first three dimensions are spatial dimensions, the last dimension is a temporal one. Here is the code:

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

I'm not sure if the implementation is correct. Could somebody verify it or does somebody find a failure?

Edit: Thanks to SleuthEye it now kind of works when I'm using bandpass_img = filter(fltr, 1, img, [], 4);. But there is still a minor problem. My images are of size 80x35x12x350, i.e. there are 350 time points. I have plotted the average time series before and after applying the bandpass filter.

Before bandpass filtering:

After bandpass filtering:

Why is this peak at the very beginning of the filtered image?

Edit 2: There is now a peak at the beginning and at the end. See:

I have made a second plot where I marked each point with a *. See:

So the first and last two time points seem to be lower.

It seems that I have to remove 2 time points at the beginning and also 2 time points at the end, so in total I will loose 4 time points.

What do you think?

Answer 1

Filtering a concatenation of all the elements using the 1-D filter function as you are doing is likely going to result in a distorted image as the smoothing makes the end of each row blend into the beginning of the next one. Unless you are trying to obtain a psychedelic rendition of your 4D data, this is unlikely to do what you are expecting it to.

Now, according to 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.

So, to smooth out the 3D images over time (which you indicated is the fourth dimension of you matrix img), you should use

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

Used as such, the signal would starts a 0 because that's the default initial state of the filter and the filter takes a few samples to ramp up. If you know the value of the initial state you can specify that with the zi argument (the 4th argument):

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

Otherwise, a typical order N linear FIR filter has a delay of N/2 samples. To get rid of the initial ramp up you can thus just discard those N/2 initial samples:

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

Similarly, if you want to see the output corresponding to the last N/2 input values, you'd have to pad your input with N/2 extra samples (zeros will do):

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