File Exchange

image thumbnail

FilterM

version 1.0.0.0 (19.9 KB) by Jan
A faster FILTER and FILTFILT: Speedup factor 2.5 to 25

40 Downloads

Updated 20 Jul 2011

View License

FilterM, FiltFiltM: Fast digital filter

These functions are compatible to MATLAB's FILTER and FILTFILT commands,
but they are faster (see screenshot):
FilterM: 30%-40% of FILTER runtime
FiltFiltM: 4%-20% of FILTFILT runtime

ADDITIONAL FEATURES:
- The dimension to operate on can be specified for FiltFiltM.
- FilterM can process the signal in backward direction. (This is the
main part of the acceleration of FiltFiltM, because it avoids to
reverse the signal two times.)
- For signals of type SINGLE, the intermediate values are stored in
DOUBLE precision to increase the accuracy. The output is converted
to SINGLE again.
- The Signal Processing Toolbox is *not* needed.

CALLING:
Y = FiltFiltM(b, a, X, Dim)
[Y, Zf] = FilterM(b, a, X, Zi, Dim, Reverse)

b, a: Filter parameters as DOUBLE vectors.
X: Signal as DOUBLE or SINGLE vector or array.
Zi, Zf: Initial and final conditions as DOUBLE or SINGLE array.
Optional, default: Zeros.
Dim: Dimension to operate on. Optional, default: 1st non-singelton.
Reverse: Flag to process the signal in reverse direction.
Optional, default: FALSE.
Y: Filtered signal, same size and type as X.
While FilterM filters in forward direction, FiltFiltM processes
the signal forward and reverse direction for a zero phase
distortion.

INTENTION:
To accelerate my FEX submission FiltFiltM, I've implemented a filter as
C-Mex, which works in reverse order. To my surprise this was faster than
running Matlab's FILTER forward, e.g. 3.7 times for a [10000 x 1] vector,
5th order Butterworth filter (Matlab 2009a, WinXP 32 bit, single core).
Therefore I've expanded the Mex such that the direction can be defined
as input. The algorithm is a direct form II transposed structure.
A future version will be mutli-threaded.

INSTALLATION:
Setup the compiler if not done before: mex -setup.
Auto-compilation: Call FilterM without inputs to start the compilation.
A pre-compiled Mex can be downloaded: http://www.n-simon.de/mex
Run the unit-tests uTest_FilterM and uTest_FiltFiltM to check validity and speed.

Tested: Matlab 6.5, 7.7, 7.8, WinXP, 32bit
Compiler: LCC2.4/3.8, BCC5.5, OWC1.8, MSVC2008
Assumed Compatibility: higher Matlab versions, Mac, Linux, 64bit

This is faster and more powerful than my former submission "FiltFiltM", which will be removed soon.

Cite As

