Solved

C++ - Lowpass FIR filtering using FFT convolution

Posted on 2012-04-08
15
2,559 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
Online Training Solution

Drastically shorten your training time with WalkMe's advanced online training solution that Guides your trainees to action. Forget about retraining and skyrocket knowledge retention rates.

 
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

SharePoint Admin?

Enable Your Employees To Focus On The Core With Intuitive Onscreen Guidance That is With You At The Moment of Need.

Question has a verified solution.

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

Suggested Solutions

Title # Comments Views Activity
No module found pypyodbc, 3 45
Outlook 13 84
Free Video Recorder 6 23
Wikipedia trigonometric functions main diagram labeled incorrectly 4 31
I have a Synology DS212+ NAS.  These are not only great for backup and normal NAS stuff, but also for delivering media throughout your home or LAN via DLNA.  I copied my whole audio collection from iTunes over to the box, but couldn't figure out how…
Go is an acronym of golang, is a programming language developed Google in 2007. Go is a new language that is mostly in the C family, with significant input from Pascal/Modula/Oberon family. Hence Go arisen as low-level language with fast compilation…
Viewers will learn how to turn a Live Set into a compressed Live Pack file, and how to install Live Packs. Make: File > Collect All And Save: File > Manage Files: Click Manage Project: Click Create Pack: Select save location: Install: Doub…
Viewers will get an overview of how to make and use Drum Racks in Ableton Live. Load new Drum Rack into empty MIDI track: Fill rack with audio samples: Re-arrange sample slots as necessary: Adjust parameters of each slot to tailor each sound a…

756 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