Look at the filter
function.
If you just need a 1-pole low-pass filter, it's
xfilt = filter(a, [1 a-1], x);
where a = T/τ, T = the time between samples, and τ (tau) is the filter time constant.
Here's the corresponding high-pass filter:
xfilt = filter([1-a a-1],[1 a-1], x);
If you need to design a filter, and have a license for the Signal Processing Toolbox, there's a bunch of functions, look at fvtool and fdatool.
A couple of issue I can spot:
when you divide the signal by its mean, you need to check that it was not zero. Otherwise the result will be NaN
.
the authors (I am following this article) used a bank of filters with frequency bands covering the entire range up to the Nyquist frequency. You are doing half of that. The normalized frequencies you pass to butter
should go all the way up to 1
(corresponds to fs/2
)
When computing the energy of each filtered signal, I think you should not divide by its mean (you have already accounted for that before). Instead simply do: E = sum(sig.^2);
for each of the filtered signals
In the last post-processing step, you should normalize to the range [0,1], and then apply the median filtering algorithm medfilt2
. The computation doesn't look right, it should be something like:
img = ( img - min(img(:)) ) ./ ( max(img(:)) - min(img(:)) );
EDIT:
With the above points in mind, I tried to rewrite the code in a vectorized way. Since you didn't post sample input images, I can't test if the result is as expected... Plus I am not sure how to interpret the final images anyway :)
%# read biospeckle images
fnames = dir( fullfile('folder','myimages*.jpg') );
fnames = {fnames.name};
N = numel(fnames); %# number of images
Fs = 1; %# sampling frequency in Hz
sz = [209 278]; %# image sizes
T = zeros([sz N],'uint8'); %# store all images
for i=1:N
T(:,:,i) = imread( fullfile('folder',fnames{i}) );
end
%# timeseries corresponding to every pixel
T = reshape(T, [prod(sz) N])'; %# columns are the signals
T = double(T); %# work with double class
%# normalize signals before filtering (avoid division by zero)
mn = mean(T,1);
T = bsxfun(@rdivide, T, mn+(mn==0)); %# divide by temporal mean
%# bank of filters
numBanks = 10;
order = 5; % butterworth filter order
fCutoff = linspace(0, Fs/2, numBanks+1)'; % lower/upper cutoff freqs
W = [fCutoff(1:end-1) fCutoff(2:end)] ./ (Fs/2); % normalized frequency bands
W(1,1) = W(1,1) + 1e-5; % adjust first freq
W(end,end) = W(end,end) - 1e-5; % adjust last freq
%# filter signals using the bank of filters
Tf = cell(numBanks,1); %# filtered signals using each filter
for i=1:numBanks
[b,a] = butter(order, W(i,:)); %# bandpass filter
Tf{i} = filter(b,a,T); %# apply filter to all signals
end
clear T %# cleanup unnecessary stuff
%# compute average energy in each signal across frequency bands
Tf = cellfun(@(x)sum(x.^2,1), Tf, 'Uniform',false);
%# normalize each to [0,1], and build corresponding images
Tf = cellfun(@(x)reshape((x-min(x))./range(x),sz), Tf, 'Uniform',false);
%# show images
for i=1:numBanks
subplot(4,3,i), imshow(Tf{i})
title( sprintf('%g - %g Hz',W(i,:).*Fs/2) )
end
colormap(gray)
(I used the image from here for the above result)
EDIT#2
Made some changes and simplified the above code a bit. This shall reduce memory footprint. For example I used cell array instead of a single multidimensional matrix to store the result. That way we don't allocate one big block of contiguous memory. I also reused same variables instead of introducing new ones at each intermediate step...
Best Answer
There are several filters that can be used, and the actual choice of the filter will depend on what you're trying to achieve. Since you mentioned Butterworth, Chebyschev and Elliptical filters, I'm assuming you're looking for IIR filters in general.
Wikipedia is a good place to start reading up on the different filters and what they do. For example, Butterworth is maximally flat in the passband and the response rolls off in the stop band. In Chebyschev, you have a smooth response in either the passband (type 2) or the stop band (type 1) and larger, irregular ripples in the other and lastly, in Elliptical filters, there's ripples in both the bands. The following image is taken from wikipedia.
So in all three cases, you have to trade something for something else. In Butterworth, you get no ripples, but the frequency response roll off is slower. In the above figure, it takes from
0.4
to about0.55
to get to half power. In Chebyschev, you get steeper roll off, but you have to allow for irregular and larger ripples in one of the bands, and in Elliptical, you get near-instant cut off, but you have ripples in both bands.The choice of filter will depend entirely on your application. Are you trying to get a clean signal with little to no losses? Then you need something that gives you a smooth response in the passband (Butterworth/Cheby2). Are you trying to kill frequencies in the stopband, and you won't mind a minor loss in the response in the passband? Then you will need something that's smooth in the stop band (Cheby1). Do you need extremely sharp cut-off corners, i.e., anything a little beyond the passband is detrimental to your analysis? If so, you should use Elliptical filters.
The thing to remember about IIR filters is that they've got poles. Unlike FIR filters where you can increase the order of the filter with the only ramification being the filter delay, increasing the order of IIR filters will make the filter unstable. By unstable, I mean you will have poles that lie outside the unit circle. To see why this is so, you can read the wiki articles on IIR filters, especially the part on stability.
To further illustrate my point, consider the following band pass filter.
Now if you look at the zero-pole diagram using
zplane(b,a)
, you'll see that there are several poles (x
) lying outside the unit circle, which makes this approach unstable.and this is evident from the fact that the frequency response is all haywire. Use
freqz(b,a)
to get the followingTo get a more stable filter with your exact design requirements, you'll need to use second order filters using the
z-p-k
method instead ofb-a
, in MATLAB. Here's how for the same filter as above:Now if you look at the characteristics of this filter, you'll see that all the poles lie inside the unit circle (hence stable) and matches the design requirements
The approach is similar for
butter
andellip
, with equivalentbuttord
andellipord
. The MATLAB documentation also has good examples on designing filters. You can build upon these examples and mine to design a filter according to what you want.To use the filter on your data, you can either do
filter(b,a,data)
orfilter(Hd,data)
depending on what filter you eventually use. If you want zero phase distortion, usefiltfilt
. However, this does not acceptdfilt
objects. So to zero-phase filter withHd
, use thefiltfilthd
file available on the Mathworks file exchange siteEDIT
This is in response to @DarenW's comment. Smoothing and filtering are two different operations, and although they're similar in some regards (moving average is a low pass filter), you can't simply substitute one for the other unless it you can be sure that it won't be of concern in the specific application.
For example, implementing Daren's suggestion on a linear chirp signal from 0-25kHz, sampled at 100kHz, this the frequency spectrum after smoothing with a Gaussian filter
Sure, the drift close to 10Hz is almost nil. However, the operation has completely changed the nature of the frequency components in the original signal. This discrepancy comes about because they completely ignored the roll-off of the smoothing operation (see red line), and assumed that it would be flat zero. If that were true, then the subtraction would've worked. But alas, that is not the case, which is why an entire field on designing filters exists.