markert / fili.js

Demo Implementation
http://markert.github.io/fili.js/
MIT License
94 stars 25 forks source link

Matlab to Fili equivalent #11

Closed rclai closed 8 years ago

rclai commented 8 years ago

Hello, I'm trying to convert this little MATLAB snippet in hopes of being able to produce the same result on your library:

dataSet = []; % Some data set
filtCutOff = 0.001;
frequency = 256;
% butterworth coefficients `a` and `b`
[b, a] = butter(1, (2*filtCutOff)/frequency, 'high');
filteredDataSet = filtfilt(b, a, dataSet);

butter documentation filtfilt documentation

My Fili code:

const SAMPLE_FREQUENCY = 256;
const HPF_CUTOFF = 0.001;
const iirCalculator = new Fili.CalcCascades();

const hpfCoeffs = iirCalculator.highpass({
  order: 1,
  characteristic: 'butterworth',
  Fs: 256,
  Fc: HPF_CUTOFF * 2 / SAMPLE_FREQUENCY,
  gain: 0,
  preGain: false,
});

const highpassFilter = new Fili.IirFilter(hpfCoeffs);
highpassFilter.simulate([/* data set */])

I am running this filter against a dataset that looks like this:

0.968589936319484
0.961296120425768
0.960164493014417
0.983637865121265
0.991536969658882
0.983507938439623
0.960591491533508
0.968489501837304
0.968526919825609
0.976329483746718

And after filtering it in MATLAB, I'm supposed to get this:

0.0228847850615585
0.0155911615511576
0.0144597265334421
0.0379332910463414
0.0458325880105157
0.0378037492430232
0.0148874948092062
0.0227856975989088
0.0228233080856626
0.0306260645201376

But after filtering it with my Fili code, I am getting this:

0.9685898049922316
0.9612957274329681
0.9601638394984774
0.9836369480529923
0.9915357847846424
0.9835064857770703
0.9605897752784536
0.9684875240260702
0.9685246793823019
0.9763269796084018

Would you happen to know where I went wrong?

markert commented 8 years ago

Hi, after a little digging into matlab manuals, I guess, I found some differences. What fili designs is a one biquad butterworth filter. A biquad is always a two stage filter. fili works with cascades of two stage filters. Another difference is that the cutoff frequencies do not match. fili gets the sampling and the cutoff frequency in Hz. So, your filter expects a 256 Hz signal and implements a high pass at 0,000007Hz. This will never do something that makes any sense. Getting close to 0 or fs/2 will make filters unstable, especially if the filter order is low. The matlab implementation on the other hand produces a 1 stage butterworth filter with a cutoff frequency of 0,000007% of the nyquist frequency. For a 256Hz signal this would be 0,0000039Hz. This filter will also behave quite strange.

rclai commented 8 years ago

Thanks for taking the time to look at this.

If you're up for further discussion, I believe it is intentional for the filter to make values close to zero. To get a better idea of what I'm doing, I am filtering out the influence of gravity on accelerometer readings. I'm essentially trying to replicate this section of code that simply tries to detect stationary periods from acceleration data. Does this make sense?

markert commented 8 years ago

No problem.

I see why you are setting the cutoff frequency close to 0. When calculating in double precision this may work. However, it will lead to results that differ by many ppm from the actual value due to rounding errors. I am responsible for the development of high precision measurement devices (weight, force, torque etc.) where this would not be acceptable. Implementing these filters as fixed point or single precision filters in an FPGA or microcontroller will also not work.

Your code snipped devides the signal spectrum into a low pass and a high pass part. Now it depends on what your latency requirements and processing power how you can implement that.

  1. IIR highpass and lowpass - low latency, potentially unstable, low complexity
  2. FIR highpass and lowpass (e.g. Kaiser-Bessel) - high latency, always stable, high calculation complexity, stationary periods will be detected
  3. FFT - Do an FFT with the signal -> set all parts you are not interested in to 0 -> do inverse FFT - high latency, high calculation complexity, always stable, frequencies in signal should not change fast
  4. wavelet transform - this really depends on what you would like to see in the signal, detecting stationary periods can be done, too