← Back to branch summary

~ubuntu-branches/ubuntu/wily/openms/wily

« back to all changes in this revision

Viewing changes to include/OpenMS/FILTERING/SMOOTHING/SavitzkyGolayFilter.h

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2012-11-12 15:58:12 UTC
  • Revision ID: package-import@ubuntu.com-20121112155812-vr15wtg9b50cuesg
Tags: upstream-1.9.0
Import upstream version 1.9.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/OpenMS/FILTERING/SMOOTHING/SavitzkyGolayFilter.h
 
1
// -*- mode: C++; tab-width: 2; -*-
 
2
// vi: set ts=2:
 
3
//
 
4
// --------------------------------------------------------------------------
 
5
//                   OpenMS Mass Spectrometry Framework
 
6
// --------------------------------------------------------------------------
 
7
//  Copyright (C) 2003-2011 -- Oliver Kohlbacher, Knut Reinert
 
8
//
 
9
//  This library is free software; you can redistribute it and/or
 
10
//  modify it under the terms of the GNU Lesser General Public
 
11
//  License as published by the Free Software Foundation; either
 
12
//  version 2.1 of the License, or (at your option) any later version.
 
13
//
 
14
//  This library is distributed in the hope that it will be useful,
 
15
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
16
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
17
//  Lesser General Public License for more details.
 
18
//
 
19
//  You should have received a copy of the GNU Lesser General Public
 
20
//  License along with this library; if not, write to the Free Software
 
21
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
22
//
 
23
// --------------------------------------------------------------------------
 
24
// $Maintainer: Alexandra Zerck $
 
25
// $Authors: Eva Lange $
 
26
// --------------------------------------------------------------------------
 
27
 
 
28
#ifndef OPENMS_FILTERING_SMOOTHING_SAVITZKYGOLAYFILTER_H
 
29
#define OPENMS_FILTERING_SMOOTHING_SAVITZKYGOLAYFILTER_H
 
30
 
 
31
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
 
32
#include <OpenMS/CONCEPT/ProgressLogger.h>
 
33
#include <OpenMS/KERNEL/MSExperiment.h>
 
34
 
 
35
#include <gsl/gsl_vector.h>
 
36
#include <gsl/gsl_matrix.h>
 
37
#include <gsl/gsl_linalg.h>
 
38
#include <gsl/gsl_permutation.h>
 
39
#include <gsl/gsl_pow_int.h>
 
40
 
 
41
namespace OpenMS
 
