C++ - Lowpass FIR filtering using FFT convolution

Hello fellow experts. I have on my spare time been trying to create a FIR Filter with FFT Convolution, with the following code:

#include "FIRFilter.h"


//Sinc function
double FIRFilter::sinc(double x)
{
	if(x == 0)
	{
		return 1;
	}
	return sin(PI*x)/(PI*x);
}

void FIRFilter::init()
{
	//Clear the filter buffer
	for(int i = 0; i < mfilterLength; i++)
	{
		this->mfilter[i] = 0;
	}

	//Clear the filterfft buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->mfilterFFT[i].im = 0;
		this->mfilterFFT[i].re = 0;
	}

	//Clear the audiofft buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->maudioFFT[i].im = 0;
		this->maudioFFT[i].re = 0;
	}

	//Clear the sortedfft buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->maudioSignal[i].im = 0;
		this->maudioSignal[i].re = 0;
	}

	//Clear the audiosignal buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->msortedFFTBuffer[i].im = 0;
		this->msortedFFTBuffer[i].re = 0;
	}
}

void FIRFilter::clean()
{
	//Clear the filter buffer
	for(int i = 0; i < this->mfilterLength; i++)
	{
		this->mfilter[i] = 0;
	}

	//Clear the filterfft buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->mfilterFFT[i].im = 0;
		this->mfilterFFT[i].re = 0;
	}

	//Clear the audiofft buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->maudioFFT[i].im = 0;
		this->maudioFFT[i].re = 0;
	}

	//Clear the sortedfft buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->maudioSignal[i].im = 0;
		this->maudioSignal[i].re = 0;
	}

	//Clear the audiosignal buffer
	for(int i = 0; i < this->mfftRes; i++)
	{
		this->msortedFFTBuffer[i].im = 0;
		this->msortedFFTBuffer[i].re = 0;
	}

	//De-allocate memory
	delete[] this->mfilter;
	delete[] this->mfilterFFT;
	delete[] this->maudioFFT;
	delete[] this->msortedFFTBuffer;
	delete[] this->maudioSignal;

}

FIRFilter::FIRFilter(int sampleRate)
{
	//Initilize varibles
	this->mSampleRate = sampleRate;
	this->mfftRes = 1024;
	this->mfilterLength = this->mfftRes;
	//this->maudioOverlapLength = this->mfftRes / 2;
	this->mSamplesNeeded = 0;

	//Create double buffers
	this->mfilter = new double[this->mfilterLength];
	//this->maudioOverlap = new double[maudioOverlapLength];

	//Initlize WDL FFT
	WDL_fft_init();

	//Create FFT Complex Buffers
	this->mfilterFFT = new WDL_FFT_COMPLEX[this->mfftRes];
	this->maudioFFT = new WDL_FFT_COMPLEX[this->mfftRes];
	this->maudioSignal = new  WDL_FFT_COMPLEX[this->mfftRes];
	this->msortedFFTBuffer = new  WDL_FFT_COMPLEX[this->mfftRes];

	this->mCutffOff = 1000;

	//Zero all the buffers
	this->init();

	//Create the Filter information
	for(int i=0; i < mfilterLength-1; i++ )
	{
		//Create lowpass filter
		double sincFilter = (2*this->mCutffOff/(double)this->mSampleRate)*sinc(2*this->mCutffOff*(i-((this->mfilterLength-1)/2.0))/(double)this->mSampleRate);
		//Apply a Blackman Window
		mfilter[i] = (0.42-0.5*cos((2*PI*i)/(double)(this->mfilterLength-1))+0.08*cos((4*PI*i)/(double)(this->mfilterLength-1))) * sincFilter;
	}

	//Copy the time domain filter data to the FFT Buffer
	for(int i = 0; i < mfilterLength; i++)
	{
		this->mfilterFFT[i].re = (WDL_FFT_REAL)this->mfilter[i];
		this->mfilterFFT[i].im = 0.; 
	}

	WDL_fft(this->mfilterFFT, this->mfftRes, false);
}

