Step by step 1-D filtering

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.

Answer 1

[Daniel Shub]

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

Answer 2

[Jan]

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.

[Daniel Shub]

Have you timed your filter function versions the built-in function? It is not obvious to me where your solution is more efficient. I am assuming that somewhere in the builtin filter function is basically code like yours. Is the savings just the overhead of the MATLAB error checking?

[Jan]

@Daniel: Of course I measured the speed before I claim it would be efficient. See the EDITED code.

[Daniel Shub]

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]

@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.