A Collection of Useful C++ Classes for Signal Processing

[toncho11]

Hi,

Your framework looks very promising!

I would like to create a Butterworth BandPass filter between 1 and 31 Hz. I am working with EEG signal form the brain usually in the range of 1Hz - 100Hz. In this case I would like to only use the signal between 1 and 31Hz.

How can I do that?

I do:

Dsp::Filter* f = new Dsp::SmoothedFilterDesign
<Dsp::Butterworth::Design::BandPass <4>, 2, Dsp::DirectFormII> (1024);

So the Center frequency should be (1+31) / 2 = 16.
And the with frequency (or the bandwidth) should be the difference between the two frequencies? So this is: 16 - 1 = 15

So in my case:

Dsp::Params params;
params[0] = 128; // sample rate
params[1] = 4; // order
params[2] = 16; // center frequency
params[3] = 15; // band width

will filter the signal letting only the frequencies between 1 and 31 with sampling rate 128Hz.

Is it correct?

Thanks in advance,

[JCJR]

I don't know anything about these C++ classes, but a typical symmetrical bandpass filter, if calculating the center frequency from the upper and lower cutoff frequencies, you would use the geometric mean.

Fc = SQRT(Fl * Fu)

So for Fl = 1 and Fu = 31, Fc = SQRT(1 * 31) = 5.568 Hz

The reason for that, is that a typical symmetrical bandpass filter is symmetrical by octaves rather than by absolute frequency. For instance, the attenuation would be identical at Fc / 1.5 and also Fc * 1.5. And the attenuation would be identical at Fc / 4 and also Fc * 4. And so forth.

With your center frequency of 5.568, the attenuation would be identical at 5.568 / 5.568 = 1, and also 5.568 * 5.568 = 31.

[toncho11]

Hi,

I would like to use bandpass filter with "A Collection of Useful C++ Classes for Digital Signal Processing" in real-time. This means that I need to be able to apply the filter sample per sample or at least per chunk. What is my best option?

I hear that IIR filters are actually better for real-time signal processing.

Can I just use Butterworth bandpass filter and apply it per sample/chunks? Some overlapping is probably required?

I will appreciate some feedabck here

I am new to DSP.

[aciddose]

Even if you can only process a fixed block at a time, it is possible to ensure you always have this size of block regardless of the chunk currently being processed by using a fixed latency (a delay) equal to the desired block-size.

That can be reported to the host which will then (assuming it can) compensate the latency for you.

If that is unacceptable then yes you need to process individual samples.

[Theo M]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

Hi,

First of all, these classes are great! Very easy to understand.

But (there is always a but...) i've got a problem and I don't know if it comes from these filters, or my method is just wrong.

What I would like, is a clean butterworth filter on the third octave center on 100Hz. But I'm working with a 48000Hz sample rate.
If I create my butterworth with this sample rate, I don't have enough point to create a clean filter.

Example :

