/*
July 17, 2015
Iowa Hills Software LLC
http://www.iowahills.com
If you find a problem with this code, please leave us a note on:
http://www.iowahills.com/feedbackcomments.html
This code is meant to serve as an example for the code needed generate
FIR coefficients via the frequency sampling method. This is not stand alone
code. We have not supplied definitions for the global variables, but they
should be self explanatory.
We show 4 pieces of code here.
1. SimulateAnalogProto() will evaluate a filter polynomial, such as a Butterworth.
It shows how to do the frequency warping and generate a low pass, high pass,
band pass, and notch response from the polynomial.
2. CalcWinPhase() will set the phase for a linear phase custom filter.
You need to define the magnitude of the response, then use this function to set the
phase appropriately for the type of filter. i.e. low pass, band pass, etc.
3. InverseFFT() is not an FFT routine, it shows how to properly fill the input array
for an Inverse FFT, set the Nyquist bin, and shows where the FIR coefficients
will be located after the Inverse FFT has been completed.
4. AdjustDelay() shows how to modify the delay of any FIR response by a fractional amount.
There is a significant amount of code you need to supply in order to use this code,
including an FFT routine and a function that defines the magnitude of your custom response.
This code assumes your FFT routine places DC in bin 0, Nyquist in bin N/2, and has the negative
frequency bins folded over around the Nyquist bin. In other words, bin 1 folds
over to bin N-1, etc. (Some FFT routines have DC in the center bin and fold around that.)
We use these 2 definitions in the code extensively.
#define iNUMSAMPLES 2048 Used for loop control
#define dNUMSAMPLES 2048.0 Used for math
These need to be a power of 2 and should be at least 1024
It does not improve this method to make them huge, but they need to be
big enough to satisfy your frequency resolution requirements. For example,
if iNUMSAMPLES = 2048, the resolution for setting the corner frequencies is Pi/1024.
The basic procedure to follow is:
1. Generate iNUMSAMPLES/2 positive frequency domain samples, either from a polynomial,
or from your custom definition.
2. Use these samples to create an input array of length iNUMSAMPLES for the Inverse FFT.
3. Perform the Inverse FFT and fill an FIR Coefficients array with the desired values.
4. Apply a window as needed. We suggest a Kaiser window for linear phase filters and a Half
cosine window for responses derived from a polynomial.
5. Adjust the group delay if needed. (This isn't typically needed.)
*/
// This function evaluates an analog filter polynomial, such as a Butterworth.
// It requires the arrays D[] and N[] be filled with the 2nd order polynomial coefficients.
// Its evaluates the polynomial H(s) at iNUMSAMPLES/2 positive frequencies and fills the complex array HofS[].
// Since the polynomial defines the phase, the CalcWinPhase() function is not used with this.
void SimulateAnalogProto(void)
{
int j, n, NPoints, NSections;
double Omega0, Q, Gain;
ComplexD s, s_sq, Numerator, Denominator, H;
NPoints = iNUMSAMPLES/2;
NSections = (Filt.NumPoles + 1) / 2; // Number of 2nd order sections.
Omega0 = 1.0 / tan(Filt.OmegaC * M_PI_2);
// Calc the Gain compensation. This is 1 for most polys.
Gain = 1.0;
for(j=0; j 1.0
Q *= Filt.OmegaC / Filt.BW;
for(j=0; j 0.5)Q = 1.0 / tan(Filt.OmegaC * M_PI_2);
Q *= Filt.OmegaC / Filt.BW;
for(j=0; j=0; j--)
{
Arg = RadPerSample * (double)(j-iNUMSAMPLES/2);
HofS[j] = WinMag[j] * ComplexD( cos(Arg), sin(Arg) );
}
break;
case NOTCH:
RightTap = ( Filt.OmegaC + Filt.BW/2.0 ) * dNUMSAMPLES / 2.0;
LeftTap = ( Filt.OmegaC - Filt.BW/2.0 ) * dNUMSAMPLES / 2.0;
CenterTap = (RightTap - LeftTap)/2.0 + LeftTap;
CenterJ = NearestInt(CenterTap); // NearestInt() rounds the float val to the nearest int val.
for(j=0; j<=CenterJ; j++)
{
Arg = RadPerSample * (double)j;
HofS[j] = WinMag[j] * ComplexD( cos(Arg), sin(Arg) );
}
for(j=iNUMSAMPLES/2; j>=CenterJ; j--)
{
Arg = RadPerSample * (double)(j-iNUMSAMPLES/2);
HofS[j] = WinMag[j] * ComplexD( cos(Arg), sin(Arg) );
}
break;
}
}
//---------------------------------------------------------------------------
// This function does this:
// 1. Fills a complex array to be used in an Inverse FFT from the HofS array filled above.
// 2. Sets the Nyquist bin appropriately, depending on the type of filter to be generated.
// 3. Calls an Inverse FFT routine (not provided here).
// 4. Gets the FIR coefficients from the output of the FFT.
// This assumes HofS[] contains iNUMSAMPLES/2 positive complex frequency samples.
void InverseFFT(void)
{
int t, f, StartT; // f and t are freq and time indexes.
// Fill the FFT's input array with the positive freq samples.
// Then fill the negative frequency bins with conjugate values.
for(f=0; f= NumTaps
dFFTSize = FFTSize; // A double var used in calculations.
FFTInput = new(std::nothrow) ComplexD[FFTSize];
if(FFTOutput == NULL)
{
ShowMessage("Unable to Allocate memory in TopForm::AdjustDelay");
return;
}
for(j=0; j