From 6c9dbb79c42c2741ef057b4f6ee17ea7e8b468c6 Mon Sep 17 00:00:00 2001 From: f4exb Date: Mon, 4 Sep 2017 01:44:36 +0200 Subject: [PATCH] Added a windowed FIR filter calculator analogous to Scipy's firwin (this one from Iowa Hills Software LLC code) --- CMakeLists.txt | 2 + sdrbase/dsp/wfir.cpp | 377 +++++++++++++++++++++++++++++++++++++++++++ sdrbase/dsp/wfir.h | 93 +++++++++++ sdrbase/sdrbase.pro | 2 + 4 files changed, 474 insertions(+) create mode 100644 sdrbase/dsp/wfir.cpp create mode 100644 sdrbase/dsp/wfir.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e6133824..b49db3af0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -136,6 +136,7 @@ set(sdrbase_SOURCES sdrbase/dsp/spectrumvis.cpp sdrbase/dsp/threadedbasebandsamplesink.cpp sdrbase/dsp/threadedbasebandsamplesource.cpp + sdrbase/dsp/wfir.cpp sdrbase/gui/aboutdialog.cpp sdrbase/gui/addpresetdialog.cpp @@ -257,6 +258,7 @@ set(sdrbase_HEADERS sdrbase/dsp/spectrumvis.h sdrbase/dsp/threadedbasebandsamplesink.h sdrbase/dsp/threadedbasebandsamplesource.h + sdrbase/dsp/wfir.h sdrbase/gui/aboutdialog.h sdrbase/gui/addpresetdialog.h diff --git a/sdrbase/dsp/wfir.cpp b/sdrbase/dsp/wfir.cpp new file mode 100644 index 000000000..14414ebe6 --- /dev/null +++ b/sdrbase/dsp/wfir.cpp @@ -0,0 +1,377 @@ +/* + July 15, 2015 + Iowa Hills Software LLC + http://www.iowahills.com + */ + +#include +#include +#include + +#include "wfir.h" + +#define M_2PI 6.28318530717958647692 // M_PI should be in the math.h file + +// This first calculates the impulse response for a rectangular window. +// It then applies the windowing function of choice to the impulse response. +void WFIR::BasicFIR(double *FirCoeff, int NumTaps, TPassTypeName PassType, + double OmegaC, double BW, TWindowType WindowType, double WinBeta) +{ + int j; + double Arg, OmegaLow, OmegaHigh; + + switch (PassType) + { + case LPF: + for (j = 0; j < NumTaps; j++) + { + Arg = (double) j - (double) (NumTaps - 1) / 2.0; + FirCoeff[j] = OmegaC * Sinc(OmegaC * Arg * M_PI); + } + break; + + case HPF: + if (NumTaps % 2 == 1) // Odd tap counts + { + for (j = 0; j < NumTaps; j++) + { + Arg = (double) j - (double) (NumTaps - 1) / 2.0; + FirCoeff[j] = Sinc(Arg * M_PI) + - OmegaC * Sinc(OmegaC * Arg * M_PI); + } + } + + else // Even tap counts + { + for (j = 0; j < NumTaps; j++) + { + Arg = (double) j - (double) (NumTaps - 1) / 2.0; + if (Arg == 0.0) + FirCoeff[j] = 0.0; + else + FirCoeff[j] = cos(OmegaC * Arg * M_PI) / M_PI / Arg + + cos(Arg * M_PI); + } + } + break; + + case BPF: + OmegaLow = OmegaC - BW / 2.0; + OmegaHigh = OmegaC + BW / 2.0; + for (j = 0; j < NumTaps; j++) + { + Arg = (double) j - (double) (NumTaps - 1) / 2.0; + if (Arg == 0.0) + FirCoeff[j] = 0.0; + else + FirCoeff[j] = (cos(OmegaLow * Arg * M_PI) + - cos(OmegaHigh * Arg * M_PI)) / M_PI / Arg; + } + break; + + case NOTCH: // If NumTaps is even for Notch filters, the response at Pi is attenuated. + OmegaLow = OmegaC - BW / 2.0; + OmegaHigh = OmegaC + BW / 2.0; + for (j = 0; j < NumTaps; j++) + { + Arg = (double) j - (double) (NumTaps - 1) / 2.0; + FirCoeff[j] = Sinc(Arg * M_PI) + - OmegaHigh * Sinc(OmegaHigh * Arg * M_PI) + - OmegaLow * Sinc(OmegaLow * Arg * M_PI); + } + break; + } + + // WindowData can be used to window data before an FFT. When used for FIR filters we set + // Alpha = 0.0 to prevent a flat top on the window and + // set UnityGain = false to prevent the window gain from getting set to unity. + WindowData(FirCoeff, NumTaps, WindowType, 0.0, WinBeta, false); + +} + +//--------------------------------------------------------------------------- + +// This gets used with the Kaiser window. +double WFIR::Bessel(double x) +{ + double Sum = 0.0, XtoIpower; + int i, j, Factorial; + for (i = 1; i < 10; i++) + { + XtoIpower = pow(x / 2.0, (double) i); + Factorial = 1; + for (j = 1; j <= i; j++) + Factorial *= j; + Sum += pow(XtoIpower / (double) Factorial, 2.0); + } + return (1.0 + Sum); +} + +//----------------------------------------------------------------------------- + +// This gets used with the Sinc window and various places in the BasicFIR function. +double WFIR::Sinc(double x) +{ + if (x > -1.0E-5 && x < 1.0E-5) + return (1.0); + return (sin(x) / x); +} + +//--------------------------------------------------------------------------- + +// These are the various windows definitions. These windows can be used for either +// FIR filter design or with an FFT for spectral analysis. +// Sourced verbatim from: ~MyDocs\Code\Common\FFTFunctions.cpp +// 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, which is at the bottom of this 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 primarily used for FIR filter design, not spectral analysis. +// In FIR filter design, Beta controls the filter's transition bandwidth and the sidelobe levels. +// The code ignores Beta except in the Kaiser, Sinc, and Sine window cases. + +// 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 will set the gain of all these windows to 1. +// Then, when the window is applied to a signal, the signal's energy content is preserved. +// 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, +// unless you have specific windowing requirements, use one of these 3 for general purpose signal analysis. +// Set UnityGain = true when using any of these windows for spectral analysis to preserve the signal's enegy level. + +void WFIR::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 == 0) + { + std::cerr + << "Failed to allocate memory in FFTFunctions::WindowFFTData() " + << std::endl; + 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 < M; j++) + { + Arg = Beta * sqrt(1.0 - pow(((double) (2 * j + 2) - dM) / dM, 2.0)); + WinCoeff[j] = Bessel(Arg) / Bessel(Beta); + } + } + + else if (WindowType == wtSINC) // Lanczos + { + for (j = 0; j < M; j++) + WinCoeff[j] = Sinc((double) (2 * j + 1 - M) / dM * M_PI); + for (j = 0; j < M; j++) + WinCoeff[j] = pow(WinCoeff[j], Beta); + } + + else if (WindowType == wtSINE) // Hanning if Beta = 2 + { + for (j = 0; j < M / 2; j++) + WinCoeff[j] = sin((double) (j + 1) * M_PI / dM); + for (j = 0; j < M / 2; j++) + WinCoeff[j] = pow(WinCoeff[j], Beta); + } + + else if (WindowType == wtHANNING) + { + for (j = 0; j < M / 2; j++) + WinCoeff[j] = 0.5 - 0.5 * cos((double) (j + 1) * M_2PI / dM); + } + + else if (WindowType == wtHAMMING) + { + for (j = 0; j < M / 2; j++) + WinCoeff[j] = 0.54 - 0.46 * cos((double) (j + 1) * M_2PI / dM); + } + + else if (WindowType == wtBLACKMAN) + { + for (j = 0; j < M / 2; j++) + { + WinCoeff[j] = 0.42 - 0.50 * cos((double) (j + 1) * M_2PI / dM) + + 0.08 * cos((double) (j + 1) * M_2PI * 2.0 / dM); + } + } + + // See: http://www.bth.se/fou/forskinfo.nsf/0/130c0940c5e7ffcdc1256f7f0065ac60/$file/ICOTA_2004_ttr_icl_mdh.pdf + else if (WindowType == wtFLATTOP) + { + for (j = 0; j <= M / 2; j++) + { + WinCoeff[j] = 1.0 + - 1.93293488969227 * cos((double) (j + 1) * M_2PI / dM) + + 1.28349769674027 + * cos((double) (j + 1) * M_2PI * 2.0 / dM) + - 0.38130801681619 + * cos((double) (j + 1) * M_2PI * 3.0 / dM) + + 0.02929730258511 + * cos((double) (j + 1) * M_2PI * 4.0 / dM); + } + } + + else if (WindowType == wtBLACKMAN_HARRIS) + { + for (j = 0; j < M / 2; j++) + { + WinCoeff[j] = 0.35875 - 0.48829 * cos((double) (j + 1) * M_2PI / dM) + + 0.14128 * cos((double) (j + 1) * M_2PI * 2.0 / dM) + - 0.01168 * cos((double) (j + 1) * M_2PI * 3.0 / dM); + } + } + + else if (WindowType == wtBLACKMAN_NUTTALL) + { + for (j = 0; j < M / 2; j++) + { + WinCoeff[j] = 0.3535819 + - 0.4891775 * cos((double) (j + 1) * M_2PI / dM) + + 0.1365995 * cos((double) (j + 1) * M_2PI * 2.0 / dM) + - 0.0106411 * cos((double) (j + 1) * M_2PI * 3.0 / dM); + } + } + + else if (WindowType == wtNUTTALL) + { + for (j = 0; j < M / 2; j++) + { + WinCoeff[j] = 0.355768 + - 0.487396 * cos((double) (j + 1) * M_2PI / dM) + + 0.144232 * cos((double) (j + 1) * M_2PI * 2.0 / dM) + - 0.012604 * cos((double) (j + 1) * M_2PI * 3.0 / dM); + } + } + + else if (WindowType == wtKAISER_BESSEL) + { + for (j = 0; j <= M / 2; j++) + { + WinCoeff[j] = 0.402 - 0.498 * cos(M_2PI * (double) (j + 1) / dM) + + 0.098 * cos(2.0 * M_2PI * (double) (j + 1) / dM) + + 0.001 * cos(3.0 * M_2PI * (double) (j + 1) / dM); + } + } + + else if (WindowType == wtTRAPEZOID) // Rectangle for Alpha = 1 Triangle for Alpha = 0 + { + int K = M / 2; + if (M % 2) + K++; + for (j = 0; j < K; j++) + WinCoeff[j] = (double) (j + 1) / (double) K; + } + + // This definition is from http://en.wikipedia.org/wiki/Window_function (Gauss Generalized normal window) + // We set their p = 2, and use Alpha in the numerator, instead of Sigma in the denominator, as most others do. + // Alpha = 2.718 puts the Gauss window response midway between the Hanning and the Flattop (basically what we want). + // It also gives the same BW as the Gauss window used in the HP 89410A Vector Signal Analyzer. + // Alpha = 1.8 puts it quite close to the Hanning. + else if (WindowType == wtGAUSS) + { + for (j = 0; j < M / 2; j++) + { + WinCoeff[j] = ((double) (j + 1) - dM / 2.0) / (dM / 2.0) * 2.7183; + WinCoeff[j] *= WinCoeff[j]; + WinCoeff[j] = exp(-WinCoeff[j]); + } + } + + else // Error. + { + std::cerr << "Incorrect window type in WindowFFTData" << std::endl; + delete[] WinCoeff; + return; + } + + // Fold the coefficients over. + for (j = 0; j < M / 2; j++) + WinCoeff[N - j - 1] = WinCoeff[j]; + + // This is the flat top if Alpha > 0. Cannot be applied to a Kaiser or Flat Top. + if (WindowType != wtKAISER && WindowType != wtFLATTOP) + { + for (j = M / 2; j < N - M / 2; j++) + WinCoeff[j] = 1.0; + } + + // This will set the gain of the window to 1. Only the Flattop window has unity gain by design. + if (UnityGain) + { + double Sum = 0.0; + for (j = 0; j < N; j++) + Sum += WinCoeff[j]; + Sum /= (double) N; + if (Sum != 0.0) + for (j = 0; j < N; j++) + WinCoeff[j] /= Sum; + } + + // Apply the window to the data. + for (j = 0; j < N; j++) + Data[j] *= WinCoeff[j]; + + delete[] WinCoeff; + +} diff --git a/sdrbase/dsp/wfir.h b/sdrbase/dsp/wfir.h new file mode 100644 index 000000000..981f80a24 --- /dev/null +++ b/sdrbase/dsp/wfir.h @@ -0,0 +1,93 @@ +/* + July 15, 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 + + Source: ~Projects\Common\BasicFIRFilterCode.cpp + + This generic FIR filter code is described in most textbooks. + e.g. Discrete Time Signal Processing, Oppenheim and Shafer + + A nice paper on this topic is: + http://dea.brunel.ac.uk/cmsp/Home_Saeed_Vaseghi/Chapter05-DigitalFilters.pdf + + This code first generates either a low pass, high pass, band pass, or notch + impulse response for a rectangular window. It then applies a window to this + impulse response. + + There are several windows available, including the Kaiser, Sinc, Hanning, + Blackman, and Hamming. Of these, the Kaiser and Sinc are probably the most useful + for FIR filters because their sidelobe levels can be controlled with the Beta parameter. + + This is a typical function call: + BasicFIR(FirCoeff, NumTaps, PassType, OmegaC, BW, wtKAISER, Beta); + BasicFIR(FirCoeff, 33, LPF, 0.2, 0.0, wtKAISER, 3.2); + 33 tap, low pass, corner frequency at 0.2, BW=0 (ignored in the low pass code), + Kaiser window, Kaiser Beta = 3.2 + + These variables should be defined similar to this: + double FirCoeff[MAXNUMTAPS]; + int NumTaps; NumTaps can be even or odd, but must be less than the FirCoeff array size. + TPassTypeName PassType; PassType is an enum defined in the header file. LPF, HPF, BPF, or NOTCH + double OmegaC 0.0 < OmegaC < 1.0 The filters corner freq, or center freq if BPF or NOTCH + double BW 0.0 < BW < 1.0 The filters band width if BPF or NOTCH + TWindowType WindowType; WindowType is an enum defined in the header to be one of these. + wtNONE, wtKAISER, wtSINC, wtHANNING, .... and others. + double Beta; 0 <= Beta <= 10.0 Beta is used with the Kaiser, Sinc, and Sine windows only. + It controls the transition BW and sidelobe level of the filters. + + + If you want to use it, Kaiser originally defined Beta as follows. + He derived its value based on the desired sidelobe level, dBAtten. + double dBAtten, Beta, Beta1=0.0, Beta2=0.0; + if(dBAtten < 21.0)dBAtten = 21.0; + if(dBAtten > 50.0)Beta1 = 0.1102 * (dBAtten - 8.7); + if(dBAtten >= 21.0 && dBAtten <= 50.0) Beta2 = 0.5842 * pow(dBAtten - 21.0, 0.4) + 0.07886 * (dBAtten - 21.0); + Beta = Beta1 + Beta2; + + */ + +#ifndef _WFIR_H_ +#define _WFIR_H_ + +class WFIR +{ +public: + enum TPassTypeName + { + LPF, HPF, BPF, NOTCH + }; + + enum TWindowType + { + wtNONE, + wtKAISER, + wtSINC, + wtHANNING, + wtHAMMING, + wtBLACKMAN, + wtFLATTOP, + wtBLACKMAN_HARRIS, + wtBLACKMAN_NUTTALL, + wtNUTTALL, + wtKAISER_BESSEL, + wtTRAPEZOID, + wtGAUSS, + wtSINE, + wtTEST + }; + + static void BasicFIR(double *FirCoeff, int NumTaps, TPassTypeName PassType, + double OmegaC, double BW, TWindowType WindowType, double WinBeta); + +private: + static void WindowData(double *Data, int N, TWindowType WindowType, + double Alpha, double Beta, bool UnityGain); + static double Bessel(double x); + static double Sinc(double x); +}; + +#endif diff --git a/sdrbase/sdrbase.pro b/sdrbase/sdrbase.pro index c546b0aaf..4054a5617 100644 --- a/sdrbase/sdrbase.pro +++ b/sdrbase/sdrbase.pro @@ -89,6 +89,7 @@ SOURCES += mainwindow.cpp\ dsp/spectrumvis.cpp\ dsp/threadedbasebandsamplesink.cpp\ dsp/threadedbasebandsamplesource.cpp\ + dsp/wfir.cpp\ gui/aboutdialog.cpp\ gui/addpresetdialog.cpp\ gui/basicchannelsettingswidget.cpp\ @@ -196,6 +197,7 @@ HEADERS += mainwindow.h\ dsp/spectrumvis.h\ dsp/threadedbasebandsamplesink.h\ dsp/threadedbasebandsamplesource.h\ + dsp/wfir.h\ gui/aboutdialog.h\ gui/addpresetdialog.h\ gui/audiodialog.h\