// 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); }