/*
By Iowa Hills Software, LLC IowaHills.com
If you find a problem with this code, please leave a note at:
http://www.iowahills.com/feedbackcomments.html
November 15, 2014
This is the FFT code contained in the IIR FIR Source Code Kit
*/
#include "FFTCode.h"
#include // For the ShowMessage function.
#include // For the new operator.
#include
//---------------------------------------------------------------------------
// This calculates the required FFT size for a given number of points.
int RequiredFFTSize(int NumPts)
{
int N = MINIMUM_FFT_SIZE;
while(N < NumPts && N < MAXIMUM_FFT_SIZE)
{
N *= 2;
}
return N;
}
//---------------------------------------------------------------------------
// This verifies that the x = 2^M and returns M
// Must be >= 8 for the Twiddle calculations
int IsValidFFTSize(int x)
{
if(x < MINIMUM_FFT_SIZE || x > MAXIMUM_FFT_SIZE || (x & (x - 1)) != 0)return(0); // x & (x - 1) ensures a power of 2
return ( (int)( log((double)x) / M_LN2 + 0.5 ) ); // return M where N = 2^M
}
//---------------------------------------------------------------------------
// Fast Fourier Transform
// This code puts DC in bin 0 and scales the output of a forward transform by 1/N.
// InputR and InputI are the real and imaginary input arrays of length N.
// The output values are returned in the Input arrays.
// TTransFormType is either FORWARD or INVERSE (defined in the header file)
// 256 pts in 50 us
void FFT(double *InputR, double *InputI, int N, TTransFormType Type)
{
int j, LogTwoOfN, *RevBits;
double *BufferR, *BufferI, *TwiddleR, *TwiddleI;
double OneOverN;
// Verify the FFT size and type.
LogTwoOfN = IsValidFFTSize(N);
if(LogTwoOfN == 0 || (Type != FORWARD && Type != INVERSE) )
{
ShowMessage("Invalid FFT type or size.");
return;
}
// Memory allocation for all the arrays.
BufferR = new(std::nothrow) double[N];
BufferI = new(std::nothrow) double[N];
TwiddleR = new(std::nothrow) double[N/2];
TwiddleI = new(std::nothrow) double[N/2];
RevBits = new(std::nothrow) int[N];
if(BufferR == NULL || BufferI == NULL ||
TwiddleR == NULL || TwiddleI == NULL || RevBits == NULL)
{
ShowMessage("FFT Memory Allocation Error");
return;
}
ReArrangeInput(InputR, InputI, BufferR, BufferI, RevBits, N);
FillTwiddleArray(TwiddleR, TwiddleI, N, Type);
Transform(InputR, InputI, BufferR, BufferI, TwiddleR, TwiddleI, N);
// The ReArrangeInput function swapped Input[] and Buffer[]. Then Transform()
// swapped them again, LogTwoOfN times. Ultimately, these swaps must be done
// an even number of times, or the pointer to Buffer gets returned.
// So we must do one more swap here, for N = 16, 64, 256, 1024, ...
OneOverN = 1.0;
if(Type == FORWARD) OneOverN = 1.0 / (double)N;
if(LogTwoOfN % 2 == 1)
{
for(j=0; j= 1)
{
for(k=0; k 0) // Loops Log2(N) times.
{
// Swap pointers, instead doing this: for(j=0; j 0.0)Mag = sqrt(Mag);
else Mag = 1.0E-12;
return(Mag);
}
//---------------------------------------------------------------------------
/*
These are the window definitions. These windows can be used for either
FIR filter design or with an FFT for spectral analysis.
For definitions, see this article: http://en.wikipedia.org/wiki/Window_function
This function has 6 inputs
Data is the array, of length N, containing the data to to be windowed.
This data is either a FIR filter sinc pulse, or the data to be analyzed by an fft.
WindowType is an enum defined in the header file.
e.g. wtKAISER, wtSINC, wtHANNING, wtHAMMING, wtBLACKMAN, ...
Alpha sets the width of the flat top.
Windows such as the Tukey and Trapezoid are defined to have a variably wide flat top.
As can be seen by its definition, the Tukey is just a Hanning window with a flat top.
Alpha can be used to give any of these windows a partial flat top, except the Flattop and Kaiser.
Alpha = 0 gives the original window. (i.e. no flat top)
To generate a Tukey window, use a Hanning with 0 < Alpha < 1
To generate a Bartlett window (triangular), use a Trapezoid window with Alpha = 0.
Alpha = 1 generates a rectangular window in all cases. (except the Flattop and Kaiser)
Beta is used with the Kaiser, Sinc, and Sine windows only.
These three windows are used primarily for FIR filter design. Then
Beta controls the filter's transition bandwidth and the sidelobe levels.
All other windows ignore Beta.
UnityGain controls whether the gain of these windows is set to unity.
Only the Flattop window has unity gain by design. The Hanning window, for example, has a gain
of 1/2. UnityGain = true sets the gain to 1, which preserves the signal's energy level
when these windows are used for spectral analysis.
Don't use this with FIR filter design however. Since most of the enegy in an FIR sinc pulse
is in the middle of the window, the window needs a peak amplitude of one, not unity gain.
Setting UnityGain = true will simply cause the resulting FIR filter to have excess gain.
If using these windows for FIR filters, start with the Kaiser, Sinc, or Sine windows and
adjust Beta for the desired transition BW and sidelobe levels (set Alpha = 0).
While the FlatTop is an excellent window for spectral analysis, don't use it for FIR filter design.
It has a peak amplitude of ~ 4.7 which causes the resulting FIR filter to have about this much gain.
It works poorly for FIR filters even if you adjust its peak amplitude.
The Trapezoid also works poorly for FIR filter design.
If using these windows with an fft for spectral analysis, start with the Hanning, Gauss, or Flattop.
When choosing a window for spectral analysis, you must trade off between resolution and amplitude
accuracy. The Hanning has the best resolution while the Flatop has the best amplitude accuracy.
The Gauss is midway between these two for both accuracy and resolution. These three were
the only windows available in the HP 89410A Vector Signal Analyzer. Which is to say, these three
are the probably the best windows for general purpose signal analysis.
*/
void WindowData(double *Data, int N, TWindowType WindowType, double Alpha, double Beta, bool UnityGain)
{
if(WindowType == wtNONE) return;
int j, M, TopWidth;
double dM, *WinCoeff;
if(WindowType == wtKAISER || WindowType == wtFLATTOP )Alpha = 0.0;
if(Alpha < 0.0)Alpha = 0.0;
if(Alpha > 1.0)Alpha = 1.0;
if(Beta < 0.0)Beta = 0.0;
if(Beta > 10.0)Beta = 10.0;
WinCoeff = new(std::nothrow) double[N+2];
if(WinCoeff == NULL)
{
ShowMessage("Failed to allocate memory in WindowData() ");
return;
}
TopWidth = (int)( Alpha * (double)N );
if(TopWidth%2 != 0)TopWidth++;
if(TopWidth > N)TopWidth = N;
M = N - TopWidth;
dM = M + 1;
// Calculate the window for N/2 points, then fold the window over (at the bottom).
// TopWidth points will be set to 1.
if(WindowType == wtKAISER)
{
double Arg;
for(j=0; j 0. Cannot be applied to a Kaiser or Flat Top.
if(WindowType != wtKAISER && WindowType != wtFLATTOP)
{
for(j=M/2; j -1.0E-5 && x < 1.0E-5)return(1.0);
return(sin(x)/x);
}
//---------------------------------------------------------------------------
/*
The header file contents.
#ifndef FFTCodeH
#define FFTCodeH
#define M_2PI 6.28318530717958647692
#define MAXIMUM_FFT_SIZE 1048576
#define MINIMUM_FFT_SIZE 8
//---------------------------------------------------------------------------
// Must retain the order on the 1st line for legacy FIR code.
enum TWindowType {wtFIRSTWINDOW, wtNONE, wtKAISER, wtSINC, wtHANNING,
wtHAMMING, wtBLACKMAN, wtFLATTOP, wtBLACKMAN_HARRIS,
wtBLACKMAN_NUTTALL, wtNUTTALL, wtKAISER_BESSEL, wtTRAPEZOID,
wtGAUSS, wtSINE, wtTEST };
enum TTransFormType {FORWARD, INVERSE};
int RequiredFFTSize(int NumPts);
int IsValidFFTSize(int x);
void FFT(double *InputR, double *InputI, int N, TTransFormType Type);
void ReArrangeInput(double *InputR, double *InputI, double *BufferR, double *BufferI, int *RevBits, int N);
void FillTwiddleArray(double *TwiddleR, double *TwiddleI, int N, TTransFormType Type);
void Transform(double *InputR, double *InputI, double *BufferR, double *BufferI, double *TwiddleR, double *TwiddleI, int N);
void DFT(double *InputR, double *InputI, int N, int Type);
void RealSigDFT(double *Samples, double *OutputR, double *OutputI, int N);
double SingleFreqDFT(double *Samples, int N, double Omega);
double Goertzel(double *Samples, int N, double Omega);
void WindowData(double *Data, int N, TWindowType WindowType, double Alpha, double Beta, bool UnityGain);
double Bessel(double x);
double Sinc(double x);
#endif
*/