WSL/SLF GitLab Repository

Commit b0af40cb authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The Butterworth filter has been replaced by a critically damped filter...

The Butterworth filter has been replaced by a critically damped filter implemented as IIRLowPass that properly corrects for the phase (ie behaves as zero-phase filter).
parent 236578fa
......@@ -11,7 +11,7 @@ SET(meteoFilters_sources
meteoFilters/FilterRate.cc
meteoFilters/FilterUnheatedPSUM.cc
meteoFilters/ProcAggregate.cc
meteoFilters/ProcButterworth.cc
meteoFilters/ProcIIRLowPass.cc
meteoFilters/ProcUnventilatedT.cc
meteoFilters/ProcShade.cc
meteoFilters/ProcPSUMDistribute.cc
......
......@@ -15,7 +15,7 @@
You should have received a copy of the GNU Lesser General Public License
along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
*/
#include <meteoio/meteoFilters/ProcButterworth.h>
#include <meteoio/meteoFilters/ProcIIRLowPass.h>
#include <meteoio/MathOptim.h>
#include <cmath>
......@@ -23,32 +23,32 @@ using namespace std;
namespace mio {
const double ProcButterworth::n = 1.; //one pass
const double ProcIIRLowPass::n = 1.; //one pass
//Butterworth
//const double ProcButterworth::g = 1.;
//const double ProcButterworth::p = sqrt(2.);
//const double ProcButterworth::c = 1. / pow( pow(2, 1./n) - 1., 1./4. ); //3dB cutoff correction
//const double ProcIIRLowPass::g = 1.;
//const double ProcIIRLowPass::p = sqrt(2.);
//const double ProcIIRLowPass::c = 1. / pow( pow(2, 1./n) - 1., 1./4. ); //3dB cutoff correction
//critically damped filter
const double ProcButterworth::g = 1.;
const double ProcButterworth::p = 2.;
const double ProcButterworth::c = 1. / sqrt( pow(2., 1./(2.*n)) - 1. ); //3dB cutoff correction
const double ProcIIRLowPass::g = 1.;
const double ProcIIRLowPass::p = 2.;
const double ProcIIRLowPass::c = 1. / sqrt( pow(2., 1./(2.*n)) - 1. ); //3dB cutoff correction
//Bessel filter
//const double ProcButterworth::g = 3.;
//const double ProcButterworth::p = 3.;
//const double ProcButterworth::c = 1./ sqrt( sqrt(pow(2., 1/n) - 3./4.) - 0.5 ) / sqrt(3.); //3dB cutoff correction
//const double ProcIIRLowPass::g = 3.;
//const double ProcIIRLowPass::p = 3.;
//const double ProcIIRLowPass::c = 1./ sqrt( sqrt(pow(2., 1/n) - 3./4.) - 0.5 ) / sqrt(3.); //3dB cutoff correction
ProcButterworth::ProcButterworth(const std::vector<std::string>& vec_args, const std::string& name)
: ProcessingBlock(name), cutoff(0.)
ProcIIRLowPass::ProcIIRLowPass(const std::vector<std::string>& vec_args, const std::string& name)
: ProcessingBlock(name), cutoff(0.), bidirectional(true)
{
parse_args(vec_args);
properties.points_before = 2;
properties.stage = ProcessingProperties::first;
}
void ProcButterworth::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
void ProcIIRLowPass::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec)
{
ovec = ivec;
......@@ -151,7 +151,7 @@ void ProcButterworth::process(const unsigned int& param, const std::vector<Meteo
//this computes the filter coefficients for a low pass filter.
//the filter parameters are computed based on the filter type (this gives the polynomial coefficients and the cutoff correction) and the number of passes.
void ProcButterworth::computeCoefficients(const double& f_s, const double& f_0, double A[3], double B[3]) const
void ProcIIRLowPass::computeCoefficients(const double& f_s, const double& f_0, double A[3], double B[3]) const
{
//using the filter polynomials, the number of passes and the cutoff correction, compute the filter coefficients
const double f_star = c * f_0 / f_s; //corrected cutoff frequency
......@@ -168,16 +168,25 @@ void ProcButterworth::computeCoefficients(const double& f_s, const double& f_0,
B[2] = 1. - (A[0] + A[1] + A[2] + B[1]);
}
void ProcButterworth::parse_args(std::vector<std::string> vec_args)
void ProcIIRLowPass::parse_args(std::vector<std::string> vec_args)
{
vector<double> filter_args;
convert_args(1, 1, vec_args, filter_args);
if (filter_args.size() != 1)
throw InvalidArgumentException("Wrong number of arguments for filter " + getName(), AT);
cutoff = filter_args[0];
bidirectional = true;
const size_t nrArgs = vec_args.size();
bool period_read = false;
for (size_t ii=0; ii<nrArgs; ii++) {
if (IOUtils::isNumeric(vec_args[ii])) {
if (period_read==true)
throw InvalidArgumentException("Cutoff period has been provided more than once", AT);
if (!IOUtils::convertString(cutoff, vec_args[ii]))
throw InvalidArgumentException("Could not parse cutoff period '"+vec_args[ii]+"'", AT);
period_read = true;
} else if (vec_args[ii]=="single_pass") {
bidirectional = false;
} else
throw InvalidArgumentException("Invalid argument \""+vec_args[ii]+"\" for filter \""+getName()+"\"", AT);
}
if (!period_read)
throw InvalidArgumentException("Please provide the cutoff period for filter " + getName(), AT);
}
}
......@@ -15,8 +15,8 @@
You should have received a copy of the GNU Lesser General Public License
along with MeteoIO. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PROCBUTTERWORTH_H
#define PROCBUTTERWORTH_H
#ifndef PROCIIRLOWPASS_H
#define PROCIIRLOWPASS_H
#include <meteoio/meteoFilters/FilterBlock.h>
#include <vector>
......@@ -25,25 +25,26 @@
namespace mio {
/**
* @class ProcButterworth
* @class ProcIIRLowPass
* @ingroup processing
* @brief Simple 2 poles Butterworth low pass filter.
* The cutoff <b>period</b> (defined as the frequency at a -3dB gain) is given in seconds as argument. The phase is removed by running
* the filter twice, first backward and then forward (this also squares the amplitude response, ie doubles the order of the filter, see
* http://www.dspguide.com/ch19/4.htm or http://unicorn.us.com/trading/allpolefilters.html)
* @brief Two poles Critically Damped low pass filter.
* The cutoff <b>period</b> (defined as the frequency at a -3dB gain) is given in seconds as argument. The phase is removed by
* bidirectional filtering, ie. running the filter twice, first backward and then forward (this also squares the amplitude response, see
* http://www.dspguide.com/ch19/4.htm or http://unicorn.us.com/trading/allpolefilters.html) with a mechanism to disable filtering for
* the points that generate overshooting (basically, it runs forward/backward as well as backward/forward and removes that points that
* are too different between the two since these are indicators of overshooting).
*
* The original paper is <i>On the Theory of Filters Amplifiers</i>, S. Butterworth, Experimental wireless & the wireless engineer,
* <b>7</b>, pp 536-541, 1930.
* @code
* VW::filter1 = BUTTERWORTH
* VW::filter1 = Low_Pass
* VW::arg1 = 10800 ;3 hours
* @endcode
*
* It is possible to disable the bidirectional filtering by adding the "single_pass" argument.
*/
class ProcButterworth : public ProcessingBlock {
class ProcIIRLowPass : public ProcessingBlock {
public:
ProcButterworth(const std::vector<std::string>& vec_args, const std::string& name);
ProcIIRLowPass(const std::vector<std::string>& vec_args, const std::string& name);
virtual void process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec);
......
......@@ -32,7 +32,7 @@
#include <meteoio/meteoFilters/FilterTukey.h>
#include <meteoio/meteoFilters/FilterMAD.h>
#include <meteoio/meteoFilters/ProcAggregate.h>
#include <meteoio/meteoFilters/ProcButterworth.h>
#include <meteoio/meteoFilters/ProcIIRLowPass.h>
#include <meteoio/meteoFilters/ProcUndercatch_WMO.h>
#include <meteoio/meteoFilters/ProcUndercatch_Forland.h>
#include <meteoio/meteoFilters/ProcUndercatch_Hamon.h>
......@@ -107,7 +107,7 @@ namespace mio {
* As well as more specific data transformations:
* - EXP_SMOOTHING: exponential smoothing of data, see ProcExpSmoothing
* - WMA_SMOOTHING: weighted moving average smoothing of data, see ProcWMASmoothing
* - BUTTERWORTH: low pass butterworth filter, see ProcButterworth
* - LOW_PASS: low pass critically damped filter, see ProcIIRLowPass
* - AGGREGATE: various data aggregation algorithms, see ProcAggregate
* - UNDERCATCH_WMO: WMO rain gauge correction for undercatch, using various correction models, see ProcUndercatch_WMO
* - UNDERCATCH_FORLAND: Forland1996 rain gauge correction for solid and liquid undercatch, using various correction models, see ProcUndercatch_Forland
......@@ -165,8 +165,8 @@ ProcessingBlock* BlockFactory::getBlock(const std::string& blockname, const std:
return new ProcExpSmoothing(vec_args, blockname);
} else if (blockname == "WMA_SMOOTHING"){
return new ProcWMASmoothing(vec_args, blockname);
} else if (blockname == "BUTTERWORTH"){
return new ProcButterworth(vec_args, blockname);
} else if (blockname == "LOW_PASS"){
return new ProcIIRLowPass(vec_args, blockname);
} else if (blockname == "AGGREGATE"){
return new ProcAggregate(vec_args, blockname);
} else if (blockname == "UNDERCATCH_WMO"){
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment