A Collection of Useful C++ Classes for Signal Processing

[antto]

thevinn: thanks, does it have the "matched z-transform" thing, like here:

Code: Select all

/* Digital filter designed by mkfilter/mkshape/gencode   A.J. Fisher
   Command line: /www/usr/fisher/helpers/mkfilter -Be -Lp -o 8 -z -a 1.2500000000e-01 0.0000000000e+00 -l */

#define NZEROS 0
#define NPOLES 8
#define GAIN   2.761491411e+00

static float xv[NZEROS+1], yv[NPOLES+1];

static void filterloop()
  { for (;;)
      { 
        xv[0] = next input value / GAIN;
        yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4]; yv[4] = yv[5]; yv[5] = yv[6]; yv[6] = yv[7]; yv[7] = yv[8]; 
        yv[8] =   xv[0]
                     + ( -0.0001374335 * yv[0]) + (  0.0022359877 * yv[1])
                     + ( -0.0169593084 * yv[2]) + (  0.0793664387 * yv[3])
                     + ( -0.2552563338 * yv[4]) + (  0.5931278473 * yv[5])
                     + ( -1.0129475617 * yv[6]) + (  1.2484472022 * yv[7]);
        next output value = yv[8];

i'm asking because this is just 8 poles now, alot easier on the CPU, yet the filter is still doing it's thing

[thevinn]

I don't include methods for matched z-transforms, and I probably won't. I can't compete with MZT because the coefficients are precalculated and incorporated as constants.

The other thing is that its a numerical method and not well suited for time varying of parameter changes. There are programs that will design those filters for you, no use reinventing the wheel.

[jphil123]

Yeah I should have posted an example. Anyway, I added the file "Documentation.cpp" to the sources and I added another template SimpleFilter. So here is the documentation:

Well, I've been running into issues with converting interleaved to non-interleaved format on iOS. Is there a possibility of passing in interleaved (LRLR) audio data? If not, any suggestions?

[thevinn]

Well, I've been running into issues with converting interleaved to non-interleaved format on iOS. Is there a possibility of passing in interleaved (LRLR) audio data? If not, any suggestions?

Deinterleave the data, process it, and then re-interleave it (if necessary) using the routines in DspFilters/Utilities.h:

Code: Select all

template <typename Td,
          typename Ts>
void deinterleave (int samples,
                   Td* left,
                   Td* right,
                   Ts const* src);

template <typename Td,
          typename Ts>
void interleave (int samples,
                 Td* dest,
                 Ts const* left,
                 Ts const* right)

[antto]

edit: solved

[Heureka]

Hi

Can someone please provide a simple recipe for using this filter design tool in a new VS project? I have DSP skills in Matlab but very little in C++ and I've struggled trying to do something very simple as described in the documentation.cpp file.

I made a CLR Console project and added the DspFilter source and header files. Then in the main method I put in the first example presented in documentation.cpp. I could build this project but when I inspected the filter pointer, f, I could only see the parameters and not the filter coefficients which is what I'm in need of.

So I suspected this had to do with the CLR Console project and I tried to use a Win32 Console project instead but then I ended up with namespace issues amongst others. I'm a newbie in terms of windows programming as you can hear.

Is there a friendly soul out there that can help out here?



Cheers

[thevinn]

1) are you using the latest set of sources?

2) can you post your sample code with an explanation of what you're trying to do?

Thanks

[Heureka]

1) are you using the latest set of sources?

2) can you post your sample code with an explanation of what you're trying to do?

Thanks

Hi,

..and thanks a lot for taking the time to reply!

Yes, I'm using the latest set of sources.

In my VS project I added the header and source files and then set the include path. Then I copy pasted an example from documentation.cpp to my main() as such:

Code: Select all

int main(array<System::String ^> ^args)
{
    // create a 2-channel Inverse Chebyshev Low Shelf of order 5 and passband ripple 0.1dB, without parameter smoothing and apply it.
    Dsp::Filter* f = new Dsp::FilterDesign <Dsp::ChebyshevII::Design::LowShelf <5>, 2>;
    Dsp::Params params;
    params[0] = 44100; // sample rate
    params[1] = 5; // order
    params[2] = 4000; // corner frequency
    params[3] = 6; // shelf gain
    params[4] = 0.1; // passband ripple
    f->setParams (params);

	char a = getchar();
    return 0;
}