void FIRFilter::lowPassFilter(double *input, double *output1, double *output2, int nFrames)
{
	for (int s = 0; s < nFrames; ++s, ++input, ++output1,++output2)
	{
		if(this->mSamplesNeeded == this->mfftRes-1)
		{
			WDL_fft(this->maudioFFT, this->mfftRes, false);

			for (int i = 0; i < this->mfftRes; i++)
			{
				int j = WDL_fft_permute(this->mfftRes, i);
				
				this->msortedFFTBuffer[i].re = this->maudioFFT[j].re;
				this->msortedFFTBuffer[i].im = this->maudioFFT[j].im;
			} 

			//Pointwise multiplication of the filter and audio data in the frequency domain
			for(int i = 0; i < this->mfftRes; i++)
			{
				int j = WDL_fft_permute(this->mfftRes, i);
				//WDL_fft_complexmul(this->msortedFFTBuffer, this->mfilterFFT, this->mfftRes);
				//Calculate real with filter
				this->maudioSignal[j].re = (this->msortedFFTBuffer[i].re * this->mfilterFFT[j].re) + (this->msortedFFTBuffer[i].im * this->mfilterFFT[j].im);
				//Calculate imaginary with filter
				this->maudioSignal[j].im = (this->msortedFFTBuffer[i].re * this->mfilterFFT[j].im) - (this->msortedFFTBuffer[i].im * this->mfilterFFT[j].re);
			}

			//Invert FFT back to audio data
			WDL_fft(this->maudioSignal, this->mfftRes, true);

			this->mSamplesNeeded = 0;
		}
		else
		{
				this->mSamplesNeeded++;
		}

		this->maudioFFT[this->mSamplesNeeded].re = (WDL_FFT_REAL)*input;
		this->maudioFFT[this->mSamplesNeeded].im = 0.; 
		//Send the process data back into the output buffer
		*output1 =  this->maudioSignal[this->mSamplesNeeded].re / this->mfftRes;
		*output2 =  this->maudioSignal[this->mSamplesNeeded].re / this->mfftRes;
	}

}

FIRFilter::~FIRFilter(void)
{
	this->clean();
}

Open in new window


Here what i do, since *input only holds one sample, is store and cout up until i have 1024 samples which i then can process with FFT Convolution.

In the code i've created a sinc filter which seems to work, since all the frequency data gets cut off at the right frequency. But somewhere when i put the convolution back together into a audio signal again and send it back to output, it gets distorted.

Here you can hear the un-processed sample file:

http://soundcloud.com/syrou/normal

and here you can hear the processed sample with my code:

http://soundcloud.com/syrou/distorted

How do i fix this?

Kind Regards
LVL 3
fjockeAsked:
Who is Participating?

[Product update] Infrastructure Analysis Tool is now available with Business Accounts.Learn More

x
I wear a lot of hats...

"The solutions and answers provided on Experts Exchange have been extremely helpful to me over the last few years. I wear a lot of hats - Developer, Database Administrator, Help Desk, etc., so I know a lot of things but not a lot about one thing. Experts Exchange gives me answers from people who do know a lot about one thing, in a easy to use platform." -Todd S.

tampnicCommented:
I'm purely an analog guy when it comes to music processing (used to design and build valve amps) so don't have much knowledge of DSP. From a musical perspective the distortion sounds like transistor clipping which is a very square cutoff - have you graphed the waveform in detail to see if the sine wave tops are being squared off harshly? Saw wave distortion is also a possibility.

In general when you look at a waveform, anywhere where there is a sharp transition in amplitude (e.g. corners of a square wave, tips of a triangular wave) can potentially sound like ringing and distortion, unlike the smooth sine wave which sounds musical. Unfortunately I couldn't help you to remove those artifacts through DSP but I hope looking at the waveforms gives you some insight into the output of your FFT.

Cheers,
   Chris
0
fjockeAuthor Commented:
Thanks for putting in some effort atleast to share your knowledge, but as you know this does not really explain what's going on.
0
tampnicCommented:
Maybe some kind of interpolation between samples needs to happen to smooth transitions out? Does increasing or reducing the resolution affect the distortion noticeably? How about doing one conversion, then do another but starting from an offset of 1/2 the sample length of the original, then synchronising and adding the two waves together?

