Solved

C++ - Lowpass FIR filtering using FFT convolution

Posted on 2012-04-08
15
2,634 Views
Last Modified: 2012-04-12
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
0
Comment
Question by:fjocke
[X]
Welcome to Experts Exchange

Add your voice to the tech community where 5M+ people just like you are talking about what matters.

  • Help others & share knowledge
  • Earn cash & points
  • Learn & ask questions
  • 6
  • 6
  • 3
15 Comments
 
LVL 7

Expert Comment

by:tampnic
ID: 37822843
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
 
LVL 3

Author Comment

by:fjocke
ID: 37823002
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
 
LVL 7

Expert Comment

by:tampnic
ID: 37823092
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
Percona Live Europe 2017 | Sep 25 - 27, 2017

The Percona Live Open Source Database Conference Europe 2017 is the premier event for the diverse and active European open source database community, as well as businesses that develop and use open source database software.

 
LVL 3

Author Comment

by:fjocke
ID: 37823683
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
 
LVL 27

Expert Comment

by:d-glitch
ID: 37824205
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
 
LVL 3

Author Comment

by:fjocke
ID: 37824623
Samplerate is 44100 and 16bit resolution
0
 
LVL 27

Expert Comment

by:d-glitch
ID: 37824752
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
 
LVL 27

Accepted Solution

by:
d-glitch earned 250 total points
ID: 37824844
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
 
LVL 27

Expert Comment

by:d-glitch
ID: 37824990
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
 
LVL 3

Author Comment

by:fjocke
ID: 37825010
I have no clue as how to create this "trapezoidal" array or where to get the weightning from.
0
 
LVL 27

Expert Comment

by:d-glitch
ID: 37825060
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
 
LVL 27

Expert Comment

by:d-glitch
ID: 37825065
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
 
LVL 7

Assisted Solution

by:tampnic
tampnic earned 250 total points
ID: 37825066
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
 
LVL 3

Author Comment

by:fjocke
ID: 37827557
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
 
LVL 3

Author Comment

by:fjocke
ID: 37836870
Yep Overlap Add solved the issue, however your knowledge learned me something new, so i'll accept your solutions with shared points :)
0

Featured Post

Get MongoDB database support online, now!

At Percona’s web store you can order your MongoDB database support needs in minutes. No hassles, no fuss, just pick and click. Pay online with a credit card. Handle your MongoDB database support now!

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

Article by: Nicole
This is a research brief on the potential colonization of humans on Mars.
The advancements in today's technology are unparalleled. Much of the technology that we have could not have been imagined twenty years ago. One of the latest additions to the list of technological advances is virtual reality. Virtual reality has an …
The viewer will learn how to user default arguments when defining functions. This method of defining functions will be contrasted with the non-default-argument of defining functions.
Viewers will learn key ranges in Sampler to make their sampled instruments sound more realistic Gather samples of various notes and drag them to Key Range panel: Set proper root key for each sample: Select all the samples with Command-A (or Ctrl…
Suggested Courses

615 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question