For some reason in this CLR Console project environment intellisense doesn't work.

During debugging the filter pointer contains these data:


[0] 44100.000000000000 double
[1] 5.0000000000000000 double
[2] 4000.0000000000000 double
[3] 6.0000000000000000 double
[4] 0.10000000000000001 double
[5] 1.9898171591677428e+149 double
[6] 5.7102397540740374e-308 double
[7] 1.8154680098451641e-314 double

I got stuck here. [5:7] can't be coefficients so they are probably pointer addresses but I don't know how to dig deeper from here.

With the other Win32 Console project looking like this:

Code: Select all

#include "stdafx.h"
#include "stdio.h"

int _tmain(int argc, _TCHAR* argv[])
{
    {
    Dsp::Filter* f = new Dsp::SmoothedFilterDesign <Dsp::RBJ::Design::LowPass, 2> (1024);
    Dsp::Params params;
    params[0] = 44100; // sample rate
    params[1] = 4000; // cutoff frequency
    params[2] = 1.25; // Q
    f->setParams (params);
  }
 
	char c = getchar();
	return 0;
}

I got errors:

Error 5 error C2061: syntax error : identifier 'SmoothedFilterDesign'
Error 3 error C2065: 'f' : undeclared identifier
Error 12 error C2065: 'f' : undeclared identifier
Error 2 error C2065: 'Filter' : undeclared identifier
Error 9 error C2065: 'params' : undeclared identifier Error 10 error C2065: 'params' : undeclared identifier
Error 11 error C2065: 'params' : undeclared identifier
Error 14 error C2065: 'params' : undeclared identifier
Error 8 error C2065: 'Params' : undeclared identifier
Error 4 error C2653: 'Dsp' : is not a class or namespace name
Error 6 error C2653: 'Dsp' : is not a class or namespace name
Error 7 error C2653: 'Dsp' : is not a class or namespace name
....

I'm too green to know what the difference is between the two types of VS projects!

BTW. I'm in need of the coefficients only (not analysis or any realizations)

sorry for the noise.

Thanks

[thevinn]

Yes, I'm using the latest set of sources.

Yes you are, good.

Then I copy pasted an example from documentation.cpp to my main():

Code: Select all

    Dsp::Filter* f = new Dsp::FilterDesign <Dsp::ChebyshevII::Design::LowShelf <5>, 2>;

Ah you copied the wrong example. You used the Design API, by going through Dsp::Filter and operator new. This is a fully featured filter with polymorphism, introspection, and the ability to process channels, none of which you need. Unfortunately, the Dsp::Filter class interface is so abstract, that it is impossible to extract coefficients! That is only possible with one of the subclasses (like FilterDesign, SmoothedFilter, SimpleFilter, or by direct instantiation of a raw filter from its namespace).

For some reason in this CLR Console project environment intellisense doesn't work.

The DspFilters library makes heavy use of templates and unfortunately, things like intellisense simply will never be able to work properly.

During debugging the filter pointer contains these data:
[0] 44100.000000000000 double
[1] 5.0000000000000000 double
[2] 4000.0000000000000 double
[3] 6.0000000000000000 double

Right, you are looking at the filter parameters. [0] is sample rate, [1] is order, [2] is cutoff frequency and [3] is ripple. A Dsp::Filter pointer is just a pointer to an abstract virtual base class. All of the goodies are in the derived class and unfortunately probably not available for inspection in the debugger.

With the other Win32 Console project looking like this:

It looks like you simply forgot to #include "DspFilters/Dsp.h".

I'm too green to know what the difference is between the two types of VS projects!

Can't help you there, I only create Win32 applications and never anything else.

BTW. I'm in need of the coefficients only (not analysis or any realizations)

sorry for the noise.