Jan (2020). FilterM (https://www.mathworks.com/matlabcentral/fileexchange/32261-filterm), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (41)

Johannes Rebling

Jan

@Jim: I stopped at unrolling 6 loops, because it was a tedious job. I have a multi-threaded version also, but I'm still not happy with the distribution to the single threads. If you want to filter a 4000x4 matrix along the 1st dimension, 4 threads are a perfect choice. But for filtering along the 2nd dimension, it would be useful to process the matrix in blocks. But how to determine a useful block size?

Jim Hokanson

@Feng Cheng, I realize this is a bit of a delayed response, but for others ... It is only slower when multiple columns are involved. My guess is that Matlab splits the columns up across multiple processing threads. This code does not have multi-threading enabled.

Jim Hokanson

Hi Jan, I noticed you unrolled up to 6 coefficients in the version I have (from a while ago). Any chance you remember if it mattered that you stopped there? I have a filter with 8 coefficients ...

Feng Cheng

OS: MAC
MATLAB: R2018b
## It looks like R2018b has updated built-in filter function because FilterM looks much slower.
```
Known answer tests:
ok: all input empty
ok: empty parameters
ok: empty signal
ok: [N x 1] vector
ok: [1 x N] vector
ok: 2D array, 1st dim
ok: 2D array, 2nd dim
ok: 3D array, 1st, 2nd, 3rd dim
ok: [1 x 1 x N] array, 1st non-singelton
ok: reverse
ok: SINGLE array
Difference to DOUBLE divided by EPS(SINGLE(1)):
FilterM: 0.249992, FILTER: 5.36918
ok: SINGLE array, reverse
ok: permutated array

Provoke problems:
ok: Too short signal in dim 2 caught.
ok: Dim out of range caught.
ok: Not integer Dim caught.

== Speed test with Butterworth filter of order 4:
Loops determined for about 2 seconds run time
1000000 x 1: 238 loops
filter: 1.951 sec
FilterM: 0.811 sec ==> 41.6%
Results are identical

100000 x 1: 1942 loops
filter: 1.662 sec
FilterM: 0.691 sec ==> 41.6%
Results are identical

10000 x 10: 1242 loops
filter: 0.413 sec
FilterM: 0.473 sec ==> 114.5%
Results are identical

1000 x 100: 1496 loops
filter: 0.336 sec
FilterM: 0.533 sec ==> 158.3%
Results are identical

100 x 100: 6206 loops
filter: 0.286 sec
FilterM: 0.358 sec ==> 125.2%
Results are identical

100 x 1000: 1432 loops
filter: 0.331 sec
FilterM: 0.522 sec ==> 157.5%
Results are identical

== FilterM passed the tests.
```

Jan

The algorithm for Matlab's FILTFILT has been changed in modern Matlab versions. The unit-test function accepted differences up to 10 eps, but 20 eps is required. Change the line 268 in uTest_FiltFilt.m to:
elseif all(abs(origy(:) - newy(:)) < 20 * eps) % 10 => 20
If this still fails, please insert this line after line 268:
fprintf(' Max difference: %g\n', max(abs(origy(:) - newy(:)) / abs(newy)));

Does the test run now? I'm posting an updated version soon.

Maria Geffen

I am getting a "Different results - DO NOT USE FiltFiltM!" error when I run the uTest_FiltFiltM function. Any suggestions as to why this is the case?

Compiled using VS 2013, run in Matlab 2018b on Win 10.

Vedran Grudenic

Hi this is really nice, there is one tiny edge case where the function doesn't work correctly and writes out of bounds: when both `b_in` and `a_in` only have 1 element. In this case, `order` becomes zero, and then the function writes into `Z[order - 1]` inside `CoreDoubleN`. So there should perhaps be a check for this edge case somewhere at the beginning.

Farhad Sedaghati

Siddhi Imming

When compiling with 'mex FilterX.c -compatibleArrayDims' as command it works fine on windows 10 x64 using MATLAB 2018a.

Lukas Leuenberger

@Michele and @Silvia
The file needs to be compiled with the switch "-compatibleArrayDims". Otherwise some conditionals won't work correctly, because the type "mwSize" is treated as unsigned instead of signed.

Michele Giugliano

on R2017a, the test suite crashes MATLAB every time it runs. :-(

Silvia Casarotto

Dear Jan,
I have succesfully used FiltFiltM.m to speed up signal filtering.
Now I would like to use the compiled version. I installed the MinGW compiler Add-on and then I run FilterM(). The compiling process terminates succesfully and a mex file is generated.
However, when I run either uTest_FilterM or uTest_FiltFiltM, Matlab crashes and automatically closes. The error message that pops-up says something like: "instruction 0x122654782 has referred to memory 0xea725374. Memory could not be written". This also happens when I try to filer random signals. I am working with Matlab 2017a 64bit and signals are stored in the workspace as double. Could you please help me with this issue?

Thank you very much,

Silvia

Jan

@Danijel: The inputs are the same as for CoreSingleN, where you find the documentation.
MX: Number of elements in the 1st dimension
NX: Number of columns, considers mutliple dimensions.

Danijel

Hi, when using CoreSingle5() in my C project, what are the parameters MX and NX? Can't figure out. My input is 64 float samples.

dhaba india

Kevin C

Hi Jan,
Do you still plan to implement multithreading in an upcoming version? If so, do you have any idea when that version will be released? It would be a big help for me.

Kevin C

Andrea Libri

Royi Avital

Hi,
Used "FilterM" in my project:
https://github.com/RoyiAvital/FastGuassianBlur

The project is about finding the most efficient method to apply a Gaussian Blur on an image.

It really made things much much easier.

Thank You.

Jan

I've received the files. Thanks! I see that I've made a mistake with the matlab.2015 address - sorry!

Michele S.

Hi Jan,
sorry to bother you, but I had to struggle quite a bit to find your e-mail address :S
Eventually I sent some files at: matlab.2010@...
Could you confirm that you received them? Thanks!

Jan

@Michele S.: The currently published version does not consider the special cases a==1 or b==1. I have a new version with a special treatment of these cases, but still struggle with the multi-threading. Please send me the data by email for some experiments (you find the address in the help section of the code). Thanks.

Michele S.

Hello Jan,
I found out that in my actual use case, FilterX performs significantly worse than the built-in filter (~ x4 slower). In my case a = 1, b is a 1445x1 double array, X is a 97280x1 double array.
I am running on r2011b, 32 bit, and the mex file is compiled with MSVC2008.
Is this to be expected due to the specific signal/coefficients sizes? Any comments or suggestions would be higly appreciated.
I can share some data and examples with you if you are willing to look further into this.
Thanks

Royi Avital

Jan,
Why not add an optional parameter of the number of threads to use?
Then you move the "Learning and Optimization" process to be the responsibility of the user.

We'd love to have a multi threaded version of this.

Thank You.

Jan

@Jim: I'm sorry that my FEX submissions makes you nervous. I attach the unit-tests functions with speed comparisons to reduce this kind of stress.
I'm gald that this function helps you to save time and energy.

@Yair: I use (sulking) camel-case (this means: some uppercase letters) on purpose to avoid name space collisions with Matlab's builtin functions.
I have some succes with the multi-threaded version, but I still struggle with a general way to determine the most efficient number of threads. If the memory transfer is the bottle neck, even two threads on a single core can be useful. For data with short columns but many rows the strategy differs substantially for machines with few or many cores. The best would be an automatic learning process during the installation, which adjust the parameters for the number of threads depending on the data size and filter length. But this is an overkill and a source of bugs.
It would be useful to have a wide-spread feedback from FEX users with real-world cases, but this is seldom: only 1 of 20'000 downloaders starts a discussion with me.

Jim Hokanson

Hi Jan,

I'll be honest I'm always a bit nervous about your performance specs given updated Matlab versions and my use cases but in my very limited testing (single channel, 1 GB data, filtfiltM only) it seems to do much better. I'm getting a runtime that is 40% of the original and even better, the memory usage is roughly around 40% of the original. Given the large memory requirements of the original this feature is quite helpful.

Thanks,
Jim

Yerlan Aitzhan

Yair Altman

Still 20-30% faster than Matlab's built-in, even after the performance boost to the built-in filter function in R2014a.

Nice documentation, unit-tests and optimization for good performance.

Suggestion: rename the functions filterm, filterx etc. (without capitals). The current FilterM file breaks because it refers to filterX instead of FilterX for example. Now that Matlab is case-sensitive, better make all functions lowercase, as is the Matlab convention.

Aside from this minor issue, a very useful utility. Well done!

Any progress with the promised multi-threaded version?

Rami Khushaba

Carlo E.D. Riboldi

Thanks for these great functions. Besides an improvement in speed, they apparently solve some instability of default filtfilt function causing NaN in the output on a randomic basis.

Royi Avital

@Jan, try even longer filters. I'd be happy if you could release an update for long FIR's even before the threaded version.
It would be great!

Thank you for this great function.

Jan

@Royi Avital: The currently published version does not consider the reduced calculations for FIR filters. I get these measurements for 2011b/64 under Windows, MSVC 2008:
x=rand(1e5, 1); B=repmat(1/21, 1, 21);
tic; for i=1:100, y=filter(B,1,x); clear('y'); end, toc
>> 0.56 sec
tic; for i=1:100, y=FilterX(B,1,x); clear('y'); end, toc
>> 0.42 sec (still faster)
The new version of FilterX considers A=1 explicitely and requires 0.27 sec only. It will be included in the next version, which is multi-threaded also. See: http://www.mathworks.com/matlabcentral/answers/24384-multi-threaded-mex-functions-in-the-fex

Royi Avital

@Jan, for FIR Filters with more than 20 coefficients FilterM lags behind MATLAB's implementation.
Do you think you could optimize for longer filters as well?
Thank you.

Jan

Sorry, currently only real values are considered. The implementation of a complex parameters is not trivial in C. It is planned later.

Royi Avital

Jan, there's an issue with complex coefficients, Could you check that?

I'm using the MEX file directly (For "Classic" use as Filter).

Thank you.

Jan

@Royi: I'm working on a multi-threaded version at the moment. Unfortunately this will work under Windows only, because I cannot test the pthread libs under Linux and MacOS.

Royi Avital

Thank you, this is great!
Do you think further optimization could be made?
Do you intentions to do so?

Jan

@Vassilis: I do not get the point. FilterM produces exactly the same results as Matlab's FILTER (except for the higher accuracy for SINGLE), therefore it leads to the same phase shift. FiltFiltM (from this submission) produces the same results as Matlab's FILTFILT. The former submission 25337 "Patch FiltFilt" is slower, depends on the SPT and does not work with modern Matlab versions - but it creates the same results as FiltFiltM.

If you get different results between FilterM and FILTER or between FiltFiltM and FILTFILT, please send me the output of unit-test functions by mail.

vkehayas

Very good addition Jan. I observed however that, as you might know, it produces different results than your FiltFiltM as well as from the built-in filtfilt function. I think this is because of the zero-padding you apply, while in the previous version you were using Gustafsson's method -- am I right? The most prominent difference is a time-lag to the signal, plus some minor distortions in the signal.

vkehayas

MATLAB Release Compatibility
Created with R2009a
Compatible with any release
Platform Compatibility
Windows macOS Linux
Acknowledgements

Inspired: filter1, Fast Guassian Blur

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!