Code: Select all (#)

/Get the buffer
memcpy(inputBuffer[0], buffer.getSampleData(0,0), size_Buffer*sizeof(float));	

//Define the band pass filter : 
Dsp::Filter* _filter = new Dsp::SmoothedFilterDesign <Dsp::Butterworth::Design::BandPass <6>, 1, Dsp::DirectFormII> (0);

//Parameters
Dsp::Params _para;
_para[0] = freq_echantillonnage; // sample rate (48000Hz in this case)
_para[1] = 6; // order
_para[2] = 100; // center frequency
_para[3] = 23; // band width   

//Process
_filter->setParams (_para);
_filter->process(size_Buffer, inputBuffer);

buffer.copyFrom(0, 0, inputBuffer[0], size_Buffer);

Here is what I got :

UglyButter.PNG

Ok, so I presume it's because of the sample rate.
To resolve this problem, I wanted to apply a downsampling on my input buffer to get a 12000Hz sample rate (the filter looks better like this).

So :

Code: Select all (#)

//Get the buffer
memcpy(inputBuffer[0], buffer.getSampleData(0,0), size_Buffer*sizeof(float));	

//Design the band pass filter : 
Dsp::Filter* _filter = new Dsp::SmoothedFilterDesign <Dsp::Butterworth::Design::BandPass <6>, 1, Dsp::DirectFormII> (0);

//Parameters
Dsp::Params _para;
_para[0] = freq_echantillonnage; // sample rate (48000Hz in this case)
_para[1] = 6; // order
_para[2] = 100; // center frequency
_para[3] = 23; // band width   

//LowPass
Dsp::SimpleFilter<Dsp::Butterworth::LowPass <18>, 1 > k;

k.setup(18, freq_echantillonnage, freq_echantillonnage/8);	//Lowpass at 6000Hz
k.process(size_Buffer, inputBuffer);

//Some temporary buffer...
float *copyIn[1];
copyIn[0] = new float [size_Buffer/4];

//I get one sample on 4 -> sample rate : 12000Hz
for (int j=0; j < size_Buffer/4; j++)
{		
       copyIn[0][j] = inputBuffer[0][4*j];
}

_filter->setParams (_para);   
_filter->process(size_Buffer/4, copyIn);

buffer.copyFrom(0, 0, copyIn[0], size_Buffer/4);

This is what I got :

BetterButter.PNG

So, the fact is that this is not my real code and in reality, my butterworth is cleaner, but, what I don't understand, this is the -12dB at the top of the filter (don't know the word in english). Why would I lose 12dB? This is way too much than I can afford...
I really don't understand what is my problem here....

Thank you very much if someone can help me with that... I'm stuck with this for way too long.

Thanks in advance.

(sorry for bad english... I hope it's understandable)

edit : Another question : Is-it possible to process a filter with any number of samples? I mean, is this possible to process 48000 sample in a single filter?
edit 2 : If it helps, I use "VST Plugin Analyser" to visualize the filter response.

[Astro11]

I have recently downloaded dspfilterscpp. Great job, it easily added to my project and is quite simple to use.
I am having a problem using Bessel filters. Looks like the cutoff frequency is not as specified, and it seems to changed when filter order is changed. Butterworth and Chebyshev filters seem correct. Any ideas?

[camsr]

Bessel filters have different gain values at the same corner.

[JCJR]

I have recently downloaded dspfilterscpp. Great job, it easily added to my project and is quite simple to use.
I am having a problem using Bessel filters. Looks like the cutoff frequency is not as specified, and it seems to changed when filter order is changed. Butterworth and Chebyshev filters seem correct. Any ideas?

I know nothing of that library, but a possibly useless comment--

A "simple specified" butterworth filter, Q = 0.7071, for instance one built with analog components-- The -3 dB frequency will be the same as the Fc of the filter.

Bessel (Q < 0.7071) and various Chebyschev's (Q > 0.7071)-- If you just set the required Q, but calculate the same Fc as with a butterworth-- It will have the expected rolloff shape, but the -3 dB point WILL NOT be at Fc.

With Q < 0.7071, the gain at Fc will be somewhat higher than -3 dB, and if you want the gain to be exactly -3 dB at a desired frequency then you have to use a lower Fc than your desired -3 dB point.

With Q > 0.7071, the gain at Fc will be less than -3 dB, and if you want the gain to be exactly -3 dB at a desired frequency then you have to use a higher Fc than your desired -3 dB point.

Some calculational methods take this into account when specifying the -3 dB frequency, so that the actual Fc is automatically scaled to make it all work out as desired. Assuming that a person definitely needs to specify an exact -3 dB point.

But simply setting Fc and Q without some kind of adjustment, will only guarantee that the -3 dB frequency agrees with the Fc, in the case of the Butterworth filter shape.

Also, at least in analog, > 2nd order filters, made by cascaded second-order and first-order sections-- To get an accurate "shape" to conform to the ideal of bessel or chebyschev, the Fc and Q of each cascaded section will usually be different for each section, and the Fc of none of the sections will be expected to be at the -3 dB point of the composite filter.

I'm just saying, unless somebody coded those smarts into a module, then perhaps one will get the desired filter shape, as long as he doesn't assume that the -3 dB point will be the same as the Fc.

I'm very ignorant of the topic, merely parroting stuff from old analog books I read long ago.

[maxthedog]

Hi,

I'm trying to make an stl vector of pointers to Dsp::Filter, but this causes xcode to freak out when trying to compile ( I don't have my laptop handy, but IIRC, it produces a lot of complaints that look reminiscent of linker errors).

For instance:

std::vector<Dsp::Filter*> filterbank;

Essentially, I'm looking for a way of creating a filter bank of bandpass filters on the fly, and then later running them as desired and fetching their results.

Any advice or alternative solutions would be greatly appreciated.

Thank you in advance for your time.

[maxthedog]

^Hmmmm... I believe I found the problem! It's me.

[olalp]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

Hi,

I downloaded this library and easily integrated it with my code but I am running into some problems and I am not sure how to debug it. Here is how I am instantiating the filter:

Code: Select all (#)

 
void DSPfilter::createFilter(dspfilt * f)
{
    if(m_bFilter)
    {
        delete filter;
        m_bFilter = false;
    }

    switch(f->type)
    {
    case HIGHPASS_BUTTERWORTH:
        filter = new Dsp::FilterDesign
              <Dsp::Butterworth::Design::HighPass <MAXORDER>, 1 >;
        m_bFilter = true;
        break;

    case HIGHPASS_BESSEL:
        filter = new Dsp::FilterDesign
              <Dsp::Bessel::Design::HighPass <MAXORDER>, 1 >;
        m_bFilter = true;
        break;

    case HIGHPASS_ELLIPTIC:
        filter = new Dsp::FilterDesign
              <Dsp::Elliptic::Design::HighPass <MAXORDER>, 1 >;
        m_bFilter = true;
        break;
    }

    Dsp::Params params;
    params[0] = f->samplingRate;    // sample rate
    if (f->order <= MAXORDER)
        params[1] = f->order;           // order
    else
        params[1] = MAXORDER;

    params[2] = f->cutoffFreq;      // cutoff frequency
    if(m_bFilter)
        filter->setParams (params);

}

Am I doing this right? The data that I want to filter is 1 channel sampled at 30kS/s. MAXORDER is set to 12 but the order that I am trying to run is 4. The cutoff frequency is usually around 400Hz. When I filter the data (with low frequency components and lots of high frequency noise), the output gets values very close to zero, nothing like what I would expect. Am I instantiating the filters right?
Also, here's the code for the actual processing:

Code: Select all (#)

if(m_pipeline.contains(PREPROC_DSP))
        {
            short int* pData[1];
            pData[0] = data;
            if(m_bDsp)
                m_dsp_f.filter->process(numSamples, pData);
        }

I should add that I added this

Code: Select all (#)

void process (int numSamples, short int* const* arrayOfChannels)

to the FilterDesign class in order to process short int data. Any ideas? I appreciate your work and any help you might provide!

[olalp]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

Hi all,

I just wanted to report a bug I found when using the following filter to do block processing of the data:

Code: Select all (#)

filter = new Dsp::FilterDesign
              <Dsp::Butterworth::Design::HighPass <MAXORDER>, 1 >;

I was getting a discontinuity at the boundary between buffers that had nothing to do with the data I was passing to it. By looking at it in more detail, I realized the filter was processing the second sample to whatever was in the adjacent memory location at the end of the array (it was off by one sample).
The problem was caused by Cascade::process(int numSamples, Sample* dest, StateType& state) const.
Instead of the line:

Code: Select all (#)

*dest++ = state.process (*dest, *this);

I modified it to

Code: Select all (#)

*dest = state.process (*dest, *this);
      ++dest;

And that fixed it.

[dinaiz]

Hi,

Nice project but I'm stumbling on a problem. The class SimpleFilter is template and has no non-template base which makes impossible to put it in a container.

I tried to make a container of Filter and use "new Dsp::SmoothedFilterDesign Dsp::RBJ::Design::LowPass,2,Dsp::DirectFormI(1024);" to populate it but my sound is completely garbled ! Does anyone knows how can I do ?

[Miles1981]

You should put either the type itself in the container, or use containers of unique_ptr<SimpleFilter>. That way, you don't try to copy elements.

[dinaiz]

Nevermind, this approach actually works. Here is the code in case it helps someone :

Code: Select all

Dsp::Filter* mFilter;
mFilter=new Dsp::SmoothedFilterDesign <Dsp::RBJ::Design::LowPass,2,Dsp::DirectFormI>(0);
	
	Dsp::Params p;
	p[0]=mSampleRate;
	p[1]=400; // Frequency
	p[2]=1.0; // Q
	mFilter->setParams(p);
        mFilter->process(buffer.getNumSamples(),buffer.getArrayOfChannels());