Solved

C++ - Lowpass FIR filtering using FFT convolution

Posted on 2012-04-08
15
2,491 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
  • 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
 
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
Highfive + Dolby Voice = No More Audio Complaints!

Poor audio quality is one of the top reasons people don’t use video conferencing. Get the crispest, clearest audio powered by Dolby Voice in every meeting. Highfive and Dolby Voice deliver the best video conferencing and audio experience for every meeting and every room.

 
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

Do You Know the 4 Main Threat Actor Types?

Do you know the main threat actor types? Most attackers fall into one of four categories, each with their own favored tactics, techniques, and procedures.

Join & Write a Comment

How to Win a Jar of Candy Corn: A Scientific Approach! I love mathematics. If you love mathematics also, you may enjoy this tip on how to use math to win your own jar of candy corn and to impress your friends. As I said, I love math, but I gu…
Technology opened people to different means of presenting information, but PowerPoint remains to be above competition. Know why PPT still works today.
The viewer will be introduced to the member functions push_back and pop_back of the vector class. The video will teach the difference between the two as well as how to use each one along with its functionality.
Viewers will learn how to include realistic velocity sensitivity to their Sampler instruments. Set the Vol<Vel parameter in the Filter/Global tab to your desired setting: Gather samples of hits of various intensity, and drag/drop into Velocity zon…

707 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

Need Help in Real-Time?

Connect with top rated Experts

18 Experts available now in Live!

Get 1:1 Help Now