• Status: Solved
  • Priority: Medium
  • Security: Public
  • Views: 2854
  • Last Modified:

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
0
fjocke
Asked:
fjocke
  • 6
  • 6
  • 3
2 Solutions
 
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
Free Tool: SSL Checker

Scans your site and returns information about your SSL implementation and certificate. Helpful for debugging and validating your SSL configuration.

One of a set of tools we are providing to everyone as a way of saying thank you for being a part of the community.

 
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
 
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
Question has a verified solution.

Are you are experiencing a similar issue? Get a personalized answer when you ask a related question.

Have a better answer? Share it in a comment.

Join & Write a Comment

Featured Post

Cloud Class® Course: SQL Server Core 2016

This course will introduce you to SQL Server Core 2016, as well as teach you about SSMS, data tools, installation, server configuration, using Management Studio, and writing and executing queries.

  • 6
  • 6
  • 3
Tackle projects and never again get stuck behind a technical roadblock.
Join Now