// October 18, 2014
// From: http://www.iowahills.com/A7ExampleCodePage.html
// If you find a problem with this code, please leave us a note on:
// http://www.iowahills.com/feedbackcomments.html
/*
The intent here is to show techniques for calculating the frequency
response of an FIR filter, not to provide ready to use code.
There are 3 functions given here. The first is a DFT, which calculates the frequency
response of an FIR filter for 0.0 < Omega < 1.0
The second function is a DFT which calculates the FIR response at one frequency.
The third function is the Goertzel algorithm, which also calculates the FIR response
at a single frequency, but unlike a DFT, it provides no phase information.
See these links for more information on Goertzel.
http://en.wikipedia.org/wiki/Goertzel_algorithm
http://www.embedded.com/design/real-world-applications/4401754/Single-tone-detection-with-the-Goertzel-algorithm
http://netwerkt.wordpress.com/2011/08/25/goertzel-filter/
*/
// FIRFreqResponse is a DFT used to calculate the frequency response of an FIR filter.
// The only reason to use this instead of an FFT is that we can evaluate the filter at
// the frequencies of our choosing, as opposed to the FFT's predetermined bin frequencies.
// It calculates the response at iNUM_PLOT_PTS frequencies. Normally, 0 <= Omega < 1
// It fills FFTOutput[] which gets used in the plotting routines.
// This function gets called hundreds of times. To save time, the twiddle factors are stored.
// They only need to be calculated on the first call, or when MaxPlotOmega changes.
// MaxPlotOmega is usually equal to 1.0, but can be any value.
// ComplexD is a complex double.
// TwiddleReal and TwiddleImag are global arrays of doubles.
// FFTOutput is a complex array used to store the response.
#define iNUM_PLOT_PTS 1024 // These can be any convenient value.
#define dNUM_PLOT_PTS 1024.0
void FIRFreqResponse(double *Samples, int N) // Samples[] = FirCoeff[] N = NumTaps
{
int j, k;
double Arg, Temp, zReal, zImag;
double SumReal, SumImag; // Real and Imag parts of the Sum.
static double PrevMaxPlotOmega = -1000.0;
// Use precalculated twiddle factors
if(PrevMaxPlotOmega != MaxPlotOmega)
{
PrevMaxPlotOmega = MaxPlotOmega;
for(j=0; j 0.0)Mag = sqrt(Mag);
else Mag = 1.0E-12; // To prevent a problem in dB()
return(Mag);
}