Not a problem whatsoever! I am here to help you through it. Of course this library is rather complex and since it makes heavy use of templates, it is probably not possible to just look at the headers and figure out how to use it, thats why I provided examples. I added another example to Documentation.cpp (and checked it in), here it is:

Code: Select all

// Set up a filter, extract the coefficients and print
// them to standard output. Note that this filter is
// not capable of processing samples, as it has no
// state. It only has coefficients.
{
  Dsp::SimpleFilter <Dsp::RBJ::LowPass> f;
  f.setup (44100, // sample rate Hz
           440, // cutoff frequency Hz
           1); // "Q" (resonance)

  std::ostringstream os;

  os << "a0 = " << f.getA0 () << "\n"
     << "a1 = " << f.getA1 () << "\n"
     << "a2 = " << f.getA2 () << "\n"
     << "b0 = " << f.getB0 () << "\n"
     << "b1 = " << f.getB1 () << "\n"
     << "b2 = " << f.getB2 () << "\n";

  std::cout << os.str();
}

As you can see in this example we create a Dsp::SimpleFilter. This means no smoothing, no polymorphism, no virtual functions, no introspection, no "params", just a basic setup() function. The rest of the code extracts the coefficients of the biquad and outputs them to cout. If you wanted to, you could give the Dsp::SimpleFilter some channel state (as in the other example in Documentation.cpp).

Dsp::SimpleFilter should be low enough for you, but we can still go one step down:

Code: Select all

// Create an instance of a raw filter. This is as low as
// it gets, any lower and we will just have either a Biquad
// or a Cascade, and you'll be setting the coefficients manually.
{
  // This is basically like eating uncooked food
  Dsp::RBJ::LowPass f;
  f.setup (44100, 440, 1);

  // calculate response at frequency 440 Hz
  Dsp::complex_t response = f.response (440./44100);
}

Dsp::RBJ::LowPass is the raw filter that Dsp::SimpleFilter wraps up into a convenient package. As an example of usage, this calculates the complex valued response at a given frequency. Note that every function in Dsp::RBJ::LowPass would also be available when wrapped in a Dsp::SimpleFilter (you can still call response on it), since it uses simple derivation to inherits its methods.

So here's a summary of the objects that expose the various APIs:

Dsp::Filter*
Polymorphic, abstract interface to a filter instance that can process channels, optionally support parameter smoothing, provide introspection into parameter values, and also friendly GUI support (as used by DspDemo)

Dsp::FilterDesign<>
Subclass of Filter that wraps up a raw filter in a convenient interface that supports the Dsp::Filter API (polymorphism, parameter introspection, etc)

Dsp::SmoothedFilterDesign<>
A subclass of FilterDesign that supports smoothing of modulated parameters.

Dsp::SimpleFilter<>
A template suitable for wrapping up a raw filter, with the option to process channels, with no introspection, polymorphism; But also, no virtual function overhead whatsoever, and the possibility of having some calls inlined.

Here are examples of the raw filters:

Dsp::RBJ::LowPass
Dsp::RBJ::HighPass
Dsp::RBJ::BandPass
Dsp::Butterworth::LowShelf
Dsp::ChebyshevI::BandStop
etc...

After considering your use-case, I added a function to Cascade that will give you direct read-only access to the stage array so you can also extract coefficients from a higher order filter if desired.

If you still have questions, don't be shy...I want this thing to be used!

[Heureka]

Hi

You are a true elaborate gentleman

I can now compute the filter coefficients so this was great support.

I'm looking forward to work with your code.

Thanks a lot for effort.

Cheers

[thevinn]

DOCUMENTATION

All symbols are in the Dsp namespace.

class Filter

This is an abstract polymorphic interface that supports any filter. The
parameters to the filter are passed in the Params structure, which is
essentially an array of floating point numbers with a hard coded size
limit (maxParameters). Each filter makes use of the Params as it sees fit.

Filter::getKind ()
Filter::getName ()
Filter::getNumParams ()
Filter::getParamInfo ()