42
{
 
43
  /**
 
44
    @brief Computes the Savitzky-Golay filter coefficients using QR decomposition.
 
45
   
 
46
    This class represents a Savitzky-Golay lowpass-filter. The idea of the Savitzky-Golay filter
 
47
    is to find filtercoefficients that preserve higher moments, which means to approximate the underlying
 
48
    function within the moving window by a polynomial of higher order (typically quadratic or quartic).
 
49
    Therefore we least-squares fit for each data point a polynomial to all points \f$ f_i \f$ in the window
 
50
    and set \f$g_i\f$ to be the value of that polynomial at position \f$ i \f$. This method is superior
 
51
    to adjacent averaging because it tends to preserve features of the data such as peak height and width, which
 
52
    are usually 'washed out' by adjacent averaging.
 
53
   
 
54
    Because of the linearity of the problem, we can reduce the work computing by fitting in advance, for fictious data
 
55
    consisiting of all zeros except for a singe 1 and then do the fits on the real data just by taking linear
 
56
    combinations. There are a particular sets of filter coefficients \f$ c_n \f$ which accomplish the process of
 
57
    polynomial least-squares fit inside a moving window. To get the symmetric coefficient-matrix
 
58
    \f$C \in R^{frameSize \times frameSize}\f$ with
 
59
    \f[ C= \left[ \begin{array}{cccc} c_{0,0} & c_{0,1} & \cdots & c_{0,frameSize-1} \\
 
60
                               \vdots  &         &         & \vdots            \\
 
61
             c_{frameSize-1,0} & c_{frameSize-1,2} & \ldots  & c_{frameSize-1,frameSize-1} \end{array} \right]\f]
 
62
    The first (last) \f$ \frac{frameSize}{2} \f$ rows of \f$ C \f$ we need to smooth the first (last)
 
63
    \f$ frameSize \f$ data points of the signal. So we use for the smoothing of the first data point the data
 
64
    point itself and the next \f$ frameSize-1 \f$ future points. For the second point we take the first datapoint,
 
65
    the data point itself and \f$ frameSize-2 \f$ of rightward data points... .
 
66
    We compute the Matrix \f$ C \f$ by solving the underlying least-squares problems with the singular value decomposition.
 
67
    Here we demonstrate the computation of the first row of a coefficient-matrix \f$ C \f$ for a Savitzky-Golay Filter
 
68
    of order=3 and frameSize=5:
 
69
    The design-matrix for the least-squares fit of a linear combination of 3 basis functions to 5 data points is:
 
70
    \f[ A=\left[ \begin{array}{ccc} x_0^0 & x_0^1 & x_0^2 \\ \\
 
71
                             x_1^0 & x_1^1 & x_1^2 \\ \\
 
72
           x_2^0 & x_2^1 & x_2^2 \\ \\
 
73
           x_3^0 & x_3^1 & x_3^2 \\ \\
 
74
           x_4^0 & x_4^1 & x_4^2 \end{array} \right]. \f]
 
75
    To smooth the first data point we have to create a design-matrix with \f$ x=[0,\ldots, frameSize-1] \f$.
 
76
    Now we have to solve the over-determined set of \f$ frameSize \f$ linear equations
 
77
    \f[ Ac=b \f]
 
78
    where \f$ b=[1,0,\ldots,0] \f$ represents the fictious data.
 
79
    Therefore we solve the normal equations of the least-squares problem
 
80
    \f[ A^TAc=A^Tb. \f]
 
81
    Now, it is possible to get
 
82
    \f[ c_n=\sum_{m=0}^8 \{(A^TA)^{-1}\}_{0,m} n^m, \f]
 
83
    with \f$ 0\le n \le 8\f$. Because we only need one row of the inverse matrix, it is possible to use LU decomposition with
 
84
    only a single backsubstitution.
 
85
    The vector \f$c=[c_0,\ldots,c_8] \f$ represents the wanted coefficients.
 
86
    Note that the solution of a least-squares problem directly from the normal equations is faster than the singular value
 
87
    decomposition but rather susceptible to roundoff error!
 
88
   
 
89
    @note This filter works only for uniform profile data!
 
90
          A polynom order of 4 is recommended.
 
91
          The bigger the frame size the smoother the signal (the more detail information get lost!). The frame size corresponds to the number
 
92
          of filter coefficients, so the width of the smoothing intervall is given by frame_size*spacing of the profile data.
 
93
    
 
94
                @note The data must be sorted according to ascending m/z!
 
95
    
 
96
                @htmlinclude OpenMS_SavitzkyGolayFilter.parameters
 
97
    
 
98
    @ingroup SignalProcessing
 
99
  */
 
100
  class OPENMS_DLLAPI SavitzkyGolayFilter 
 
101
        : public ProgressLogger, 
 
102
                public DefaultParamHandler
 
103
  {
 
104
    public:
 
105
      /// Constructor
 
106
      SavitzkyGolayFilter();
 
107
 
 
108
      /// Destructor
 
109
      virtual ~SavitzkyGolayFilter();
 
110
 
 
111
      /** 
 
112
        @brief Removed the noise from an MSSpectrum containing profile data.
 
113
      */
 
114
      template <typename PeakType>
 
115
      void filter(MSSpectrum<PeakType>& spectrum)
 
116
      {
 
117
        UInt n = (UInt)spectrum.size();
 
118
        
 
119
        typename MSSpectrum<PeakType>::iterator first = spectrum.begin();
 
120
        typename MSSpectrum<PeakType>::iterator last = spectrum.end();
 
121
        
 
122
        //copy the data AND META DATA to the output container
 
123
        MSSpectrum<PeakType> output = spectrum;
 
124
        
 
125
        if (frame_size_ > n)
 
126
        {
 
127
          return;
 
128
        }
 
129
 
 
130
        int i;
 
131
        UInt j;
 
132
        int mid=(frame_size_/2);
 
133
        double help;
 
134
 
 
135
        typename MSSpectrum<PeakType>::iterator it_forward;
 
136
        typename MSSpectrum<PeakType>::iterator it_help;
 
137
        typename MSSpectrum<PeakType>::iterator out_it = output.begin();
 
138
 
 
139
        // compute the transient on
 
140
        for (i=0; i <= mid; ++i)
 
141
        {
 
142
          it_forward=(first-i);
 
143
          help=0;
 
144
 
 
145
          for (j=0; j < frame_size_; ++j)
 
146
          {
 
147
            help+=it_forward->getIntensity()*coeffs_[(i+1)*frame_size_-1-j];
 
148
            ++it_forward;
 
149
          }
 
150
 
 
151
 
 
152
          out_it->setPosition(first->getPosition());
 
153
          out_it->setIntensity(std::max(0.0,help));
 
154
          ++out_it;
 
155
          ++first;
 
156
        }
 
157
 
 
158
        // compute the steady state output
 
159
        it_help=(last-mid);
 
160
        while (first!=it_help)
 
161
        {
 
162
          it_forward=(first-mid);
 
163
          help=0;
 
164
 
 
165
          for (j=0; j < frame_size_; ++j)
 
166
          {
 
167
            help+=it_forward->getIntensity()*coeffs_[mid*frame_size_+j];
 
168
            ++it_forward;
 
169
          }
 
170
 
 
171
 
 
172
          out_it->setPosition(first->getPosition());
 
173
          out_it->setIntensity(std::max(0.0,help));
 
174
          ++out_it;
 
175
          ++first;
 
176
        }
 
177
 
 
178
        // compute the transient off
 
179
        for (i=(mid-1); i >= 0; --i)
 
180
        {
 
181
          it_forward=(first-(frame_size_-i-1));
 
182
          help=0;
 
183
 
 
184
          for (j=0; j < frame_size_; ++j)
 
185
          {
 
186
            help+=it_forward->getIntensity()*coeffs_[i*frame_size_+j];
 
187
            ++it_forward;
 
188
          }
 
189
 
 
190
          out_it->setPosition(first->getPosition());
 
191
          out_it->setIntensity(std::max(0.0,help));
 
192
          ++out_it;
 
193
          ++first;
 
194
        }
 
195
        
 
196
        spectrum = output;
 
197
      }
 
198
 
 
199
      /** 
 
200
        @brief Removed the noise from an MSExperiment containing profile data.
 
201
      */
 
202
      template <typename PeakType>
 
203
      void filterExperiment(MSExperiment<PeakType>& map)
 
204
      {
 
205
        startProgress(0,map.size(),"smoothing data");
 
206
        for (Size i = 0; i < map.size(); ++i)
 
207
        {
 
208
          filter(map[i]);
 
209
          setProgress(i);
 
210
        }
 
211
        endProgress();
 
212
      }
 
213
          
 
214
    protected:
 
215
                        /// Coefficients
 
216
                        std::vector<DoubleReal> coeffs_;
 
217
                        /// UInt of the filter kernel (number of pre-tabulated coefficients)
 
218
      UInt frame_size_;
 
219
      /// The order of the smoothing polynomial.
 
220
      UInt order_;
 
221
                        // Docu in base class
 
222
      virtual void updateMembers_();
 
223
      
 
224
  };
 
225
 
 
226
} // namespace OpenMS
 
227
#endif // OPENMS_FILTERING_SMOOTHING_SAVITZKYGOLAYFILTER_H
 
228