libpappsomspp
Library for mass spectrometry
savgolfilter.cpp
Go to the documentation of this file.
1 /* BEGIN software license
2  *
3  * msXpertSuite - mass spectrometry software suite
4  * -----------------------------------------------
5  * Copyright(C) 2009,...,2018 Filippo Rusconi
6  *
7  * http://www.msxpertsuite.org
8  *
9  * This file is part of the msXpertSuite project.
10  *
11  * The msXpertSuite project is the successor of the massXpert project. This
12  * project now includes various independent modules:
13  *
14  * - massXpert, model polymer chemistries and simulate mass spectrometric data;
15  * - mineXpert, a powerful TIC chromatogram/mass spectrum viewer/miner;
16  *
17  * This program is free software: you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation, either version 3 of the License, or
20  * (at your option) any later version.
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program. If not, see <http://www.gnu.org/licenses/>.
29  *
30  * END software license
31  */
32 
33 
34 #include <qmath.h>
35 
36 
37 #include <QVector>
38 #include <QDebug>
39 
40 
41 #include "savgolfilter.h"
42 
43 
44 namespace pappso
45 {
46 
47 
48 #define SWAP(a, b) \
49  tempr = (a); \
50  (a) = (b); \
51  (b) = tempr
52 
53 
55  int nL, int nR, int m, int lD, bool convolveWithNr)
56 {
57  m_nL = nL;
58  m_nR = nR;
59  m_m = m;
60  m_lD = lD;
61  m_convolveWithNr = convolveWithNr;
62 }
63 
64 
66 {
67  // This function only copies the parameters, not the data.
68 
69  m_nL = other.m_nL;
70  m_nR = other.m_nR;
71  m_m = other.m_m;
72  m_lD = other.m_lD;
74 }
75 
76 
78 {
79 }
80 
83 {
84  if(&other == this)
85  return *this;
86 
87  // This function only copies the parameters, not the data.
88 
89  m_nL = other.m_nL;
90  m_nR = other.m_nR;
91  m_m = other.m_m;
92  m_lD = other.m_lD;
94 
95  return *this;
96 }
97 
98 
99 Trace &
101 {
102  // Initialize data:
103 
104  // We want the filter to stay constant so we create a local copy of the data.
105 
106  int data_point_count = data_points.size();
107 
108  pappso_double *x_data_p = dvector(1, data_point_count);
109  pappso_double *y_initial_data_p = dvector(1, data_point_count);
110  pappso_double *y_filtered_data_p = nullptr;
111 
112  if(m_convolveWithNr)
113  y_filtered_data_p = dvector(1, 2 * data_point_count);
114  else
115  y_filtered_data_p = dvector(1, data_point_count);
116 
117  for(int iter = 0; iter < data_point_count; ++iter)
118  {
119  x_data_p[iter] = data_points.at(iter).x;
120  y_initial_data_p[iter] = data_points.at(iter).y;
121  }
122 
123  // Now run the filter.
124 
125  runFilter(y_initial_data_p, y_filtered_data_p, data_point_count);
126 
127  // Put back the modified y values into the trace.
128  auto iter_yf = y_filtered_data_p;
129  for(auto &data_point : data_points)
130  {
131  data_point.y = *iter_yf;
132  iter_yf++;
133  }
134 
135  return data_points;
136 }
137 
138 
141 {
143 }
144 
145 
146 //! Return a string with the textual representation of the configuration data.
147 QString
149 {
150  return QString::asprintf(
151  "Savitzy-Golay filter parameters:\n"
152  "nL: %d ; nR: %d ; m: %d ; lD: %d ; convolveWithNr : %s",
153  m_nL,
154  m_nR,
155  m_m,
156  m_lD,
157  (m_convolveWithNr ? "true" : "false"));
158 }
159 
160 
161 int *
162 FilterSavitzkyGolay::ivector(long nl, long nh) const
163 {
164  int *v;
165  v = (int *)malloc((size_t)((nh - nl + 2) * sizeof(int)));
166  if(!v)
167  {
168  qFatal("Error: Allocation failure.");
169  }
170  return v - nl + 1;
171 }
172 
173 
175 FilterSavitzkyGolay::dvector(long nl, long nh) const
176 {
177  pappso_double *v;
178  long k;
179  v = (pappso_double *)malloc((size_t)((nh - nl + 2) * sizeof(pappso_double)));
180  if(!v)
181  {
182  qFatal("Error: Allocation failure.");
183  }
184  for(k = nl; k <= nh; k++)
185  v[k] = 0.0;
186  return v - nl + 1;
187 }
188 
189 
190 pappso_double **
191 FilterSavitzkyGolay::dmatrix(long nrl, long nrh, long ncl, long nch) const
192 {
193  long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
194  pappso_double **m;
195  m = (pappso_double **)malloc((size_t)((nrow + 1) * sizeof(pappso_double *)));
196  if(!m)
197  {
198  qFatal("Error: Allocation failure.");
199  }
200  m += 1;
201  m -= nrl;
202  m[nrl] = (pappso_double *)malloc(
203  (size_t)((nrow * ncol + 1) * sizeof(pappso_double)));
204  if(!m[nrl])
205  {
206  qFatal("Error: Allocation failure.");
207  }
208  m[nrl] += 1;
209  m[nrl] -= ncl;
210  for(i = nrl + 1; i <= nrh; i++)
211  m[i] = m[i - 1] + ncol;
212  return m;
213 }
214 
215 
216 void
217 FilterSavitzkyGolay::free_ivector(int *v, long nl, [[maybe_unused]] long nh) const
218 {
219  free((char *)(v + nl - 1));
220 }
221 
222 
223 void
224 FilterSavitzkyGolay::free_dvector(pappso_double *v, long nl, [[maybe_unused]] long nh) const
225 {
226  free((char *)(v + nl - 1));
227 }
228 
229 
230 void
232  pappso_double **m, long nrl, [[maybe_unused]] long nrh, long ncl, [[maybe_unused]] long nch) const
233 {
234  free((char *)(m[nrl] + ncl - 1));
235  free((char *)(m + nrl - 1));
236 }
237 
238 
239 void
241  int n,
242  int *indx,
243  pappso_double b[]) const
244 {
245  int i, ii = 0, ip, j;
247 
248  for(i = 1; i <= n; i++)
249  {
250  ip = indx[i];
251  sum = b[ip];
252  b[ip] = b[i];
253  if(ii)
254  for(j = ii; j <= i - 1; j++)
255  sum -= a[i][j] * b[j];
256  else if(sum)
257  ii = i;
258  b[i] = sum;
259  }
260  for(i = n; i >= 1; i--)
261  {
262  sum = b[i];
263  for(j = i + 1; j <= n; j++)
264  sum -= a[i][j] * b[j];
265  b[i] = sum / a[i][i];
266  }
267 }
268 
269 
270 void
272  int n,
273  int *indx,
274  pappso_double *d) const
275 {
276  int i, imax = 0, j, k;
277  pappso_double big, dum, sum, temp;
278  pappso_double *vv;
279 
280  vv = dvector(1, n);
281  *d = 1.0;
282  for(i = 1; i <= n; i++)
283  {
284  big = 0.0;
285  for(j = 1; j <= n; j++)
286  if((temp = fabs(a[i][j])) > big)
287  big = temp;
288  if(big == 0.0)
289  {
290  qFatal("Error: Singular matrix found in routine ludcmp().");
291  }
292  vv[i] = 1.0 / big;
293  }
294  for(j = 1; j <= n; j++)
295  {
296  for(i = 1; i < j; i++)
297  {
298  sum = a[i][j];
299  for(k = 1; k < i; k++)
300  sum -= a[i][k] * a[k][j];
301  a[i][j] = sum;
302  }
303  big = 0.0;
304  for(i = j; i <= n; i++)
305  {
306  sum = a[i][j];
307  for(k = 1; k < j; k++)
308  sum -= a[i][k] * a[k][j];
309  a[i][j] = sum;
310  if((dum = vv[i] * fabs(sum)) >= big)
311  {
312  big = dum;
313  imax = i;
314  }
315  }
316  if(j != imax)
317  {
318  for(k = 1; k <= n; k++)
319  {
320  dum = a[imax][k];
321  a[imax][k] = a[j][k];
322  a[j][k] = dum;
323  }
324  *d = -(*d);
325  vv[imax] = vv[j];
326  }
327  indx[j] = imax;
328  if(a[j][j] == 0.0)
329  a[j][j] = std::numeric_limits<pappso_double>::epsilon();
330  if(j != n)
331  {
332  dum = 1.0 / (a[j][j]);
333  for(i = j + 1; i <= n; i++)
334  a[i][j] *= dum;
335  }
336  }
337  free_dvector(vv, 1, n);
338 }
339 
340 
341 void
342 FilterSavitzkyGolay::four1(pappso_double data[], unsigned long nn, int isign)
343 {
344  unsigned long n, mmax, m, j, istep, i;
345  pappso_double wtemp, wr, wpr, wpi, wi, theta;
346  pappso_double tempr, tempi;
347 
348  n = nn << 1;
349  j = 1;
350  for(i = 1; i < n; i += 2)
351  {
352  if(j > i)
353  {
354  SWAP(data[j], data[i]);
355  SWAP(data[j + 1], data[i + 1]);
356  }
357  m = n >> 1;
358  while(m >= 2 && j > m)
359  {
360  j -= m;
361  m >>= 1;
362  }
363  j += m;
364  }
365  mmax = 2;
366  while(n > mmax)
367  {
368  istep = mmax << 1;
369  theta = isign * (6.28318530717959 / mmax);
370  wtemp = sin(0.5 * theta);
371  wpr = -2.0 * wtemp * wtemp;
372  wpi = sin(theta);
373  wr = 1.0;
374  wi = 0.0;
375  for(m = 1; m < mmax; m += 2)
376  {
377  for(i = m; i <= n; i += istep)
378  {
379  j = i + mmax;
380  tempr = wr * data[j] - wi * data[j + 1];
381  tempi = wr * data[j + 1] + wi * data[j];
382  data[j] = data[i] - tempr;
383  data[j + 1] = data[i + 1] - tempi;
384  data[i] += tempr;
385  data[i + 1] += tempi;
386  }
387  wr = (wtemp = wr) * wpr - wi * wpi + wr;
388  wi = wi * wpr + wtemp * wpi + wi;
389  }
390  mmax = istep;
391  }
392 }
393 
394 
395 void
397  pappso_double data2[],
398  pappso_double fft1[],
399  pappso_double fft2[],
400  unsigned long n)
401 {
402  unsigned long nn3, nn2, jj, j;
403  pappso_double rep, rem, aip, aim;
404 
405  nn3 = 1 + (nn2 = 2 + n + n);
406  for(j = 1, jj = 2; j <= n; j++, jj += 2)
407  {
408  fft1[jj - 1] = data1[j];
409  fft1[jj] = data2[j];
410  }
411  four1(fft1, n, 1);
412  fft2[1] = fft1[2];
413  fft1[2] = fft2[2] = 0.0;
414  for(j = 3; j <= n + 1; j += 2)
415  {
416  rep = 0.5 * (fft1[j] + fft1[nn2 - j]);
417  rem = 0.5 * (fft1[j] - fft1[nn2 - j]);
418  aip = 0.5 * (fft1[j + 1] + fft1[nn3 - j]);
419  aim = 0.5 * (fft1[j + 1] - fft1[nn3 - j]);
420  fft1[j] = rep;
421  fft1[j + 1] = aim;
422  fft1[nn2 - j] = rep;
423  fft1[nn3 - j] = -aim;
424  fft2[j] = aip;
425  fft2[j + 1] = -rem;
426  fft2[nn2 - j] = aip;
427  fft2[nn3 - j] = rem;
428  }
429 }
430 
431 
432 void
433 FilterSavitzkyGolay::realft(pappso_double data[], unsigned long n, int isign)
434 {
435  unsigned long i, i1, i2, i3, i4, np3;
436  pappso_double c1 = 0.5, c2, h1r, h1i, h2r, h2i;
437  pappso_double wr, wi, wpr, wpi, wtemp, theta;
438 
439  theta = 3.141592653589793 / (pappso_double)(n >> 1);
440  if(isign == 1)
441  {
442  c2 = -0.5;
443  four1(data, n >> 1, 1);
444  }
445  else
446  {
447  c2 = 0.5;
448  theta = -theta;
449  }
450  wtemp = sin(0.5 * theta);
451  wpr = -2.0 * wtemp * wtemp;
452  wpi = sin(theta);
453  wr = 1.0 + wpr;
454  wi = wpi;
455  np3 = n + 3;
456  for(i = 2; i <= (n >> 2); i++)
457  {
458  i4 = 1 + (i3 = np3 - (i2 = 1 + (i1 = i + i - 1)));
459  h1r = c1 * (data[i1] + data[i3]);
460  h1i = c1 * (data[i2] - data[i4]);
461  h2r = -c2 * (data[i2] + data[i4]);
462  h2i = c2 * (data[i1] - data[i3]);
463  data[i1] = h1r + wr * h2r - wi * h2i;
464  data[i2] = h1i + wr * h2i + wi * h2r;
465  data[i3] = h1r - wr * h2r + wi * h2i;
466  data[i4] = -h1i + wr * h2i + wi * h2r;
467  wr = (wtemp = wr) * wpr - wi * wpi + wr;
468  wi = wi * wpr + wtemp * wpi + wi;
469  }
470  if(isign == 1)
471  {
472  data[1] = (h1r = data[1]) + data[2];
473  data[2] = h1r - data[2];
474  }
475  else
476  {
477  data[1] = c1 * ((h1r = data[1]) + data[2]);
478  data[2] = c1 * (h1r - data[2]);
479  four1(data, n >> 1, -1);
480  }
481 }
482 
483 
484 char
486  unsigned long n,
487  pappso_double respns[],
488  unsigned long m,
489  int isign,
490  pappso_double ans[])
491 {
492  unsigned long i, no2;
493  pappso_double dum, mag2, *fft;
494 
495  fft = dvector(1, n << 1);
496  for(i = 1; i <= (m - 1) / 2; i++)
497  respns[n + 1 - i] = respns[m + 1 - i];
498  for(i = (m + 3) / 2; i <= n - (m - 1) / 2; i++)
499  respns[i] = 0.0;
500  twofft(data, respns, fft, ans, n);
501  no2 = n >> 1;
502  for(i = 2; i <= n + 2; i += 2)
503  {
504  if(isign == 1)
505  {
506  ans[i - 1] =
507  (fft[i - 1] * (dum = ans[i - 1]) - fft[i] * ans[i]) / no2;
508  ans[i] = (fft[i] * dum + fft[i - 1] * ans[i]) / no2;
509  }
510  else if(isign == -1)
511  {
512  if((mag2 = ans[i - 1] * ans[i - 1] + ans[i] * ans[i]) == 0.0)
513  {
514  qDebug("Attempt of deconvolving at zero response in convlv().");
515  return (1);
516  }
517  ans[i - 1] =
518  (fft[i - 1] * (dum = ans[i - 1]) + fft[i] * ans[i]) / mag2 / no2;
519  ans[i] = (fft[i] * dum - fft[i - 1] * ans[i]) / mag2 / no2;
520  }
521  else
522  {
523  qDebug("No meaning for isign in convlv().");
524  return (1);
525  }
526  }
527  ans[2] = ans[n + 1];
528  realft(ans, n, -1);
529  free_dvector(fft, 1, n << 1);
530  return (0);
531 }
532 
533 
534 char
536  pappso_double c[], int np, int nl, int nr, int ld, int m) const
537 {
538  int imj, ipj, j, k, kk, mm, *indx;
539  pappso_double d, fac, sum, **a, *b;
540 
541  if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m)
542  {
543  qDebug("Inconsistent arguments detected in routine sgcoeff.");
544  return (1);
545  }
546  indx = ivector(1, m + 1);
547  a = dmatrix(1, m + 1, 1, m + 1);
548  b = dvector(1, m + 1);
549  for(ipj = 0; ipj <= (m << 1); ipj++)
550  {
551  sum = (ipj ? 0.0 : 1.0);
552  for(k = 1; k <= nr; k++)
553  sum += pow((pappso_double)k, (pappso_double)ipj);
554  for(k = 1; k <= nl; k++)
555  sum += pow((pappso_double)-k, (pappso_double)ipj);
556  mm = (ipj < 2 * m - ipj ? ipj : 2 * m - ipj);
557  for(imj = -mm; imj <= mm; imj += 2)
558  a[1 + (ipj + imj) / 2][1 + (ipj - imj) / 2] = sum;
559  }
560  ludcmp(a, m + 1, indx, &d);
561  for(j = 1; j <= m + 1; j++)
562  b[j] = 0.0;
563  b[ld + 1] = 1.0;
564  lubksb(a, m + 1, indx, b);
565  for(kk = 1; kk <= np; kk++)
566  c[kk] = 0.0;
567  for(k = -nl; k <= nr; k++)
568  {
569  sum = b[1];
570  fac = 1.0;
571  for(mm = 1; mm <= m; mm++)
572  sum += b[mm + 1] * (fac *= k);
573  kk = ((np - k) % np) + 1;
574  c[kk] = sum;
575  }
576  free_dvector(b, 1, m + 1);
577  free_dmatrix(a, 1, m + 1, 1, m + 1);
578  free_ivector(indx, 1, m + 1);
579  return (0);
580 }
581 
582 
583 //! Perform the Savitzky-Golay filtering process.
584 /*
585  The results are in the \c y_filtered_data_p C array of pappso_double values.
586  */
587 char
589  double *y_filtered_data_p,
590  int data_point_count) const
591 {
592  int np = m_nL + 1 + m_nR;
593  pappso_double *c;
594  char retval;
595 
596 #if CONVOLVE_WITH_NR_CONVLV
597  c = dvector(1, data_point_count);
598  retval = sgcoeff(c, np, m_nL, m_nR, m_lD, m_m);
599  if(retval == 0)
600  convlv(y_data_p, data_point_count, c, np, 1, y_filtered_data_p);
601  free_dvector(c, 1, data_point_count);
602 #else
603  int j;
604  long int k;
605  c = dvector(1, m_nL + m_nR + 1);
606  retval = sgcoeff(c, np, m_nL, m_nR, m_lD, m_m);
607  if(retval == 0)
608  {
609  qDebug() << __FILE__ << __LINE__ << __FUNCTION__ << "()"
610  << "retval is 0";
611 
612  for(k = 1; k <= m_nL; k++)
613  {
614  for(y_filtered_data_p[k] = 0.0, j = -m_nL; j <= m_nR; j++)
615  {
616  if(k + j >= 1)
617  {
618  y_filtered_data_p[k] +=
619  c[(j >= 0 ? j + 1 : m_nR + m_nL + 2 + j)] * y_data_p[k + j];
620  }
621  }
622  }
623  for(k = m_nL + 1; k <= data_point_count - m_nR; k++)
624  {
625  for(y_filtered_data_p[k] = 0.0, j = -m_nL; j <= m_nR; j++)
626  {
627  y_filtered_data_p[k] +=
628  c[(j >= 0 ? j + 1 : m_nR + m_nL + 2 + j)] * y_data_p[k + j];
629  }
630  }
631  for(k = data_point_count - m_nR + 1; k <= data_point_count; k++)
632  {
633  for(y_filtered_data_p[k] = 0.0, j = -m_nL; j <= m_nR; j++)
634  {
635  if(k + j <= data_point_count)
636  {
637  y_filtered_data_p[k] +=
638  c[(j >= 0 ? j + 1 : m_nR + m_nL + 2 + j)] * y_data_p[k + j];
639  }
640  }
641  }
642  }
643 
644  free_dvector(c, 1, m_nR + m_nL + 1);
645 #endif
646 
647  return (retval);
648 }
649 
650 
651 } // namespace pappso
pappso::FilterSavitzkyGolay::runFilter
char runFilter(double *y_data_p, double *y_filtered_data_p, int data_point_count) const
Perform the Savitzky-Golay filtering process.
Definition: savgolfilter.cpp:588
pappso::pappso_double
double pappso_double
A type definition for doubles.
Definition: types.h:69
pappso::FilterSavitzkyGolay::convlv
char convlv(pappso_double data[], unsigned long n, pappso_double respns[], unsigned long m, int isign, pappso_double ans[])
Definition: savgolfilter.cpp:485
pappso::FilterSavitzkyGolay::~FilterSavitzkyGolay
virtual ~FilterSavitzkyGolay()
Definition: savgolfilter.cpp:77
savgolfilter.h
pappso
Definition: aa.cpp:38
pappso::FilterSavitzkyGolay::sgcoeff
char sgcoeff(pappso_double c[], int np, int nl, int nr, int ld, int m) const
Definition: savgolfilter.cpp:535
pappso::FilterSavitzkyGolay::lubksb
void lubksb(pappso_double **a, int n, int *indx, pappso_double b[]) const
Definition: savgolfilter.cpp:240
pappso::PeptideIonNter::a
@ a
pappso::FilterSavitzkyGolay::m_nR
int m_nR
number of data points on the right of the filtered point
Definition: savgolfilter.h:152
pappso::FilterSavitzkyGolay::ludcmp
void ludcmp(pappso_double **a, int n, int *indx, pappso_double *d) const
Definition: savgolfilter.cpp:271
pappso::FilterSavitzkyGolay::free_dvector
void free_dvector(pappso_double *v, long nl, long nh) const
Definition: savgolfilter.cpp:224
pappso::FilterSavitzkyGolay::free_ivector
void free_ivector(int *v, long nl, long nh) const
Definition: savgolfilter.cpp:217
pappso::FilterSavitzkyGolay::getParameters
SavGolParams getParameters() const
Definition: savgolfilter.cpp:140
pappso::FilterSavitzkyGolay::m_lD
int m_lD
Definition: savgolfilter.h:157
pappso::FilterSavitzkyGolay::toString
QString toString() const
Return a string with the textual representation of the configuration data.
Definition: savgolfilter.cpp:148
pappso::FilterSavitzkyGolay::realft
void realft(pappso_double data[], unsigned long n, int isign)
Definition: savgolfilter.cpp:433
pappso::Trace
A simple container of DataPoint instances.
Definition: trace.h:126
pappso::SavGolParams
Parameters for the Savitzky-Golay filter.
Definition: savgolfilter.h:48
pappso::FilterSavitzkyGolay::m_nL
int m_nL
number of data points on the left of the filtered point
Definition: savgolfilter.h:150
pappso::PeptideIonNter::c
@ c
pappso::FilterSavitzkyGolay::dvector
pappso_double * dvector(long nl, long nh) const
Definition: savgolfilter.cpp:175
pappso::FilterSavitzkyGolay::filter
Trace & filter(Trace &data_points) const override
Definition: savgolfilter.cpp:100
pappso::FilterSavitzkyGolay::m_convolveWithNr
bool m_convolveWithNr
set to false for best results
Definition: savgolfilter.h:160
pappso::FilterSavitzkyGolay::FilterSavitzkyGolay
FilterSavitzkyGolay(int nL, int nR, int m, int lD, bool convolveWithNr=false)
Definition: savgolfilter.cpp:54
pappso::FilterSavitzkyGolay::twofft
void twofft(pappso_double data1[], pappso_double data2[], pappso_double fft1[], pappso_double fft2[], unsigned long n)
Definition: savgolfilter.cpp:396
pappso::XicExtractMethod::sum
@ sum
sum of intensities
pappso::FilterSavitzkyGolay::operator=
FilterSavitzkyGolay & operator=(const FilterSavitzkyGolay &other)
Definition: savgolfilter.cpp:82
pappso::PeptideIonNter::b
@ b
pappso::FilterSavitzkyGolay::free_dmatrix
void free_dmatrix(pappso_double **m, long nrl, long nrh, long ncl, long nch) const
Definition: savgolfilter.cpp:231
pappso::FilterSavitzkyGolay::four1
void four1(pappso_double data[], unsigned long nn, int isign)
Definition: savgolfilter.cpp:342
pappso::FilterSavitzkyGolay
uses Savitsky-Golay filter on trace
Definition: savgolfilter.h:110
pappso::FilterSavitzkyGolay::m_m
int m_m
Definition: savgolfilter.h:154
pappso::FilterSavitzkyGolay::ivector
int * ivector(long nl, long nh) const
Definition: savgolfilter.cpp:162
SWAP
#define SWAP(a, b)
Definition: savgolfilter.cpp:48
pappso::FilterSavitzkyGolay::dmatrix
pappso_double ** dmatrix(long nrl, long nrh, long ncl, long nch) const
Definition: savgolfilter.cpp:191