Through the use of these functions, the caller can determine the meaning
of each indexed filter parameter at run-time. The ParamInfo structure
contains methods that describe information about an individual parameter,
including convenience functions to map a filter parameter to a "control
value" in the range 0...1, suitable for presentation by a GUI element such
as a knob or scrollbar.

Filter::getDefaultParams ()
Filter::getParams ()
Filter::getParam ()
Filter::setParam ()
Filter::findParamId ()
Filter::setParamById ()
Filter::setParams ()
Filter::copyParamsFrom ()

These methods allow the caller to inspect the values of the parameters,
and set the filter parameters in various ways. When parameters are changed
they take effect on the filter immediately.

Filter::getPoleZeros ()
Filter::response ()

For analysis, these routines provide insight into the pole/zero arrangement
in the z-plane, and the complex valued response at a given normalized
frequency in the range (0..nyquist = 0.5]. From the complex number the
magnitude and phase can be calculated.

Filter::getNumChannels()
Filter::reset()
Filter::process()

These functions are for applying the filter to channels of data. If the
filter was not created with channel state (i.e. Channels==0 in the derived
class template) then they will throw an exception.

To create a Filter object, use operator new on a subclass template with
appropriate parameters based on the type of filter you want. Here are the
subclasses.



template <class DesignClass, int Channels = 0, class StateType = DirectFormII>
class FilterDesign : public Filter

This subclass of Filter takes a DesignClass (explained below) representing
a filter, an optional parameter for the number of channels of data to
process, and an optional customizable choice of which state realization
to use for processing samples. Channels may be zero, in which case the
object can only be used for analysis.

Because the DesignClass is a member and not inherited, it is in general
not possible to call members of the DesignClass directly. You must go
through the Filter interface.



template <class DesignClass, int Channels, class StateType = DirectFormII>
class SmoothedFilterDesign : public Filter

This subclass of FilterDesign implements a filter of the given DesignClass,
and also performs smoothing of parameters over time. Specifically, when
one or more filter parameters (such as cutoff frequency) are changed, the
class creates a transition over a given number of samples from the original
values to the new values. This process is invisible and seamless to the
caller, except that the constructor takes an additional parameter that
indicates the duration of transitions when parameters change.



template <class FilterClass, int Channels = 0, class StateType = DirectFormII>
class SimpleFilter : public FilterClass

This is a simple wrapper around a given raw FilterClass (explained below).
It uses inheritance so all of the members of the FilterClass are available
to instances of this object. The simple wrapper provides state information
for processing channels in the given form.

The wrapper does not support introspection, parameter smoothing, or the
Params style of applying filter settings. Instead, it uses the interface
of the given FilterClass, which is typically a function called setup()
that takes a list of arguments representing the parameters.

The use of this class bypasses the virtual function overhead of going
through a Filter object. It is not practical to change filter parameters
of a SimpleFilter, unless you are re-using the filter for a brand new
stream of data in which case reset() should be called immediately before
or after changing parameters, to clear the state and prevent audible
artifacts.



Filter family namespaces

Each family of filters is given its own namespace. Currently these namespaces
include:

RBJ: Filters from the RBJ Cookbook
Butterworth: Filters with Butterworth response
ChebyshevI: Filters using Chebyshev polynomials (ripple in the passband)
ChebyshevII: "Inverse Chebyshev" filters (ripple in the stopband)
Elliptic: Filters with ripple in both the passband and stopband
Bessel: Uses Bessel polynomials, theoretically with linear phase
Legendre: "Optimum-L" filters with steepest transition and monotonic passband.
Custom: Simple filters that allow poles and zeros to be specified directly

<class FilterClass>

Within each namespace we have a set of "raw filters" (each one is an example
of a FilterClass). For example, the raw filters in the Butterworth namespace are:

Butterworth::LowPass
Butterworth::HighPass
Butterworth::BandPass
Butterworth::BandStop
Butterworth::LowShelf
Butterworth::HighShelf
Butterworth::BandShelf

