WSL/SLF GitLab Repository

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

Since the IIR filter implementation can do both High and Low pass, it has been...

Since the IIR filter implementation can do both High and Low pass, it has been renamed to IIR and options have been added to use it as either low or high pass filter. The documentation has been improved and stability constraints now lead to a warning message if the filter would start to be unstable.
parent b0af40cb
......@@ -11,7 +11,7 @@ SET(meteoFilters_sources
meteoFilters/FilterRate.cc
meteoFilters/FilterUnheatedPSUM.cc
meteoFilters/ProcAggregate.cc
meteoFilters/ProcIIRLowPass.cc
meteoFilters/ProcIIR.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/ProcIIRLowPass.h>
#include <meteoio/meteoFilters/ProcIIR.h>
#include <meteoio/MathOptim.h>
#include <cmath>
......@@ -23,32 +23,16 @@ using namespace std;
namespace mio {
const double ProcIIRLowPass::n = 1.; //one pass
//Butterworth
//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 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 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
ProcIIRLowPass::ProcIIRLowPass(const std::vector<std::string>& vec_args, const std::string& name)
: ProcessingBlock(name), cutoff(0.), bidirectional(true)
ProcIIR::ProcIIR(const std::vector<std::string>& vec_args, const std::string& name)
: ProcessingBlock(name), cutoff(0.), g(0.), p(0.), c(0.), bidirectional(true), low_pass(true)
{
parse_args(vec_args);
properties.points_before = 2;
properties.stage = ProcessingProperties::first;
getFilterParameters(CRITICALLY_DAMPED, low_pass, 2., g, p, c);
}
void ProcIIRLowPass::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
void ProcIIR::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec)
{
ovec = ivec;
......@@ -57,106 +41,81 @@ void ProcIIRLowPass::process(const unsigned int& param, const std::vector<MeteoD
const double days = ivec.back().date.getJulian() - ivec.front().date.getJulian();
const size_t nr_data_pts = ivec.size();
const double sampling_rate = static_cast<double>(nr_data_pts-1) / (days*24.*3600.); //in Hz
double A[3], B[3];
computeCoefficients(sampling_rate, 1./cutoff, A, B);
std::vector<double> X(3, IOUtils::nodata), Y(3, IOUtils::nodata);
std::vector<double> X(3, IOUtils::nodata), Y(3, IOUtils::nodata);
if (!bidirectional) {
//only forward filter
for (size_t ii=0; ii<ovec.size(); ++ii){
const double& raw_val = ivec[ii](param);
//propagate in X and Y
X[2] = X[1]; X[1] = X[0]; X[0] = raw_val;
Y[2] = Y[1]; Y[1] = Y[0]; Y[0] = raw_val; //Y[0] will be overwritten but in case of nodata we still propagate a value
if (X[2]==IOUtils::nodata || X[1]==IOUtils::nodata || X[0]==IOUtils::nodata) continue;
if (Y[2]==IOUtils::nodata || Y[1]==IOUtils::nodata) continue;
Y[0] = A[0]*X[0] + A[1]*X[1] + A[2]*X[2] + B[1]*Y[1] + B[2]*Y[2];
ovec[ii](param) = Y[0];
}
for (size_t ii=0; ii<ovec.size(); ++ii)
ovec[ii](param) = filterPoint(ivec[ii](param), A, B, X, Y);
} else { //bidirectional filtering
//because bidirectional filtering introduces some overshoot, we do it twice: backward/forward, then forward/backward and
//only accept the values that returned close values between the two. Large difference between these two indicate that
//there is some heavy overshoot going on and therefore we keep the original data.
std::vector<double> vecTmp(nr_data_pts), vecRes1(nr_data_pts);
//backward filter
for (size_t ii=ovec.size(); ii--> 0;){
const double& raw_val = ivec[ii](param);
//propagate in X and Y
X[2] = X[1]; X[1] = X[0]; X[0] = raw_val;
Y[2] = Y[1]; Y[1] = Y[0]; Y[0] = raw_val; //Y[0] will be overwritten but in case of nodata we still propagate a value
if (X[2]==IOUtils::nodata || X[1]==IOUtils::nodata || X[0]==IOUtils::nodata) continue;
if (Y[2]==IOUtils::nodata || Y[1]==IOUtils::nodata) continue;
Y[0] = A[0]*X[0] + A[1]*X[1] + A[2]*X[2] + B[1]*Y[1] + B[2]*Y[2];
vecTmp[ii] = Y[0];
}
//forward filter
for (size_t ii=0; ii<ovec.size(); ++ii){
const double& raw_val = vecTmp[ii];
//propagate in X and Y
X[2] = X[1]; X[1] = X[0]; X[0] = raw_val;
Y[2] = Y[1]; Y[1] = Y[0]; Y[0] = raw_val; //Y[0] will be overwritten but in case of nodata we still propagate a value
if (X[2]==IOUtils::nodata || X[1]==IOUtils::nodata || X[0]==IOUtils::nodata) continue;
if (Y[2]==IOUtils::nodata || Y[1]==IOUtils::nodata) continue;
Y[0] = A[0]*X[0] + A[1]*X[1] + A[2]*X[2] + B[1]*Y[1] + B[2]*Y[2];
vecRes1[ii] = Y[0];
}
std::vector<double> vecTmp(nr_data_pts);
for (size_t ii=ovec.size(); ii--> 0;)
vecTmp[ii] = filterPoint(ivec[ii](param), A, B, X, Y);
std::vector<double> vecRes2(nr_data_pts);
//forward filter
for (size_t ii=0; ii<ovec.size(); ++ii){
const double& raw_val = ivec[ii](param);
//propagate in X and Y
X[2] = X[1]; X[1] = X[0]; X[0] = raw_val;
Y[2] = Y[1]; Y[1] = Y[0]; Y[0] = raw_val; //Y[0] will be overwritten but in case of nodata we still propagate a value
if (X[2]==IOUtils::nodata || X[1]==IOUtils::nodata || X[0]==IOUtils::nodata) continue;
if (Y[2]==IOUtils::nodata || Y[1]==IOUtils::nodata) continue;
Y[0] = A[0]*X[0] + A[1]*X[1] + A[2]*X[2] + B[1]*Y[1] + B[2]*Y[2];
vecTmp[ii] = Y[0];
}
//backward filter
for (size_t ii=ovec.size(); ii--> 0;){
const double& raw_val = vecTmp[ii];
//propagate in X and Y
X[2] = X[1]; X[1] = X[0]; X[0] = raw_val;
Y[2] = Y[1]; Y[1] = Y[0]; Y[0] = raw_val; //Y[0] will be overwritten but in case of nodata we still propagate a value
if (X[2]==IOUtils::nodata || X[1]==IOUtils::nodata || X[0]==IOUtils::nodata) continue;
if (Y[2]==IOUtils::nodata || Y[1]==IOUtils::nodata) continue;
std::fill(X.begin(), X.end(), IOUtils::nodata);
for (size_t ii=0; ii<ovec.size(); ++ii)
ovec[ii](param) = filterPoint(vecTmp[ii], A, B, X, Y);
}
}
Y[0] = A[0]*X[0] + A[1]*X[1] + A[2]*X[2] + B[1]*Y[1] + B[2]*Y[2];
vecRes2[ii] = Y[0];
}
void ProcIIR::parse_args(std::vector<std::string> vec_args)
{
const size_t nrArgs = vec_args.size();
bool period_read = false;
for(size_t ii=0; ii<ovec.size(); ++ii) {
const double filt1 = vecRes1[ii];
const double filt2 = vecRes2[ii];
const double ref = std::max(abs(filt1), abs(filt2));
const double threshold = (ref>1e-9)? ref*1e-2 : 1e-4;
if (abs((filt1-filt2)) < threshold)
ovec[ii](param) = filt1;
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 {
const std::string arg( IOUtils::strToUpper(vec_args[ii]) );
if (arg=="SINGLE_PASS")
bidirectional = false;
else if (arg=="LP")
low_pass = true;
else if (arg=="HP")
low_pass=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);
}
double ProcIIR::filterPoint(const double& raw_val, const double A[3], const double B[3], std::vector<double> &X, std::vector<double> &Y)
{
//propagate in X and Y
X[2] = X[1]; X[1] = X[0]; X[0] = raw_val;
Y[2] = Y[1]; Y[1] = Y[0]; Y[0] = raw_val; //Y[0] will be overwritten but in case of nodata we still propagate a value
if (X[2]==IOUtils::nodata || X[1]==IOUtils::nodata || X[0]==IOUtils::nodata) return Y[0];
if (Y[2]==IOUtils::nodata || Y[1]==IOUtils::nodata) return Y[0];
Y[0] = A[0]*X[0] + A[1]*X[1] + A[2]*X[2] + B[1]*Y[1] + B[2]*Y[2];
return Y[0];
}
//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 ProcIIRLowPass::computeCoefficients(const double& f_s, const double& f_0, double A[3], double B[3]) const
//the filter parameters are computed based on the filter polynomial coefficients and the cutoff correction.
void ProcIIR::computeCoefficients(const double& fs, const double& f0, 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
const double f_star = (low_pass)? c*f0/fs : 0.5 - c*f0/fs; //corrected cutoff frequency
const double w_0 = tan(Cst::PI*f_star); //warp cutoff frequency
if ((low_pass && f_star>=0.25) || (!low_pass && f_star<=0.25)) {
std::cerr << "[W] in the '" << getName() << "' filter, the chosen cutoff frequency is incompatible with the sampling rate (unstable behavior): ";
std::cerr << "f* = " << f_star << " 1/f0 = " << 1./f0 << " 1/fs = " << 1./fs << "\n";
}
const double K1 = p * w_0;
const double K2 = g * w_0*w_0;
......@@ -166,27 +125,31 @@ void ProcIIRLowPass::computeCoefficients(const double& f_s, const double& f_0, d
B[1] = 2*A[0] * (1./K2 - 1.);
B[2] = 1. - (A[0] + A[1] + A[2] + B[1]);
if (!low_pass) { //some signs are different for High Pass
A[1] = -A[1];
B[1] = -B[1];
}
}
void ProcIIRLowPass::parse_args(std::vector<std::string> vec_args)
void ProcIIR::getFilterParameters(const IIR_Type& type, const bool& isLowPass, const double& n, double &g, double &p, double &c)
{
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);
if (type==BUTTERWORTH) {
g = 1.;
p = sqrt(2.);
c = 1. / pow( pow(2, 1./n) - 1., 1./4. ); //3dB cutoff correction
} else if (type==CRITICALLY_DAMPED) {
g = 1.;
p = 2.;
c = sqrt( pow(2., 1./(2.*n)) - 1. ); //3dB cutoff correction
} else if (type==BESSEL) {
g = 3.;
p = 3.;
c = sqrt( sqrt(pow(2., 1/n) - 3./4.) - 0.5 ) / sqrt(3.); //3dB cutoff correction
} else
throw UnknownValueException("The requested IIR filter type has not been defined", AT);
if (isLowPass) c = 1. / c;
}
}
......@@ -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 PROCIIRLOWPASS_H
#define PROCIIRLOWPASS_H
#ifndef PROCIIR_H
#define PROCIIR_H
#include <meteoio/meteoFilters/FilterBlock.h>
#include <vector>
......@@ -25,38 +25,47 @@
namespace mio {
/**
* @class ProcIIRLowPass
* @class ProcIIR
* @ingroup processing
* @brief Two poles Critically Damped low pass filter.
* @brief Infinite Impulse Response (IIR) filter.
* This filter can either be used as a low pass or high pass filter. It is based on a Critically Damped, 2 poles filter (considering that
* it is better avoid overshooting even at the cost of a gentler falloff). It is possible to use it as a Low Pass (**LP**) or High Pass (**HP**)
*
* 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).
* bidirectional filtering, ie. running the filter twice, first backward and then forward (this also squares the amplitude response). But
* it is possible to disable this bidirectional filtering by adding the "single_pass" argument.
*
* @code
* VW::filter1 = Low_Pass
* VW::arg1 = 10800 ;3 hours
* HS::filter1 = IIR
* HS::arg1 = LP 10800 ;3 hours
* @endcode
*
* It is possible to disable the bidirectional filtering by adding the "single_pass" argument.
* To know more: http://unicorn.us.com/trading/allpolefilters.html and http://www.dspguide.com/ch19/4.htm.
*/
class ProcIIRLowPass : public ProcessingBlock {
class ProcIIR : public ProcessingBlock {
public:
ProcIIRLowPass(const std::vector<std::string>& vec_args, const std::string& name);
ProcIIR(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);
private:
void computeCoefficients(const double& samplerate, const double& f_cutoff, double A[3], double B[3]) const;
typedef enum IIR_TYPE {
BUTTERWORTH,
CRITICALLY_DAMPED,
BESSEL
} IIR_Type;
static void getFilterParameters(const IIR_Type& type, const bool& isLowPass, const double& n, double &g, double &p, double &c);
static double filterPoint(const double& raw_val, const double A[3], const double B[3], std::vector<double> &X, std::vector<double> &Y);
void computeCoefficients(const double& fs, const double& f0, double A[3], double B[3]) const;
void parse_args(std::vector<std::string> vec_args);
double cutoff;
bool bidirectional;
static const double n, p, g, c; ///< filter definition: number of passes, polynomial coefficients, 3dB cutoff correction
double g, p, c; ///< filter definition: number of passes, polynomial coefficients, 3dB cutoff correction
bool bidirectional, low_pass;
};
} //end namespace
......
......@@ -32,7 +32,7 @@
#include <meteoio/meteoFilters/FilterTukey.h>
#include <meteoio/meteoFilters/FilterMAD.h>
#include <meteoio/meteoFilters/ProcAggregate.h>
#include <meteoio/meteoFilters/ProcIIRLowPass.h>
#include <meteoio/meteoFilters/ProcIIR.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
* - LOW_PASS: low pass critically damped filter, see ProcIIRLowPass
* - IIR: Low Pass or High Pass critically damped filter, see ProcIIR
* - 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 == "LOW_PASS"){
return new ProcIIRLowPass(vec_args, blockname);
} else if (blockname == "IIR"){
return new ProcIIR(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