savgolfilter.h

Go to the documentation of this file.

/* BEGIN software license
 *
 * msXpertSuite - mass spectrometry software suite
 * -----------------------------------------------
 * Copyright(C) 2009,...,2018 Filippo Rusconi
 *
 * http://www.msxpertsuite.org
 *
 * This file is part of the msXpertSuite project.
 *
 * The msXpertSuite project is the successor of the massXpert project. This
 * project now includes various independent modules:
 *
 * - massXpert, model polymer chemistries and simulate mass spectrometric data;
 * - mineXpert, a powerful TIC chromatogram/mass spectrum viewer/miner;
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 *
 * END software license
 */
 
 
#pragma once
 
 
#include <QObject>
 
#include "../../trace/trace.h"
#include "pappsomspp/export-import-config.h"
#include "filternameinterface.h"
 
 
namespace pappso
{
 
 
//! Parameters for the Savitzky-Golay filter
struct PMSPP_LIB_DECL SavGolParams
{
  int nL = 15;
  //!< number of data points on the left of the filtered point
  int nR = 15;
  //!< number of data points on the right of the filtered point
  int m = 4;
  //!< order of the polynomial to use in the regression analysis leading to the
  //! Savitzky-Golay coefficients (typically between 2 and 6)
  int lD = 0;
  //!< specifies the order of the derivative to extract from the Savitzky-Golay
  //! smoothing algorithm (for regular smoothing, use 0)
  bool convolveWithNr = false;
  //!< set to false for best results
 
  SavGolParams(){};
 
  SavGolParams(const SavGolParams &other)
    : nL{other.nL},
      nR{other.nR},
      m{other.m},
      lD{other.lD},
      convolveWithNr{other.convolveWithNr}
  {
  }
 
  SavGolParams(
    int nLParam, int nRParam, int mParam, int lDParam, bool convolveWithNrParam)
  {
    nL             = nLParam;
    nR             = nRParam;
    m              = mParam;
    lD             = lDParam;
    convolveWithNr = convolveWithNrParam;
  }
 
  void
  initialize(
    int nLParam, int nRParam, int mParam, int lDParam, bool convolveWithNrParam)
  {
    nL             = nLParam;
    nR             = nRParam;
    m              = mParam;
    lD             = lDParam;
    convolveWithNr = convolveWithNrParam;
  }
 
  void
  initialize(const SavGolParams &other)
  {
    nL             = other.nL;
    nR             = other.nR;
    m              = other.m;
    lD             = other.lD;
    convolveWithNr = other.convolveWithNr;
  }
 
  void
  initialize(const QString &parameters)
  {
    // Typical string: "15;15;4;0;false"
 
    QStringList params_list = parameters.split(";");
 
    nL             = params_list.at(0).toInt();
    nR             = params_list.at(1).toInt();
    m              = params_list.at(2).toInt();
    lD             = params_list.at(3).toInt();
    convolveWithNr = (params_list.at(4) == "true" ? true : false);
  }
 
  QString
  toString() const
  {
    // Typical string: "15;15;4;0;false"
    return QString("%1;%2;%3;%4;%5")
      .arg(QString::number(nL))
      .arg(QString::number(nR))
      .arg(QString::number(m))
      .arg(QString::number(lD))
      .arg(convolveWithNr ? "true" : "false");
  }
};
 
 
class FilterSavitzkyGolay;
 
typedef std::shared_ptr<FilterSavitzkyGolay> FilterSavitzkyGolaySPtr;
typedef std::shared_ptr<const FilterSavitzkyGolay> FilterSavitzkyGolayCstSPtr;
 
 
/**
 * @brief uses Savitsky-Golay filter on trace
 */
class PMSPP_LIB_DECL FilterSavitzkyGolay : public FilterNameInterface
{
  public:
   FilterSavitzkyGolay();
  /** Construct a FilterSavitzkyGolay instance using the Savitzky-Golay
   parameters \param nL number of data point left of the point being filtered
   \param nR number of data point right of the point being filtered
   \param m order of the polynomial to use in the regression analysis
   \param lD order of the derivative to extract
   \param convolveWithNr set to false
   */
  FilterSavitzkyGolay(
    int nL, int nR, int m, int lD, bool convolveWithNr = false);
 
  FilterSavitzkyGolay(const SavGolParams sav_gol_params);
 
  FilterSavitzkyGolay(const QString &parameters);
 
  /**
   * Copy constructor
   *
   * @param other TODO
   */
  FilterSavitzkyGolay(const FilterSavitzkyGolay &other);
 
  /**
   * Destructor
   */
  virtual ~FilterSavitzkyGolay();
 
  FilterSavitzkyGolay &operator=(const FilterSavitzkyGolay &other);
 
  Trace &filter(Trace &data_points) const override;
 
  SavGolParams getParameters() const;
 
  char runFilter(double *y_data_p,
                 double *y_filtered_data_p,
                 int data_point_count) const;
 
  void filteredData(std::vector<pappso_double> &data);
 
  QString name() const override;
 
  // Utility function
 
  QString toString() const override;
 
  protected:
  void buildFilterFromString(const QString &strBuildParams) override;
 
  private:
  ///// Parameters to configure the Savitzky-Golay filter algorithm
  // The default values in m_params work fine for mass spectra.
  SavGolParams m_params;
 
  ///// Data used for running the algorithm.
 
  //! C array of keys of the Trace
  pappso_double *m_x;
 
  //! C array of raw values of the Trace
  pappso_double *m_yr;
 
  //! C array of filtered values after the computation has been performed
  pappso_double *m_yf;
 
 
  ///// Functions that actually implement the algorithm.
  int *ivector(long nl, long nh) const;
  pappso_double *dvector(long nl, long nh) const;
  pappso_double **dmatrix(long nrl, long nrh, long ncl, long nch) const;
  void free_ivector(int *v, long nl, long nh) const;
  void free_dvector(pappso_double *v, long nl, long nh) const;
  void
  free_dmatrix(pappso_double **m, long nrl, long nrh, long ncl, long nch) const;
  void lubksb(pappso_double **a, int n, int *indx, pappso_double b[]) const;
  void ludcmp(pappso_double **a, int n, int *indx, pappso_double *d) const;
  void four1(pappso_double data[], unsigned long nn, int isign);
  void twofft(pappso_double data1[],
              pappso_double data2[],
              pappso_double fft1[],
              pappso_double fft2[],
              unsigned long n);
  void realft(pappso_double data[], unsigned long n, int isign);
  char convlv(pappso_double data[],
              unsigned long n,
              pappso_double respns[],
              unsigned long m,
              int isign,
              pappso_double ans[]);
  char sgcoeff(pappso_double c[], int np, int nl, int nr, int ld, int m) const;
};
} // namespace pappso