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
218  long nl,
219  [[maybe_unused]] long nh) const
220 {
221  free((char *)(v + nl - 1));
222 }
223 
224 
225 void
227  long nl,
228  [[maybe_unused]] long nh) const
229 {
230  free((char *)(v + nl - 1));
231 }
232 
233 
234 void
236  long nrl,
237  [[maybe_unused]] long nrh,
238  long ncl,
239  [[maybe_unused]] long nch) const
240 {
241  free((char *)(m[nrl] + ncl - 1));
242  free((char *)(m + nrl - 1));
243 }
244 
245 
246 void
248  int n,
249  int *indx,
250  pappso_double b[]) const
251 {
252  int i, ii = 0, ip, j;
254 
255  for(i = 1; i <= n; i++)
256  {
257  ip = indx[i];
258  sum = b[ip];
259  b[ip] = b[i];
260  if(ii)
261  for(j = ii; j <= i - 1; j++)
262  sum -= a[i][j] * b[j];
263  else if(sum)
264  ii = i;
265  b[i] = sum;
266  }
267  for(i = n; i >= 1; i--)
268  {
269  sum = b[i];
270  for(j = i + 1; j <= n; j++)
271  sum -= a[i][j] * b[j];
272  b[i] = sum / a[i][i];
273  }
274 }
275 
276 
277 void
279  int n,
280  int *indx,
281  pappso_double *d) const
282 {
283  int i, imax = 0, j, k;
284  pappso_double big, dum, sum, temp;
285  pappso_double *vv;
286 
287  vv = dvector(1, n);
288  *d = 1.0;
289  for(i = 1; i <= n; i++)
290  {
291  big = 0.0;
292  for(j = 1; j <= n; j++)
293  if((temp = fabs(a[i][j])) > big)
294  big = temp;
295  if(big == 0.0)
296  {
297  qFatal("Error: Singular matrix found in routine ludcmp().");
298  }
299  vv[i] = 1.0 / big;
300  }
301  for(j = 1; j <= n; j++)
302  {
303  for(i = 1; i < j; i++)
304  {
305  sum = a[i][j];
306  for(k = 1; k < i; k++)
307  sum -= a[i][k] * a[k][j];
308  a[i][j] = sum;
309  }
310  big = 0.0;
311  for(i = j; i <= n; i++)
312  {
313  sum = a[i][j];
314  for(k = 1; k < j; k++)
315  sum -= a[i][k] * a[k][j];
316  a[i][j] = sum;
317  if((dum = vv[i] * fabs(sum)) >= big)
318  {
319  big = dum;
320  imax = i;
321  }
322  }
323  if(j != imax)
324  {
325  for(k = 1; k <= n; k++)
326  {
327  dum = a[imax][k];
328  a[imax][k] = a[j][k];
329  a[j][k] = dum;
330  }
331  *d = -(*d);
332  vv[imax] = vv[j];
333  }
334  indx[j] = imax;
335  if(a[j][j] == 0.0)
336  a[j][j] = std::numeric_limits<pappso_double>::epsilon();
337  if(j != n)
338  {
339  dum = 1.0 / (a[j][j]);
340  for(i = j + 1; i <= n; i++)
341  a[i][j] *= dum;
342  }
343  }
344  free_dvector(vv, 1, n);
345 }
346 
347 
348 void
349 FilterSavitzkyGolay::four1(pappso_double data[], unsigned long nn, int isign)
350 {
351  unsigned long n, mmax, m, j, istep, i;
352  pappso_double wtemp, wr, wpr, wpi, wi, theta;
353  pappso_double tempr, tempi;
354 
355  n = nn << 1;
356  j = 1;
357  for(i = 1; i < n; i += 2)
358  {
359  if(j > i)
360  {
361  SWAP(data[j], data[i]);
362  SWAP(data[j + 1], data[i + 1]);
363  }
364  m = n >> 1;
365  while(m >= 2 && j > m)
366  {
367  j -= m;
368  m >>= 1;
369  }
370  j += m;
371  }
372  mmax = 2;
373  while(n > mmax)
374  {
375  istep = mmax << 1;
376  theta = isign * (6.28318530717959 / mmax);
377  wtemp = sin(0.5 * theta);
378  wpr = -2.0 * wtemp * wtemp;
379  wpi = sin(theta);
380  wr = 1.0;
381  wi = 0.0;
382  for(m = 1; m < mmax; m += 2)
383  {
384  for(i = m; i <= n; i += istep)
385  {
386  j = i + mmax;
387  tempr = wr * data[j] - wi * data[j + 1];
388  tempi = wr * data[j + 1] + wi * data[j];
389  data[j] = data[i] - tempr;
390  data[j + 1] = data[i + 1] - tempi;
391  data[i] += tempr;
392  data[i + 1] += tempi;
393  }
394  wr = (wtemp = wr) * wpr - wi * wpi + wr;
395  wi = wi * wpr + wtemp * wpi + wi;
396  }
397  mmax = istep;
398  }
399 }
400 
401 
402 void
404  pappso_double data2[],
405  pappso_double fft1[],
406  pappso_double fft2[],
407  unsigned long n)
408 {
409  unsigned long nn3, nn2, jj, j;
410  pappso_double rep, rem, aip, aim;
411 
412  nn3 = 1 + (nn2 = 2 + n + n);
413  for(j = 1, jj = 2; j <= n; j++, jj += 2)
414  {
415  fft1[jj - 1] = data1[j];
416  fft1[jj] = data2[j];
417  }
418  four1(fft1, n, 1);
419  fft2[1] = fft1[2];
420  fft1[2] = fft2[2] = 0.0;
421  for(j = 3; j <= n + 1; j += 2)
422  {
423  rep = 0.5 * (fft1[j] + fft1[nn2 - j]);
424  rem = 0.5 * (fft1[j] - fft1[nn2 - j]);
425  aip = 0.5 * (fft1[j + 1] + fft1[nn3 - j]);
426  aim = 0.5 * (fft1[j + 1] - fft1[nn3 - j]);
427  fft1[j] = rep;
428  fft1[j + 1] = aim;
429  fft1[nn2 - j] = rep;
430  fft1[nn3 - j] = -aim;
431  fft2[j] = aip;
432  fft2[j + 1] = -rem;
433  fft2[nn2 - j] = aip;
434  fft2[nn3 - j] = rem;
435  }
436 }
437 
438 
439 void
440 FilterSavitzkyGolay::realft(pappso_double data[], unsigned long n, int isign)
441 {
442  unsigned long i, i1, i2, i3, i4, np3;
443  pappso_double c1 = 0.5, c2, h1r, h1i, h2r, h2i;
444  pappso_double wr, wi, wpr, wpi, wtemp, theta;
445 
446  theta = 3.141592653589793 / (pappso_double)(n >> 1);
447  if(isign == 1)
448  {
449  c2 = -0.5;
450  four1(data, n >> 1, 1);
451  }
452  else
453  {
454  c2 = 0.5;
455  theta = -theta;
456  }
457  wtemp = sin(0.5 * theta);
458  wpr = -2.0 * wtemp * wtemp;
459  wpi = sin(theta);
460  wr = 1.0 + wpr;
461  wi = wpi;
462  np3 = n + 3;
463  for(i = 2; i <= (n >> 2); i++)
464  {
465  i4 = 1 + (i3 = np3 - (i2 = 1 + (i1 = i + i - 1)));
466  h1r = c1 * (data[i1] + data[i3]);
467  h1i = c1 * (data[i2] - data[i4]);
468  h2r = -c2 * (data[i2] + data[i4]);
469  h2i = c2 * (data[i1] - data[i3]);
470  data[i1] = h1r + wr * h2r - wi * h2i;
471  data[i2] = h1i + wr * h2i + wi * h2r;
472  data[i3] = h1r - wr * h2r + wi * h2i;
473  data[i4] = -h1i + wr * h2i + wi * h2r;
474  wr = (wtemp = wr) * wpr - wi * wpi + wr;
475  wi = wi * wpr + wtemp * wpi + wi;
476  }
477  if(isign == 1)
478  {
479  data[1] = (h1r = data[1]) + data[2];
480  data[2] = h1r - data[2];
481  }
482  else
483  {
484  data[1] = c1 * ((h1r = data[1]) + data[2]);
485  data[2] = c1 * (h1r - data[2]);
486  four1(data, n >> 1, -1);
487  }
488 }
489 
490 
491 char
493  unsigned long n,
494  pappso_double respns[],
495  unsigned long m,
496  int isign,
497  pappso_double ans[])
498 {
499  unsigned long i, no2;
500  pappso_double dum, mag2, *fft;
501 
502  fft = dvector(1, n << 1);
503  for(i = 1; i <= (m - 1) / 2; i++)
504  respns[n + 1 - i] = respns[m + 1 - i];
505  for(i = (m + 3) / 2; i <= n - (m - 1) / 2; i++)
506  respns[i] = 0.0;
507  twofft(data, respns, fft, ans, n);
508  no2 = n >> 1;
509  for(i = 2; i <= n + 2; i += 2)
510  {
511  if(isign == 1)
512  {
513  ans[i - 1] =
514  (fft[i - 1] * (dum = ans[i - 1]) - fft[i] * ans[i]) / no2;
515  ans[i] = (fft[i] * dum + fft[i - 1] * ans[i]) / no2;
516  }
517  else if(isign == -1)
518  {
519  if((mag2 = ans[i - 1] * ans[i - 1] + ans[i] * ans[i]) == 0.0)
520  {
521  qDebug("Attempt of deconvolving at zero response in convlv().");
522  return (1);
523  }
524  ans[i - 1] =
525  (fft[i - 1] * (dum = ans[i - 1]) + fft[i] * ans[i]) / mag2 / no2;
526  ans[i] = (fft[i] * dum - fft[i - 1] * ans[i]) / mag2 / no2;
527  }
528  else
529  {
530  qDebug("No meaning for isign in convlv().");
531  return (1);
532  }
533  }
534  ans[2] = ans[n + 1];
535  realft(ans, n, -1);
536  free_dvector(fft, 1, n << 1);
537  return (0);
538 }
539 
540 
541 char
543  pappso_double c[], int np, int nl, int nr, int ld, int m) const
544 {
545  int imj, ipj, j, k, kk, mm, *indx;
546  pappso_double d, fac, sum, **a, *b;
547 
548  if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m)
549  {
550  qDebug("Inconsistent arguments detected in routine sgcoeff.");
551  return (1);
552  }
553  indx = ivector(1, m + 1);
554  a = dmatrix(1, m + 1, 1, m + 1);
555  b = dvector(1, m + 1);
556  for(ipj = 0; ipj <= (m << 1); ipj++)
557  {
558  sum = (ipj ? 0.0 : 1.0);
559  for(k = 1; k <= nr; k++)
560  sum += pow((pappso_double)k, (pappso_double)ipj);
561  for(k = 1; k <= nl; k++)
562  sum += pow((pappso_double)-k, (pappso_double)ipj);
563  mm = (ipj < 2 * m - ipj ? ipj : 2 * m - ipj);
564  for(imj = -mm; imj <= mm; imj += 2)
565  a[1 + (ipj + imj) / 2][1 + (ipj - imj) / 2] = sum;
566  }
567  ludcmp(a, m + 1, indx, &d);
568  for(j = 1; j <= m + 1; j++)
569  b[j] = 0.0;
570  b[ld + 1] = 1.0;
571  lubksb(a, m + 1, indx, b);
572  for(kk = 1; kk <= np; kk++)
573  c[kk] = 0.0;
574  for(k = -nl; k <= nr; k++)
575  {
576  sum = b[1];
577  fac = 1.0;
578  for(mm = 1; mm <= m; mm++)
579  sum += b[mm + 1] * (fac *= k);
580  kk = ((np - k) % np) + 1;
581  c[kk] = sum;
582  }
583  free_dvector(b, 1, m + 1);
584  free_dmatrix(a, 1, m + 1, 1, m + 1);
585  free_ivector(indx, 1, m + 1);
586  return (0);
587 }
588 
589 
590 //! Perform the Savitzky-Golay filtering process.
591 /*
592  The results are in the \c y_filtered_data_p C array of pappso_double values.
593  */
594 char
596  double *y_filtered_data_p,
597  int data_point_count) const
598 {
599  int np = m_nL + 1 + m_nR;
600  pappso_double *c;
601  char retval;
602 
603 #if CONVOLVE_WITH_NR_CONVLV
604  c = dvector(1, data_point_count);
605  retval = sgcoeff(c, np, m_nL, m_nR, m_lD, m_m);
606  if(retval == 0)
607  convlv(y_data_p, data_point_count, c, np, 1, y_filtered_data_p);
608  free_dvector(c, 1, data_point_count);
609 #else
610  int j;
611  long int k;
612  c = dvector(1, m_nL + m_nR + 1);
613  retval = sgcoeff(c, np, m_nL, m_nR, m_lD, m_m);
614  if(retval == 0)
615  {
616  qDebug() << __FILE__ << __LINE__ << __FUNCTION__ << "()"
617  << "retval is 0";
618 
619  for(k = 1; k <= m_nL; k++)
620  {
621  for(y_filtered_data_p[k] = 0.0, j = -m_nL; j <= m_nR; j++)
622  {
623  if(k + j >= 1)
624  {
625  y_filtered_data_p[k] +=
626  c[(j >= 0 ? j + 1 : m_nR + m_nL + 2 + j)] * y_data_p[k + j];
627  }
628  }
629  }
630  for(k = m_nL + 1; k <= data_point_count - m_nR; k++)
631  {
632  for(y_filtered_data_p[k] = 0.0, j = -m_nL; j <= m_nR; j++)
633  {
634  y_filtered_data_p[k] +=
635  c[(j >= 0 ? j + 1 : m_nR + m_nL + 2 + j)] * y_data_p[k + j];
636  }
637  }
638  for(k = data_point_count - m_nR + 1; k <= data_point_count; k++)
639  {
640  for(y_filtered_data_p[k] = 0.0, j = -m_nL; j <= m_nR; j++)
641  {
642  if(k + j <= data_point_count)
643  {
644  y_filtered_data_p[k] +=
645  c[(j >= 0 ? j + 1 : m_nR + m_nL + 2 + j)] * y_data_p[k + j];
646  }
647  }
648  }
649  }
650 
651  free_dvector(c, 1, m_nR + m_nL + 1);
652 #endif
653 
654  return (retval);
655 }
656 
657 
658 } // 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:595
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:492
pappso::FilterSavitzkyGolay::~FilterSavitzkyGolay
virtual ~FilterSavitzkyGolay()
Definition: savgolfilter.cpp:77
savgolfilter.h
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
pappso::FilterSavitzkyGolay::sgcoeff
char sgcoeff(pappso_double c[], int np, int nl, int nr, int ld, int m) const
Definition: savgolfilter.cpp:542
pappso::FilterSavitzkyGolay::lubksb
void lubksb(pappso_double **a, int n, int *indx, pappso_double b[]) const
Definition: savgolfilter.cpp:247
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:278
pappso::FilterSavitzkyGolay::free_dvector
void free_dvector(pappso_double *v, long nl, long nh) const
Definition: savgolfilter.cpp:226
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:440
pappso::Trace
A simple container of DataPoint instances.
Definition: trace.h:132
pappso::SavGolParams
Parameters for the Savitzky-Golay filter.
Definition: savgolfilter.h:49
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:403
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:235
pappso::FilterSavitzkyGolay::four1
void four1(pappso_double data[], unsigned long nn, int isign)
Definition: savgolfilter.cpp:349
pappso::FilterSavitzkyGolay
uses Savitsky-Golay filter on trace
Definition: savgolfilter.h:111
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