WSL/SLF GitLab Repository

Commit 526dc163 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A low pass Butterworth filter has been implemented as a tentative answer to...

A low pass Butterworth filter has been implemented as a tentative answer to issue 308. Unfortunately, this filter has a non-null phase (ie the signal is shifted in time), meaning that for our applications, it is mostly worthless...
parent dc224f06
......@@ -12,6 +12,7 @@ SET(meteofilters_sources
meteofilters/FilterTukey.cc
meteofilters/FilterRate.cc
meteofilters/FilterUnheatedHNW.cc
meteofilters/ProcButterworth.cc
meteofilters/ProcUnventilatedT.cc
meteofilters/ProcUndercatch_WMO.cc
meteofilters/ProcUndercatch_Hamon.cc
......
/***********************************************************************************/
/* Copyright 2009 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
/***********************************************************************************/
/* This file is part of MeteoIO.
MeteoIO is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
MeteoIO is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
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/MathOptim.h>
#include <cmath>
using namespace std;
namespace mio {
ProcButterworth::ProcButterworth(const std::vector<std::string>& vec_args)
: ProcessingBlock("BUTTERWORTH"), X(3, IOUtils::nodata), Y(3, IOUtils::nodata), cutoff(0.)
{
parse_args(vec_args);
properties.points_before = 2;
properties.stage = ProcessingProperties::first;
}
void ProcButterworth::process(const unsigned int& param, const std::vector<MeteoData>& ivec,
std::vector<MeteoData>& ovec)
{
ovec = ivec;
if(ivec.size()<2 || cutoff==0.) return;
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);
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];
}
}
void ProcButterworth::computeCoefficients(const double& samplerate, const double& f_cutoff, double A[3], double B[3]) const
{
const double QcRaw = (2. * Cst::PI * f_cutoff) / samplerate; // Find cutoff frequency in [0..PI]
const double QcWarp = tan(QcRaw); // Warp cutoff frequency
const double gain = 1. / (1.+Cst::Sqrt2/QcWarp + 2./Optim::pow2(QcWarp));
B[0] = 1.;
B[1] = (2. - 2. * 2./Optim::pow2(QcWarp)) * gain;
B[2] = (1. - Cst::Sqrt2/QcWarp + 2./Optim::pow2(QcWarp)) * gain;
A[0] = 1. * gain;
A[1] = 2. * gain;
A[2] = 1. * gain;
}
void ProcButterworth::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];
}
}
/***********************************************************************************/
/* Copyright 2009 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
/***********************************************************************************/
/* This file is part of MeteoIO.
MeteoIO is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
MeteoIO is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
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__
#include <meteoio/meteofilters/FilterBlock.h>
#include <vector>
#include <string>
namespace mio {
/**
* @class ProcButterworth
* @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. This filter introduces a phase
* (ie the data has a delay compared to the original data). Future implementation might be non-causal in order to remove the phase.
*
* 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::arg1 = 10800 ;3 hours
* @endcode
*
*/
class ProcButterworth : public ProcessingBlock {
public:
ProcButterworth(const std::vector<std::string>& vec_args);
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;
void parse_args(std::vector<std::string> vec_args);
std::vector<double> X, Y;
double cutoff;
};
} //end namespace
#endif
......@@ -27,6 +27,7 @@
#include <meteoio/meteofilters/FilterUnheatedHNW.h>
#include <meteoio/meteofilters/FilterTukey.h>
#include <meteoio/meteofilters/FilterMAD.h>
#include <meteoio/meteofilters/ProcButterworth.h>
#include <meteoio/meteofilters/ProcUndercatch_WMO.h>
#include <meteoio/meteofilters/ProcUndercatch_Hamon.h>
#include <meteoio/meteofilters/ProcUnventilatedT.h>
......@@ -83,7 +84,8 @@ namespace mio {
*
* A few data transformations are also supported besides filtering:
* - EXP_SMOOTHING: exponential smoothing of data, see ProcExpSmoothing
* - WMA_SMOOTHING weighted moving average smoothing of data, see ProcWMASmoothing
* - WMA_SMOOTHING: weighted moving average smoothing of data, see ProcWMASmoothing
* - BUTTERWORTH: low pass butterworth filter, see ProcButterworth
* - MEDIAN_AVG: running median average over a given window, see FilterMedianAvg
* - MEAN_AVG: running mean average over a given window, see FilterMeanAvg
* - WIND_AVG: vector average over a given window, see FilterWindAvg (currently, getting both vw AND dw is broken)
......@@ -110,6 +112,7 @@ bool BlockFactory::initStaticData()
availableBlocks.insert("RATE");
availableBlocks.insert("TUKEY");
availableBlocks.insert("MAD");
availableBlocks.insert("BUTTERWORTH");
availableBlocks.insert("UNHEATED_RAINGAUGE");
availableBlocks.insert("UNDERCATCH_WMO");
availableBlocks.insert("UNDERCATCH_HAMON");
......@@ -147,6 +150,8 @@ ProcessingBlock* BlockFactory::getBlock(const std::string& blockname, const std:
return new FilterTukey(vec_args);
} else if (blockname == "MAD"){
return new FilterMAD(vec_args);
} else if (blockname == "BUTTERWORTH"){
return new ProcButterworth(vec_args);
} else if (blockname == "UNHEATED_RAINGAUGE"){
return new FilterUnheatedHNW(vec_args);
} else if (blockname == "UNDERCATCH_WMO"){
......
......@@ -27,12 +27,6 @@ namespace mio {
class ProcessingProperties {
public:
ProcessingProperties() : time_before(0., 0.), time_after(0., 0.),
points_before(0), points_after(0),
stage(first) {}
const std::string toString() const;
typedef enum PROC_STAGE { none, ///< never activate this block
first, ///< activate at first stage
second, ///< activate at second stage
......@@ -40,6 +34,12 @@ class ProcessingProperties {
//once ///< activate at stage one or two, but only once
} proc_stage;
ProcessingProperties() : time_before(0., 0.), time_after(0., 0.),
points_before(0), points_after(0),
stage(first) {}
const std::string toString() const;
Duration time_before;
Duration time_after;
......
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