When a class template (such as SimpleFilter) requires a FilterClass, it is
expecting an identifier of a raw filter. For example, Legendre::LowPass. The
raw filters do not have any support for introspection or the Params style of
changing filter settings. All they offer is a setup() function for updating
the IIR coefficients to a given set of parameters.

<class DesignClass>

Each filter family namespace also has the nested namespace "Design". Inside
this namespace we have all of the raw filter names repeated, except that
these classes additional provide the Design interface, which adds
introspection, polymorphism, the Params style of changing filter settings,
and in general all of the features necessary to interoperate with the Filter
virtual base class and its derived classes. For example, the design filters
from the Butterworth namespace are:

Butterworth::Design::LowPass
Butterworth::Design::HighPass
Butterworth::Design::BandPass
Butterworth::Design::BandStop
Butterworth::Design::LowShelf
Butterworth::Design::HighShelf
Butterworth::Design::BandShelf

For any class template that expects a DesignClass, you must pass a suitable
object from the Design namespace of the desired filter family. For example,
ChebyshevI::Design::BandPass.

[Heureka]

Hi Again,

I can extract the RBJ filter coefficients but I can't seem to access the higher order filter coefficients like for example a 10th order Butterworth filter. This is what I'm doing:

Code: Select all

case LowPass:
{
 Dsp::Butterworth::LowPass <int(ButterworthMaxFilterOrder)> f_LP;
 f_LP.setup(int(args[0]),1.0,0.1);

  for(int i=0;i<numRows;i++)
  {
   array2D[i][0] = f_LP.m_stages[i].getB0();
   array2D[i][1] = f_LP.m_stages[i].getB1();
   array2D[i][2] = f_LP.m_stages[i].getB2();
   array2D[i][3] = f_LP.m_stages[i].getA0();
   array2D[i][4] = f_LP.m_stages[i].getA1();
   array2D[i][5] = f_LP.m_stages[i].getA2();
  }
  break;
...

..the arrays seems to be private. I could make them accessible by changing your code but if this is not intended I would leave this to you? I can't access "numStages" either.


I really fancy your JUCE application but I did notice that the group-delay changes oddly at some mid frequency. This even happens when the phase is smooth and well unwrapped so maybe there is an error?

Thanks again for sharing this.

[thevinn]

Are you working from the latest sources, or from the zip archive? The latest sources has public access to the stage array via operator [].

Regardless, I added getNumStages() to Cascade and also added another example in Documentation.cpp that does what you want, please do a SVN Update and you should have something that works. Let me know.

[thevinn]

Taking a look at your example code, it seems you overlooked the array operator overload. This is a working example (from the latest SVN trunk)

Code: Select all

  // Extract coefficients from a Cascade
{
  Dsp::SimpleFilter <Dsp::Butterworth::HighPass <3> > f;
  f.setup (3, 44100, 2000);

  std::ostringstream os;

  os << "numStages = " << f.getNumStages() << "\n"
     << "a0[0] = " << f[0].getA0 () << "\n"
     << "a1[0] = " << f[0].getA1 () << "\n"
     << "a2[0] = " << f[0].getA2 () << "\n"
     << "b0[0] = " << f[0].getB0 () << "\n"
     << "b1[0] = " << f[0].getB1 () << "\n"
     << "b2[0] = " << f[0].getB2 () << "\n"
     << "a0[1] = " << f[1].getA0 () << "\n"
     << "a1[1] = " << f[1].getA1 () << "\n"
     << "a2[1] = " << f[1].getA2 () << "\n"
     << "b0[1] = " << f[1].getB0 () << "\n"
     << "b1[1] = " << f[1].getB1 () << "\n"
     << "b2[1] = " << f[1].getB2 () << "\n"
    ;

  std::cout << os.str();
}

Note the use of the array index operator [] applied directly to f.

[Heureka]

Are you working from the latest sources, or from the zip archive? The latest sources has public access to the stage array via operator [].

Regardless, I added getNumStages() to Cascade and also added another example in Documentation.cpp that does what you want, please do a SVN Update and you should have something that works. Let me know.

Hi

Sorry for this delayed response.

I now have the code working. This is perfect.

Thank you again