C / C++ FFT (Fast Fourier Transform) Based Signal Filter
| IT | No Comments
The following design is a FFT (Fast Fourier Transform) based signal filter developed in C / C++.
You can use this type of filter to amplify or dampen very specific bands.
Alternatively you could use it as a band pass, low pass, or high pass filter by simply setting coefficient ranges to zero.
And yet another great application is to use it as a pitch shifter, for example, with digital audio, and keep the timing the same. I.e. change the pitch without playing the sound faster or slower.
//**************************** class DSP_FFT_Filter ***************************
// Defines are filter based on Fast Fourier Transform
// transferFunctionCoef and numberOfElements define the transfer function.
// these coefficients are multiplied with the complex vectors returned by FFT
// Note: coeff vectors appear in reverse order, ie.
// transferFunctionCoef[numberOfElements-1] is the base frequency,
// transferFunctionCoef[numberOfElements-2] is first harmonic.
class DSP_FFT_Filter : public DSP_ProcessorBase<double, double>
{
public:
DSP_FFT_Filter(float *transferFunctionCoef, int numberOfElements, int blockSize);
virtual ~DSP_FFT_Filter()
{
if (m_pTransferFunctionCoef) delete[] m_pTransferFunctionCoef;
#ifdef _DEBUG
if (m_pFFTCoefficients) delete[] m_pFFTCoefficients;
#endif
}
virtual char* GetProcessorName(void){return "DSP_FFT_Filter";}
virtual bool ProcessSamples(const DSP_Samples &inputSamples);
virtual bool ProcessFFTCoeff(float *data, unsigned long nn);
float *m_pTransferFunctionCoef;
int m_NumberOfElementsInTransferFunction;
#ifdef _DEBUG
float *m_pFFTCoefficients;
int m_NumberOfCoeffients;
#endif
};
The header above is simple, implementation follows below:
//************************ class DSP_FFT_Filter *******************************
// returns true if number is a power of 2
bool IsPowerOfTwo(unsigned number)
{
bool result = false;
unsigned mask = 1;
for (int i = 0; i < 32; i++)
{
if ((number & mask) == number) // if one bit is set only, return true
{
result = true;
break;
}
else
mask <<= 1; // shift mask and try again...
}
return result;
}
// returns bit number of MSB, ex: number = 8 -> retVal = 3
int GetHighestBit(unsigned number)
{
unsigned mask = 0x80000000;
for (int i = 31; i >= 0; i--)
{
if (number & mask)
return i;
else
mask >>= 1;
}
return 0;
}
// bool four1(float data[], unsigned long nn, int isign);
// Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as 1; or replaces
// data[1..2*nn] by nn times its inverse discrete Fourier transform, if isign is input as - 1.
// data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn MUST
// be an integer power of 2 (this is not checked for!).
// data[0] is NOT used
// data[1] contains the real part of element #1
// data[2] contains the imaginary part of element #1
// and so on ...
// in case of real data, nn MUST BE A POWER OF TWO!!!
// returns true if check was ok
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
bool four1(float data[], unsigned long nn, int isign)
{
bool result = false;
assert(IsPowerOfTwo(nn));
IF_NLOG(IsPowerOfTwo(nn), DSPLIB_ET_FFT_INPUT_MUST_BE_POWER_OF_TWO, "four1")
{
unsigned long n, mmax, m, j, istep, i;
double wtemp, wr, wpr, wpi, wi, theta;
float tempr, tempi;
n = nn << 1;
j = 1;
for (i = 1; i < n; i += 2)
{ // This is the bit-reversal section of the routine.
if (j > i)
{
SWAP(data[j],data[i]); //Exchange the two complex numbers.
SWAP(data[j+1],data[i+1]);
}
m=nn;
while (m >= 2 && j > m)
{
j -=m;
m >>= 1;
}
j +=m;
}
//Here begins the Danielson-Lanczos section of the routine.
mmax = 2;
while (n > mmax)
{ // Outer loop executed log 2 nn times.
istep = mmax << 1;
theta = isign * (6.28318530717959/mmax); //Initialize the trigonometric recurrence.
wtemp = sin (0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
for (m=1;m<mmax;m+=2)
{ //Here are the two nested inner loops.
for (i=m;i<=n;i+=istep)
{
j = i + mmax; //This is the Danielson-Lanczos for-mula:
tempr = float(wr * double(data[j]) - wi * double(data[j+1]));
tempi = float(wr * double(data[j+1]) + wi * double(data[j]));
data[j] = data[i] - tempr;
data[j+1] = data[i+1] - tempi;
data[i] += tempr;
data[i+1] += tempi;
}
wr = (wtemp = wr) * wpr - wi * wpi + wr;// Trigonometric recurrence.
wi = wi * wpr + wtemp * wpi + wi;
}
mmax=istep;
}
result = true;
}
return result;
}
DSP_FFT_Filter::DSP_FFT_Filter(float *transferFunctionCoef, int numberOfElements, int blockSize)
: DSP_ProcessorBase<double, double>(blockSize)
{
m_pTransferFunctionCoef = NULL;
m_NumberOfElementsInTransferFunction = 0;
// copy over transferFunction coefficients
IF_NLOG(transferFunctionCoef, DSPLIB_ET_FFT_FILTER_TRANSFER_FUNC_COEFS_MISSING, "")
{
m_pTransferFunctionCoef = new float[numberOfElements];
m_NumberOfElementsInTransferFunction = numberOfElements;
for (int i = 0; i < numberOfElements; i++)
{
m_pTransferFunctionCoef[i] = transferFunctionCoef[i];
}
}
#ifdef _DEBUG
m_pFFTCoefficients = NULL;
m_NumberOfCoeffients = 0;
#endif
}
bool DSP_FFT_Filter::ProcessSamples(const DSP_Samples &inputSamples)
{
bool result = false;
// ensure input is power of two...if not we will need to add some extra
// padding
// so figure out if input length is power of two, if not compute extra space needed
int extraSpaceRequired = 0;
if (IsPowerOfTwo(inputSamples.m_NumberOfSamples) == false)
{
extraSpaceRequired = (1 << (GetHighestBit(inputSamples.m_NumberOfSamples) + 1)) - inputSamples.m_NumberOfSamples;
}
DSP_Samples *outputSamples = (DSP_Samples*)m_pOutputSamples;
IF_NLOG(outputSamples, DSPLIB_ET_OUTPUTSAMPLE_OBJ_MISSING, "")
{
outputSamples->ReallocWithoutCopy(inputSamples.m_NumberOfSamples + extraSpaceRequired);
assert(outputSamples->m_BlockSize >= inputSamples.m_NumberOfSamples + extraSpaceRequired);
IF_NLOG(outputSamples->m_BlockSize >= inputSamples.m_NumberOfSamples + extraSpaceRequired, DSPLIB_ET_OUTPUTSAMPLE_OBJ_BLOCKSIZE_TOO_SMALL, "")
{
outputSamples->m_NumberOfSamples = inputSamples.m_NumberOfSamples + extraSpaceRequired;
outputSamples->m_SampleRate = inputSamples.m_SampleRate;
assert(IsPowerOfTwo(inputSamples.m_NumberOfSamples + extraSpaceRequired));
IF_NLOG(IsPowerOfTwo(inputSamples.m_NumberOfSamples + extraSpaceRequired), DSPLIB_ET_FFT_INPUT_MUST_BE_POWER_OF_TWO, "")
{
int originalNumberOfSamples = inputSamples.m_NumberOfSamples + extraSpaceRequired;
int numberOfSamples = inputSamples.m_NumberOfSamples;
// allocate temp buffer
float *tmpData = new float[(numberOfSamples + extraSpaceRequired) * 2 + 2];
memset(tmpData, 0, sizeof(float) * ((numberOfSamples + extraSpaceRequired) * 2 + 2));
IF_NLOG(tmpData, ET_OUT_OF_MEMORY, "")
{
// convert samples into floats first...
for (int i = 0; i < (numberOfSamples + extraSpaceRequired) && i < inputSamples.m_NumberOfSamples; i++)
{
tmpData[1 + i * 2] = float(inputSamples.m_pSamples[i]);
}
// run FFT
four1(tmpData, numberOfSamples + extraSpaceRequired, 1);
#ifdef _DEBUG
// if in debug mode, store FFT results for later review
if (m_pFFTCoefficients) delete[] m_pFFTCoefficients;
m_NumberOfCoeffients = numberOfSamples + extraSpaceRequired;
m_pFFTCoefficients = new float[(numberOfSamples + extraSpaceRequired) * 2];
for(int jj = 0; jj < numberOfSamples + extraSpaceRequired; jj++)
{
m_pFFTCoefficients[jj * 2] = tmpData[1 + jj*2];
m_pFFTCoefficients[jj * 2 + 1] = tmpData[1 + jj*2 + 1];
}
#endif
IF_NLOG(ProcessFFTCoeff(tmpData, numberOfSamples + extraSpaceRequired), DSPLIB_ET_FFT_FILTER_PROCESS_COEFS_FAILED, "")
{
// run IFFT
four1(tmpData, numberOfSamples + extraSpaceRequired, -1);
// convert samples back to doubles and normalize...and we're done!
for (i = 0; i < numberOfSamples + extraSpaceRequired; i++)
{
outputSamples->m_pSamples[i] = double(tmpData[1+i*2]) / double(originalNumberOfSamples);
}
outputSamples->m_NumberOfSamples = originalNumberOfSamples;
result = true;
}
delete[] tmpData;
}
}
}
}
return result;
}
bool DSP_FFT_Filter::ProcessFFTCoeff(float *tmpData, unsigned long nn)
{
bool result = true;
// apply transfer function coef by MULTIPLYING (modify filter attenuation)
for (int i = 0; i < m_NumberOfElementsInTransferFunction && i < int(nn); i++)
{
// and write back result
tmpData[1 + i*2] = tmpData[1 + i*2] * m_pTransferFunctionCoef[i];
tmpData[1 + i*2 + 1] = tmpData[1 + i*2 + 1] * m_pTransferFunctionCoef[i];
}
return result;
}