Cheers,
   Chris
0
Python 3 Fundamentals

This course will teach participants about installing and configuring Python, syntax, importing, statements, types, strings, booleans, files, lists, tuples, comprehensions, functions, and classes.

fjockeAuthor Commented:
What you are talking about does sound something like the method Overlap add/save http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method which i've heard is what this code needs, however i can't see why it's needed.

And neither how to implement it.

About increasing resolution, yes it does indeed seem to affect it very much, for the worse.
0
d-glitchCommented:
What are the sample rate and bit resolution?

I think tampnic is probably on to something.  The high frequencies you are filtering out may be smoothing out the transitions between your sample sets.
If you filter them out you get the clicking.  You can check this by plotting the analog data at the transitions before and after filtering.

I think the usual way to handle this is to use overlapping trapezoidal filters for each sample set.  This allows you to make a gradual transition between samples rather than an abrupt one.
This another version of what tampnic suggested.  You probably can get by with much less than 50% overlap.
0
fjockeAuthor Commented:
Samplerate is 44100 and 16bit resolution
0
d-glitchCommented:
44100/1024  ==>  43 FFT sample sets per second

That can mean 43 clicks per second if you don't take care of the transitions between FFT sample sets.
And that is likely the noise you are hearing on the distorted sample.
0
d-glitchCommented:
To employ a trapezoidal filter you would take your array of 1024 samples and multiply it by another array that has a trapezoidal weighting before you do the FFT.

   [ 0/32   1/32   2/32   ...   31/32   1   1   1   1   ...   1   1   31/32   30/32   ...   2/32   1/32   0/32 ]

Note how the amplitude at the beginning and end of the FFT array are zero, so you shouldn't have any clicks.

You would have to let the arrays overlap by 32 samples (or 64, or 100, or 128).  You have to experiment to see how much overlap you need.
This increases the processing time slightly, but it should fix your problem.
0

Experts Exchange Solution brought to you by

Your issues matter to us.

Facing a tech roadblock? Get the help and guidance you need from experienced professionals who care. Ask your question anytime, anywhere, with no hassle.

Start your 7-day free trial
d-glitchCommented:
On the array multiplication above:  You only have to modify the first 31 and last 31 elements of each segment.

The effect of the trapezoidal multiplications is that one segment fades out while the next one is fading in.  This should eliminate the abrupt transitions/clicks.
0
fjockeAuthor Commented:
I have no clue as how to create this "trapezoidal" array or where to get the weightning from.
0
d-glitchCommented:
These are the values:  

 [ 0/32   1/32   2/32   ...   31/32   1   1   1   1   ...   1   1   31/32   30/32   ...   2/32   1/32   0/32 ]

The first 31 elements are      0/32 = 0.0000      1/32 =0 .0313     2/32 = 0.0625   ...   31/32 = 0.9688

The next 963 elements are   1.0000  <==  This means don't change the elements of the signal array.

The last 31 elements are     31/32 = 0.9688    30/32 = 0.9375   ...   1/32 =0 .0313     0/32 = 0.0000
0
d-glitchCommented:
The weighting means ramp from zero to one at the beginning ==>  fade in

and  ramp from one to zero at the end ==> fade out.
0
tampnicCommented:
Why not push your output through a software spectrum analyser - if you see peaks in the 43Hz, 86Hz, 172Hz (maybe also higher harmonics of 43Hz) regions then the distortion can be attributed to the transition between samples and dealt with using d-glitch's method.

Cheers,
   Chris
0
fjockeAuthor Commented:
I'll be testing overlap add out before, after i read some theory it seems like overlap add is what's really missing here. I'll get back to you with the results, and accepting the solutions
0
fjockeAuthor Commented:
Yep Overlap Add solved the issue, however your knowledge learned me something new, so i'll accept your solutions with shared points :)
0
It's more than this solution.Get answers and train to solve all your tech problems - anytime, anywhere.Try it for free Edge Out The Competitionfor your dream job with proven skills and certifications.Get started today Stand Outas the employee with proven skills.Start learning today for free Move Your Career Forwardwith certification training in the latest technologies.Start your trial today
C++

From novice to tech pro — start learning today.