DSPFilters C++ library help

[makisce]

Hello all,

I have a vibration analysis software for iOS and Android that is used to collect vibration data from accelerometers and signal conditioners that use USB Audio as an interface. My software works fine on Android but on iOS i have some differences in the overall RMS values along with high ski slope on high-frequency vibrations due to bad filters that I was using. To cut the story short I did some time trying to resolve this issue on iOS and I found the only library that was not for "music studio" and on-screen piano software and this library is the DSPFilters C++ library.
I have been trying to implement the DSPFilters library (https://github.com/vinniefalco/DSPFilters) on iOS and in one of the files this forum was mentioned for further questions and discussion.

The library can't be used on modern XCode versions and Swift and therefore I have imported the code and used wrapper classes in order to make it work. I have done it successfully but I have one issue with the Butterowrth filters that are crashing without any human-readable error. I have used the exact piece of code as it was written in the examples but it does not work.

Can someone please help me resolve this issue?

Regards,

[matt42]

I dont know so much about that library, but it seems that the author doesn't post here anymore and I'm unsure if there is active support or updates.

What specifications are you looking for in your filter? There are probably much more portable solutions out there for you to use. You could check out http://jaggedplanet.com/iir/iir-explorer.asp for example.

[makisce]

Thank you for your answer. The author did not have any activity on GitHub either and the documentation itself is a disaster. The only example that I needed (Butterworth) crashes and so it says in the documentation, that that example produces an error (for demonstration I suppose).
I need filters that will calculate their coefficients automatically because when I make measurements I read the measurement configuration parameters from an application that is used to create routes which user then downloads to a device and goes on-site to take vibration measurements on machines. Those parameters are pre-configured per machine by a vibration specialist and after data is collected he performs off-site analysis. The data collection parameters have multiple filtering options and the user can select different cut-off frequencies, therefore filter input parameters should consist of only the sampling rate and the cut-off frequency and all coefficients should be automatically calculated. But the most important thing is that I need a high-pass filter that supports frequencies less than 10 Hz (5, 2, 1.5, 0.5). I have a beautiful library in Java (TarsosDSP) that has all that I need but for iOS is almost no existing because all filters that I managed to find had a much higher range.

[Music Engineer]

yes, i think, vinnie has moved on to other work. this here is the main thread in this forum about the library:
viewtopic.php?f=33&t=249926
...maybe you can find some help and useful information there

[mystran]

A continuous time (ie. "analog") Butterworth LP (and HP) of order N with normalised cutoff (we'll fix that below) can be decomposed into N/2 biquads (and a extra one-pole, if N is odd). All of these share the same cutoff (=1 for the normalised version), so we only have to compute Q factors that place the poles such that they form a half-circle (the comments below explain the math):

Code: Select all

// so we have a circle around 0
// at the normalized frequency
//  we have i^(pi/8) and i^(3*pi/8)

// then the pole pairs are simple..
//  (s - pp)*(s - conj(pp))
// so .. s^2 + 2*s*pp.re + pp.re^2 + pp.im^2
//  .. but we know cabs(pp) = 1, by definition
// so we can simply use sin here..
for(int i = 0; i < order / 2; ++i)
{
    q[i] = 1.0 / (2*sin((1+2*i)*pi/(2*order)));
}

Converting the LP prototype into a HP is a matter of placing 2 zeroes at DC (for Butterworth anyway) and this happens to be equivalent to converting each of the LP biquads into HP biquads, so we can actually handle it by simply taking a different output below [read: it's free in this particular case].

Now we have some unit-cutoff analog biquads with varying Q factors. For digital implementation, we could apply the bilinear transform directly, but as it turns out direct form filters have pretty poor numerical behaviour, especially at low cutoff values. Fortunately we can use trapezoidal integration for the classic analog state variable directly, which (after some optimisations) has roughly the same CPU cost (well, at least close enough that it's unlikely to matter). We just need to warp the cutoff:

Code: Select all

// coeffs, computed in setup
double f, r, g;
// state, initialise to zero
double z1 = 0, z2 = 0;

// fc = cutoff / sampleRate
void setup(double fc, double Q)
{
  f = tan(pi*fc);
  // convert Q to damping factor
  r = 1 / Q;
  // this is precomputed division for what is essentially
  // Gaussian elimination below.
  g = 1 / (f * (f+r) + 1);
}

double stepHP(double in)
{
  // compute the three "outputs" of the basic analog SVF
  double hp = (in - (f + r)*z1 - z2) * g;
  double bp = z1 + f*hp;
  double lp = z2 + f*bp;
  // integrate the new states
  z1 = 2*bp - z1;
  z2 = 2*lp - z2;
  // we wanted high-pass, so..
  return hp;
}

Now.. if you wrap this into a class and use one for each of the N/2 biquads with different Q factors (but the same cutoff for each) as computed above (filtering with each of these in series), you'll have Butterworth HP filters of arbitrary even orders (and odd orders if we add a one-pole, which I won't type here in this terrible editor until I know it's actually necessary).

ps. taking the reciprocal of the sine twice like this to turn a damping factor into Q and back into a damping factor is obviously pointless, but I just wanted to try to keep it more clear what is happening

pps. for the cutoff frequencies mentioned, at typical audio sampling rates, this might still work well enough with single-precision floats (since trapezoidal is a lot more tolerant than direct forms), but I would still suggest using doubles just to be safe

pps: note that taking the 2-pole "BP" output like this won't give you a higher-order Butterworth band-pass (but can still be used for the 2nd order case), since those need quite a bit more complicated math..

[mystran]

One more thing I want to note though: while Butterworth filters are "maximally flat" in spectral domain, they do oscillate (or "ring") somewhat in time-domain, especially for higher orders, and more importantly there is some non-linear phase-shift in the pass-band. I'm pointing this out, because if you're doing time-domain analysis of the results and want to remove low-frequency noise (eg. varying DC offsets), then using a more gradual design (eg. Bessel or some Gaussian approximation) can potentially preserve the waveform shapes better (even if you don't want to use linear-phase FIR filters), at the cost of more gradual transition around the cutoff. No idea if this is a concern for you, but just figured it might be worth pointing out.

[makisce]

One more thing I want to note though: while Butterworth filters are "maximally flat" in spectral domain, they do oscillate (or "ring") somewhat in time-domain, especially for higher orders, and more importantly there is some non-linear phase-shift in the pass-band. I'm pointing this out, because if you're doing time-domain analysis of the results and want to remove low-frequency noise (eg. varying DC offsets), then using a more gradual design (eg. Bessel or some Gaussian approximation) can potentially preserve the waveform shapes better (even if you don't want to use linear-phase FIR filters), at the cost of more gradual transition around the cutoff. No idea if this is a concern for you, but just figured it might be worth pointing out.

I do apologize for my late response.
The problem I'm facing is in finding appropriate filters for iOS, for audio signal that is generated by a piezoelectric accelerometer (see Digiducer 333D01 https://digiducer.com/). On Android, I use Tarsos DSP library which does the job perfectly for the purpose (results are mathematically validated and confirmed using appropriate vibration generator/shaker and confirmed with third-party hardware vibration analysis devices). We use the mobile app for collecting vibration parameters along with some other information which is then sent to our web application where vibration analysts can perform condition monitoring and troubleshooting.
For iOS I was unable to find proper filters because everything they have is for music and sometimes we have frequencies below 2Hz. When having really strong vibrations the higher frequencies can appear as a ski-slope after the signal is resampled. Therefore we filter the signal using Low-Pass filter with cutoff frequency equal to Fmax and a High-pass filter which can have values as low as 0.5 to 10 Hz, and then we do a resamplping and do the rest of the signal processing before calculating the overall RMS and the spectrum.
When I find out about this library I eliminated the ski-slope problem for Acceleration and Velocity by using HP and LP filters by using this code and the Q factor equal to 0.707:

For the Low-pass filter:

Code: Select all

float* audioData[1];
audioData[0] = audioSamples;
{
      Dsp::Filter* f = new Dsp::SmoothedFilterDesign
        <Dsp::RBJ::Design::LowPass, 1> (1024);
      Dsp::Params params;
      params[0] = samplingRate; // sample rate
      params[1] = cutOffFrtequency; // cutoff frequency
      params[2] = q; // Q
      f->setParams (params);
      f->process (numSamples, audioData);
}

For the High-pass filter:

Code: Select all

float* audioData[1];
audioData[0] = audioSamples;
{
    float* audioData[1];
    audioData[0] = audioSamples;
    {
        Dsp::SmoothedFilterDesign <Dsp::RBJ::Design::HighPass, 1> f (1024);
        Dsp::Params params;
        params[0] = samplingRate; // sample rate
        params[1] = cutOffFrtequency; // cutoff frequency
        params[2] = q; // Q
        f.setParams (params);
        f.process (numSamples, audioData);
    }
}

But I have issues with Demodulation (envelope). I have an audio file which I use for a reference and I run it on both Android and iOS and when I compare the results for Acceleration and Velocity I get almost identical values in both overall RMS and in top 10 peaks which are identical in both frequency and amplitude. For demodulation I have a big difference in iOS.
The demodulation algoirthm at one point has a Bandpass filter. As I'm not an audio specialist nor I know about filters, I can say that in the Android application I use Bandpass IIR filter whih takes the following parameters:
- Center frequency
- Bandwidth
- Sampling rate.

In order to calculate the center frequency (Fc) we take the lower cutoff (f1) and the upper cutoff frequency (f2).
** If f1/f2 >= 1.1 => Fc = sqrt(f1 * f2).
** If f1/f2 < 1.1 => fc = (f1 * f2)/2
The bandwidth is always f2 - f1.

And the filter is that simple. When I apply the same logic for this:

Code: Select all

    float* audioData[1];
    audioData[0] = audioSamples;
    {
      Dsp::Filter* f = new Dsp::SmoothedFilterDesign
        <Dsp::Butterworth::Design::BandPass <4>, 1, Dsp::DirectFormII> (1024);
      Dsp::Params params;
      params[0] = samplingRate; // sample rate
      params[1] = 4; // order
      params[2] = centerFrequency; // center frequency
      params[3] = bandWidth; // band width
      f->setParams (params);
      f->process (numSamples, audioData);
    }

I get results that do not have sense.

I forgot to mention that also the HP and LP filters are IIR filters and are taking only the cutoff frequency and the sample rate as parameters.

Is there anything that I can use to fix this issue?

[makisce]

And another thing. We use Hanning window and averaging in the signal processing. We save the normal waveform of the last average for further analysis.
And our web application is www.vib.cloud
Android: https://play.google.com/store/apps/deta ... bcloud_pro

[mystran]

Are you using the GPL version of Tarsos DSP that's on Github? If that's the case, then your application would also fall under the GPL license. If this is the case, you could just give me a link to the source code (since GPL requires this to be available anyway) and I could take a look at what the existing working code does.

[Ivan_C]

There is a filter design class in JUCE now btw