Step by step 1-D filtering
7 views (last 30 days)
Show older comments
Hello all.
I am in need of a 1-D lowpass filter in Matlab. My problem though is that I do not have all the measurements available. These are created in a step by step way and I need to filter after every step that generates a measurement. Furthermore I have no access to the signal processing toolbox and all its functions. I am aware that there is a possibility to use Zf and Zi in the filter function and I guess that can be used to do the step by step behaviour that I need. But I am unsure of how to achive a lowpass filter with that. Any help or directions would be most welcome.
0 Comments
Accepted Answer
Daniel Shub
on 29 Jun 2011
Step by step 1-D filtering is requires you to have a filter defined in terms of b and a. A very simple low pass filter is
b = [1, 1];
a = 1;
You will probably want to use something more complicated, but you haven't provided any details. Once the filter is designed, then you can do the following without the signal processing toolbox to filer step by step.
npts = 100;
nblocks = 5;
y = zeros(npts*nblocks, 1);
x = randn(npts, 1);
[y(1:npts), zi] = filter(b, a, x);
for iblock = 2:nblocks
x = randn(npts, 1);
[y((1:npts)+(npts*(iblock-1))), zi] = filter(b, a, x, zi);
end
0 Comments
More Answers (1)
Jan
on 11 Jul 2011
In this thread an M-file is posted, which emulates FILTER: Answers: use-filter-constants-to-hard-code-filter. This can be adjusted to solve your problem efficiently.
[EDITED]: Daniel asked, if this would be really efficient. Yes. If you work with a fixed filter, you can untroll the loops inside the filter function and get this for the b=[1,1], a=1 filter:
function y = StepFilter_Daniel(data)
npts = length(data);
nblocks = 5000;
y = zeros(npts*nblocks, 1);
b = [1, 1];
a = 1;
zi = 0;
for iblock = 1:nblocks
x = data;
[y((1:npts)+(npts*(iblock-1))), zi] = filter(b, a, x, zi);
end
return;
function y = StepFilter_Jan(data)
npts = length(data);
nblocks = 5000;
y = zeros(npts * nblocks, 1);
v = zeros(npts, 1);
z1 = 0;
for iblock = 1:nblocks
x = data;
for m = 1:npts
Xm = x(m);
v(m) = Xm + z1;
z1 = Xm;
end
y((1:npts)+(npts*(iblock-1))) = v;
end
return;
For a chunk length of 100 points, StepFilter_Jan is 45% faster than StepFilter_Daniel. The shorter the chunks, the larger is the advantage of the M-implementation without calling FILTER. For the original question for a step-by-step filtering the loop over m can be omitted an the M-version is 60% faster than calling filter.
4 Comments
Daniel Shub
on 27 Jul 2011
It always amazes me how inefficient some (possibly most) MATLAB functions are. I am blessed with having access to substantial computing power and computationally simple tasks, so I rarely care how long something takes.
Jan
on 27 Jul 2011
@Daniel: I've implemented a FILTER as MEX, which operates on the signal in backward direction. My implementation is naive (no multi-threading, no SSE, no assembler, but some unrolled loops) and equals the direct form II transposed structure as shown in "doc filter". I needed this backward processing for a faster FILTFILT to avoid reversing the signal twice. To my deep surprise my C-mex runs much faster than FILTER (at least of Matlab 2009a). Even after inserting the handling of arbitrary dimensions, checking inputs, working in both directions and accepting DOUBLEs and SINGLEs, the naive approach is still 2.5 times faster (see: http://www.mathworks.com/matlabcentral/fileexchange/32261-filterm ), but results are identical. Now I'm very curious how it is possible to write a *slower* function without inserting some SLEEP commands...
I've offered the code to TMW, but after waiting for some years without a reply, I decided to publish it in the FEX.
See Also
Categories
Find more on Multirate Signal Processing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!