FIR Filter for decimation

[andy_k]

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 KVR community!
I wish to all of you a happy 2018!

I am trying to use an FIR filter for reducing aliasing in JUCE
The main idea is based on this info.
http://www.willpirkle.com/Downloads/AN- ... yPhase.pdf (http://www.willpirkle.com/Downloads/AN-1MultiRateandPolyPhase.pdf)
and part of the implementation was based on this topic

https://forum.juce.com/t/basic-fir-filtering/24477/8 (https://forum.juce.com/t/basic-fir-filtering/24477/8)

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"
https://discourse-cdn-sjc1.com/business ... b009a6.png (https://discourse-cdn-sjc1.com/business/uploads/juce/original/2X/d/d64b82f867485f24e00519d14cae03506eb009a6.png)

Ande here is with it:

https://discourse-cdn-sjc1.com/business ... c724d3.png (https://discourse-cdn-sjc1.com/business/uploads/juce/original/2X/f/f73dc92b72ab0117b38bc7f05a7842a4bec724d3.png)

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

JUCE CODE:

FIR_Filter.h

Code: Select all (#)

#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

Code: Select all (#)

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

Code: Select all (#)

   CFIR_Filter LPF_Filter;

In PluginProcessor.cpp

Code: Select all (#)

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;
           
            
        }
    }
}

[matt42]

Your "oversampled" wave form isn't oversampled at all. All you have done is filter the high frequency content from an aliased wave form.

Code: Select all

angleIncrement= (2 * M_PI * frequency ) / (samplingRate*nOversamplingRatio) ;
//......
 for (int i=0; i<nOversamplingRatio; i++ )
currentAngle+=angleIncrement;

Can you see it? currentAngle increments in exactly the same way whatever nOversamplingRatio. In other words your ratio has no effect. You are generating the same wave form in all cases.

and this

Code: Select all

float *monoBuffer = new float[ buffer.getNumSamples()];//Where is the extra space for an oversampled waveform?

To oversample you need to do something like:

Code: Select all

float *monoBuffer = new float[ buffer.getNumSamples()*nOversampleRatio];

Then you calculate the saw wave at the new sample rate:

Code: Select all

for(int sample = 0; sample < buffer.getNumSamples()*nOverSampleRatio; ++sample)//taking oversampling ratio into account
        
            {
            
                monoBuffer [sample] = (float) 1 - (1  / M_PI * currentAngle ); //Saw

The remember to fix the current angle calculation:

Code: Select all

angleIncrement= (2 * M_PI * frequency ) / (samplingRate*nOversamplingRatio) ;
//......
 //for (int i=0; i<nOversamplingRatio; i++ ) don't want this line
currentAngle+=angleIncremen;

and you'll need to make sure you're FIR is set to remove frequencies above 16000, or whatever cut off you chose at the oversampled sample rate.

And then you throw away the unwanted samples, for example every second sample if you upsamples by x2.

It seems you don't really understand what oversampling is. You might want to spend some time swatting up and understand what you are trying to do.

Also, allocation of buffers in the process loop is a bad idea due to unbounded execution time, though you can get away with it if you are just testing things out. Not deleting the buffer and then reallocating every process call is even worse!

Hope that helps

Edit: Oversampling a saw will work, but to push the aliasing noise down to low levels requires very high oversampling ratio. A more efficient method would one of the blit/blep methods, like polyblep for example.

[S0lo]

In an nut shell,

1. Calculate every thing at a higher rate. like 2X, 4X or 8X etc..... Lets call it nX
2. Filter with a FIR at half the original sample rate. NOT half the oversampled sample rate (nX).
3. Decimate. that is THROW AWAY every n-1 samples.

One note here. Do not use "zero stuffing" or any other interpolation method. These methods are meant for sample rate conversion of already none aliased signals. Not for oversampling to remove aliasing.

Thats as far as I know.

[aciddose]

I've said this before, it seems to be a common misunderstanding for beginners in DSP:

Oversampling is always the worst possible solution to any aliasing related problem due to discontinuities with infinite harmonic spectra. In other words it's the worst performing, lowest quality and most expensive method.

Oversampling is somewhat of a misnomer which is what leads to this confusion: you should always sample at the rate required to bound the entire band of the signal being sampled. So for example if you apply a non-linear process which increases the bandwidth such as "ringmod" which is 2nd order (quadratic) = double the bandwidth required, sampling at 2x to effectively sample the complete resulting band is the solution. So the "over" refers to the base or target rate which implies a filtering step before decimation to the target rate. Otherwise the proper term is simply "sampling".

This technique doesn't apply when you have a signal with an infinite number of harmonic partials. In order to reduce the level of aliased harmonics you may need to sample at an infinite rate in the case where you have amplitude = 1. In cases such as a first order discontinuity with 1/N harmonic amplitudes you'll only gain 6.02 dB every time you double the sample rate. In such cases you must find a way to efficiently filter or otherwise eliminate or reduce those partials with some other method before the signal is sampled.

[andy_k]

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 and thank you for the answers.
I know that I can calculate the saw using the oversampling ratio like this:

Code: Select all (#)

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 ) / (samplingRatenOversamplingRatio) ;
currentAngle+=angleIncrement;
if(currentAngle>(2 * M_PI)) currentAngle-= (2 * M_PI);
}

This is very inefficient since it needs to:

Calculate the Saw at a high sample rate.
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 ... or—part-3/ (http://www.earlevel.com/main/2012/05/09/a-wavetable-oscillator%E2%80%94part-3/)

[matt42]

Will Pirkle is talking about resampling an existing signal by zero padding then filtering. This method isn't applicable to what you're doing. You are not resampling a band limited saw, after all, but trying to generate one in the first place.

Edit: the quote you mention actually does refer to decimation. The theoretical zeroes you mention though have nothing to do with the decimation as described in Will's PDF. The efficiency of his decimation simply comes down to the fact that with a FIR we only need to calculate the filter response for the retained output samples, rather than every output sample. For 2x decimation we only run the filter on every second sample for example. No point calculating samples that will be discarded after all. This point might be less than obvious and perhaps be obfuscated by description of the polyphase implementation. You still need to calculate the wave table at the higher sample rate for this to work.

A band limited wave table will most likely be a lot more efficient than an over sampled saw wave. You should give that a shot.

Also you might find some polyblep source code if you google it. I think Will Pirkle describes this method in his synth book if I recall correctly

[matt42]

Some of what Will Pickle says could actually be very misleading. With regards the decimation:

The result is
that we will decimate before we filter, running our polyphase filters at the original (slower) sample rate.

This is really misleading the samples are not discarded before filtering. Look at the algorithm:

Input sequence on the PDF is a0 a1 a2 b1 b2 c1 c2

The output is described as:

Code: Select all

y(n) = a0h(0) + b0h(3) + c0h(6)
 + a1h(1) + b1h(4) + c1h(7) 
 + a2h(2) + b2h(5) + c2h(8)

Or In polyphase form

Code: Select all

y(n) = y0(n) + y1(n) + y2(n)
where:
y0(n) = a0h(0) + b0h(3) + c0h(6)
y1(n) = a1h(1) + b1h(4) + c1h(7)
y2(n) = a2h(0) + b2h(5) + c2h(8)

so calculating y(n) requires all inputs. They haven't simply been decimated before filtering. I think it's terrible wording on his part.

In this case 3x down sample the signal is decimated before filtering as he says, but decimated into 3 signals which between them contain all the original samples. It isn't a magic trick to just discard information and then filter

[matt42]

. But the combination of oversampling factor and filter order makes it unusable for larger values (e.g. 16x oversampling, filter order 640).

While 16x will be inefficient a 640 order filter is huge! Is that a typo?

[S0lo]

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

Code: Select all

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 ) / (samplingRatenOversamplingRatio) ;
currentAngle+=angleIncrement;
if(currentAngle>(2 * M_PI)) currentAngle-= (2 * M_PI);
}

This is very inefficient since it needs to:

You are right, it is very inefficient for waveform generation, Google for Minblep and polybleb.

[andy_k]

While 16x will be inefficient a 640 order filter is huge! Is that a typo?

It is not a typo. I was using only 1/16 of the coefficients (40 coeff).