FIR Filter for decimation


#1

Hi JUCE community!
I wish to all of you a happy 2018!

I am trying to use an FIR filter for reducing aliasing.
The main idea is based on this info.

http://www.willpirkle.com/Downloads/AN-1MultiRateandPolyPhase.pdf

and part of the implementation was based on this topic

Here is the problem. The filter seems to cut the desired frequencies, but the aliasing is not reduced at all.

Here is a sawtooth at 440 Hz without "oversampling"


Ande here is with it:

Any ideas will be very appreciated!!!
Thanks a lot!

FIR_Filter.h

#ifndef FIR_Filter_h
#define FIR_Filter_h
#include <iostream>

using namespace std;

#endif /* FIR_Filter_h */

class CFIR_Filter
{
public:
    CFIR_Filter (void);
    virtual ~CFIR_Filter(void);
   
    
    void init ();
    
    float processFilterBlock (float x);
    int getOverSamplingRatio() {        return m_nOverSamplingRatio;} ;
    
    inline void setOverSamplingRatio(int m_nOverSamplingRatio) {this->m_nOverSamplingRatio = m_nOverSamplingRatio;
         }
    
    inline int getOffset() { return m_Offset;} 
    
    inline void setOffset(int m_Offset) {this->m_Offset = m_Offset;}
    
    inline int getIRLength() {
     
        return m_nIRLength;} 
    
    inline void setIRLength(int m_nIRLength) {this->m_nIRLength = m_nIRLength;}
    
    inline void setDelayLength() {this->m_nDelayLineSize = this->m_nIRLength /this->m_nOverSamplingRatio;
       }


    inline void setDelayLineIndex() {this->DelayLineIndex = this->m_nDelayLineSize-1;
      }


    
private:
    // Set OverSampling ratio
    int m_nOverSamplingRatio = 16;
    // Set the IRLenght here
    int m_nIRLength = 640;
    
    // The IR Array here
    float *m_FilterImpulseResponse;
    
    // The Delay Line Array here
    int m_nIRIndex = 0;
    int m_nDelayLineSize = m_nIRLength/m_nOverSamplingRatio;
    int m_Offset = m_nOverSamplingRatio-1;
    
    
    float *m_DelayLine;
    int DelayLineIndex= 0;
    int DelayLineLoopCounter=0;
   
    
};

FIR_Filter.cpp

CFIR_Filter::CFIR_Filter (void)
{
    //init array impulse response
    m_FilterImpulseResponse = new float[m_nIRLength];
    memset(m_FilterImpulseResponse, 0, m_nIRLength*sizeof(float));
    
    // Add your impulse response here
    
	 m_FilterImpulseResponse[0] = -0.0000000396927433;
     m_FilterImpulseResponse[1] = -0.0000000578565071;
     m_FilterImpulseResponse[2] = -0.0000000987668898;
.....
	 m_FilterImpulseResponse[637] = -0.0000000987668898;
     m_FilterImpulseResponse[638] = -0.0000000578565071;
     m_FilterImpulseResponse[639] = -0.0000000396927433;
    
    
    
    //Init array Delay Line, fill with zeros
    m_DelayLine = new float[m_nIRLength/m_nOverSamplingRatio];
    memset(m_DelayLine, 0, m_nIRLength*sizeof(float));
    for (int i=0; i<m_nIRLength; i++) m_DelayLine[i] = 0.0;
  
    
}

  CFIR_Filter::~CFIR_Filter(void)
{
    if(m_DelayLine) delete  m_DelayLine;
    if(m_FilterImpulseResponse) delete m_FilterImpulseResponse;
}
void CFIR_Filter::init()
{
    
}

float CFIR_Filter::processFilterBlock(float x)
{
    

    m_DelayLine[DelayLineIndex] = x;
    DelayLineLoopCounter = DelayLineIndex;
    float y = 0;

    for (int i = m_Offset; i < m_nIRLength; i+=m_nOverSamplingRatio)  // for each tap
    {
        y+= m_DelayLine[DelayLineLoopCounter] * m_FilterImpulseResponse[i];
            if (--DelayLineLoopCounter < 0)
            DelayLineLoopCounter += m_nDelayLineSize;              // wrap buffer index
    }
    if (++DelayLineIndex >= m_nDelayLineSize)
        DelayLineIndex -= m_nDelayLineSize;              // wrap buffer index
 
    return y*m_nOverSamplingRatio;

}

In PluginProcessor.h

   	CFIR_Filter LPF_Filter;

In PluginProcessor.cpp

OscAudioProcessor::OscAudioProcessor()
#ifndef JucePlugin_PreferredChannelConfigurations
     : AudioProcessor (BusesProperties()
                     #if ! JucePlugin_IsMidiEffect
                      #if ! JucePlugin_IsSynth
                       .withInput  ("Input",  AudioChannelSet::stereo(), true)
                      #endif
                       .withOutput ("Output", AudioChannelSet::stereo(), true)
                     #endif
                       )
#endif

{
    
    LPF_Filter.setOverSamplingRatio(nOversamplingRatio);
    LPF_Filter.setOffset(0);
    LPF_Filter.setIRLength(640);
    LPF_Filter.setDelayLength();
    
}

void OscAudioProcessor::processBlock (AudioSampleBuffer& buffer, MidiBuffer& midiMessages)
{
 
    const int totalNumInputChannels  = getTotalNumInputChannels();
    const int totalNumOutputChannels = getTotalNumOutputChannels();
    
 for (int i = totalNumInputChannels; i < totalNumOutputChannels; ++i)
        buffer.clear (i, 0, buffer.getNumSamples());

    // This is the place where you'd normally do the guts of your plugin's
    // audio processing...
   
    
    float *monoBuffer = new float[ buffer.getNumSamples()];
    for(int sample = 0; sample < buffer.getNumSamples(); ++sample)
        
            {
            
                monoBuffer [sample] = (float) 1 - (1  / M_PI * currentAngle ); //Saw
            if(aliasProcessor)
            {
                monoBuffer [sample] =LPF_Filter.processFilterBlock((float) monoBuffer [sample]);
            }
     
            // Oversample
             angleIncrement= (2 * M_PI * frequency ) / (samplingRate*nOversamplingRatio) ;
            

            if (oscSwitch)
            {
                frequency+=0.01;
                if(frequency> (samplingRate/2))
                {
                    frequency=0;
                    oscSwitch=false;
                }
            }
     #if 0
            currentAngle+=angleIncrement*nOversamplingRatio;
#else
            for (int i=0; i<nOversamplingRatio; i++ )
                currentAngle+=angleIncrement;
#endif
        
        if(currentAngle>(2 * M_PI)) currentAngle-= (2 * M_PI);
        }  
    for (int channel = 0; channel < totalNumInputChannels; ++channel)
    //*
    {
        float* channelData = buffer.getWritePointer (channel);
              // ..do something to the data..
        for (int sample = 0; sample < (buffer.getNumSamples()); sample++)
        {
            channelData[sample] =  monoBuffer [sample]*onOff;
        }
    }
}


#2

Can you fix the code formatting so we can get an easier look at things? A code block should have three " ``` " (back ticks) at the beginning and end to do the whole block, i.e.

\```
paste unmodified code here
\```

(ignore the backslashes)


#3

Done! And thank you Jonathonracz!


#4

The second chart shows the band limiting, but neither shows the evidence of aliasing - I can’t see any inharmonic frequencies introduced. Can you zoom in on the higher frequencies and/or share some audio files?


#5

Hi!
Here are some sound files:

440 Saw without the filter

440 Saw with the filter

Sweep 0.1Hz to 22050 Hz without the filter

Sweep 0.1Hz to 22050 Hz with the filter

Thanks again!
Andy


#6

OK, that’s definitely aliasing heaps and the filter isn’t doing anything to prevent it.

I just noticed that you are synthesising your saw before you upsample. Think about that for a moment - you’re creating the aliases inside your wanted audio band before you even hit the oversampling code.

What you need to do is synthesise at the oversampled rate and then downsample the result. So if you are 16x oversampling, then generate 16 times as many samples and then decimate. You might also need to stop multiplying the angleIncrement by nOversamplingRatio.


#7

Hi and thank you for the answer.
I know that I can calculate the saw using the oversampling ratio like this:

```
float *monoBuffer = new float[ buffer.getNumSamples()*nOversamplingRatio];
for(int sample = 0; sample < buffer.getNumSamples()nOversamplingRatio; ++sample)
{
monoBuffer [sample] = (float) 1 - (1 / M_PI * currentAngle ); //Saw
angleIncrement= (2 * M_PI * frequency ) / (samplingRate
nOversamplingRatio) ;
currentAngle+=angleIncrement;
if(currentAngle>(2 * M_PI)) currentAngle-= (2 * M_PI);
}
```

This is very inefficient since it needs to:

  1. Calculate the Saw at a high sample rate.
  2. Filter using all the coefficients of the FIR filter.

By changing my code to this and using all the coefficients I got a lot of alias reduction. But the combination of oversampling factor and filter order makes it unusable for larger values (e.g. 16x oversampling, filter order 640).

As quoted by Will Pirkle:

“You can see this is horribly inefficient since we convolve on each iteration of the loop, even though we discard M-1 of these convolutions. We can certainly fix that by only convolving once after M samples have been acquired.”

My attempt was to crack down the filter in subfilters and just use one of this parts for the final convolution. That is why I said my wave was “oversampled”. In fact, was like the original wave with theoretical zeros in between.

Anyway, I am not happy with the results of the experiment. Maybe I want to try something like this:
http://www.earlevel.com/main/2012/05/09/a-wavetable-oscillator—part-3/


#8

Haha - never said it would be efficient :slight_smile:

You might also want to check out BLITs and BLEPs if you haven’t already


#9

Thanks Andrew!
I already tried this!
He also has this specific project for PolyBLEP


#10

Well, how did it turn out? Got rid of the aliasing?


#11

That’s a short FIR filter (9 taps) and my guess is that it’s unlikely to be sufficient for reducing aliasing significantly, especially since the spectrum of the saw falls off slowly.

You can use the scipy.signal library to design FIR filters with specific responses using the remez function. The problem with a longer filter, however, is that it’s a longer filter and that introduces latency, especially if you’re going with a linear phase FIR filter (as shown in your code which uses L/R symmetrical coefficients). Scipy.signal does have tools for converting a linear phase filter to a minimum phase filter which would help with the latency, but this is still overkill for an oscillator and you’ll be better off going with one of the established antialiasing methods.


#12

The FIR filter is now reducing considerably the aliasing. The PolyBleep worked well in WDL-OL (not tried this yet in JUCE). The proposed method in http://www.earlevel.com/main/2012/05/09/a-wavetable-oscillator—part-3/3 works very well.
Now I am trying minBlep…


#13

Thank you very much for the replies!