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 #include "../../exception/exceptionnotrecognized.h"
41 
42 #include "savgolfilter.h"
43 
44 
45 namespace pappso
46 {
47 
48 
49 #define SWAP(a, b) \
50  tempr = (a); \
51  (a) = (b); \
52  (b) = tempr
53 
54 
56  int nL, int nR, int m, int lD, bool convolveWithNr)
57 {
58  m_params.nL = nL;
59  m_params.nR = nR;
60  m_params.m = m;
61  m_params.lD = lD;
62  m_params.convolveWithNr = convolveWithNr;
63 }
64 
65 
67 {
68  m_params.nL = sav_gol_params.nL;
69  m_params.nR = sav_gol_params.nR;
70  m_params.m = sav_gol_params.m;
71  m_params.lD = sav_gol_params.lD;
72  m_params.convolveWithNr = sav_gol_params.convolveWithNr;
73 }
74 
75 
77 {
78  // This function only copies the parameters, not the data.
79 
80  m_params.nL = other.m_params.nL;
81  m_params.nR = other.m_params.nR;
82  m_params.m = other.m_params.m;
83  m_params.lD = other.m_params.lD;
85 }
86 
87 
89 {
90 }
91 
94 {
95  if(&other == this)
96  return *this;
97 
98  // This function only copies the parameters, not the data.
99 
100  m_params.nL = other.m_params.nL;
101  m_params.nR = other.m_params.nR;
102  m_params.m = other.m_params.m;
103  m_params.lD = other.m_params.lD;
105 
106  return *this;
107 }
108 
109 
111 {
112  buildFilterFromString(parameters);
113 }
114 
115 
116 void
118 {
119  // Typical string: "Savitzky-Golay|15;15;4;0;false"
120  if(parameters.startsWith(QString("%1|").arg(name())))
121  {
122  QStringList params = parameters.split("|").back().split(";");
123 
124  m_params.nL = params.at(0).toInt();
125  m_params.nR = params.at(1).toInt();
126  m_params.m = params.at(2).toInt();
127  m_params.lD = params.at(3).toInt();
128  m_params.convolveWithNr = (params.at(4) == "true" ? true : false);
129  }
130  else
131  {
133  QString("Building of FilterSavitzkyGolay from string %1 failed")
134  .arg(parameters));
135  }
136 }
137 
138 
139 Trace &
141 {
142  // Initialize data:
143 
144  // We want the filter to stay constant so we create a local copy of the data.
145 
146  int data_point_count = data_points.size();
147 
148  pappso_double *x_data_p = dvector(1, data_point_count);
149  pappso_double *y_initial_data_p = dvector(1, data_point_count);
150  pappso_double *y_filtered_data_p = nullptr;
151 
153  y_filtered_data_p = dvector(1, 2 * data_point_count);
154  else
155  y_filtered_data_p = dvector(1, data_point_count);
156 
157  for(int iter = 0; iter < data_point_count; ++iter)
158  {
159  x_data_p[iter] = data_points.at(iter).x;
160  y_initial_data_p[iter] = data_points.at(iter).y;
161  }
162 
163  // Now run the filter.
164 
165  runFilter(y_initial_data_p, y_filtered_data_p, data_point_count);
166 
167  // Put back the modified y values into the trace.
168  auto iter_yf = y_filtered_data_p;
169  for(auto &data_point : data_points)
170  {
171  data_point.y = *iter_yf;
172  iter_yf++;
173  }
174 
175  return data_points;
176 }
177 
178 
181 {
182  return SavGolParams(
184 }
185 
186 
187 //! Return a string with the textual representation of the configuration data.
188 QString
190 {
191  return QString("%1|%2").arg(name()).arg(m_params.toString());
192 }
193 
194 
195 QString
197 {
198  return "Savitzky-Golay";
199 }
200 
201 
202 int *
203 FilterSavitzkyGolay::ivector(long nl, long nh) const
204 {
205  int *v;
206  v = (int *)malloc((size_t)((nh - nl + 2) * sizeof(int)));
207  if(!v)
208  {
209  qFatal("Error: Allocation failure.");
210  }
211  return v - nl + 1;
212 }
213 
214 
216 FilterSavitzkyGolay::dvector(long nl, long nh) const
217 {
218  pappso_double *v;
219  long k;
220  v = (pappso_double *)malloc((size_t)((nh - nl + 2) * sizeof(pappso_double)));
221  if(!v)
222  {
223  qFatal("Error: Allocation failure.");
224  }
225  for(k = nl; k <= nh; k++)
226  v[k] = 0.0;
227  return v - nl + 1;
228 }
229 
230 
231 pappso_double **
232 FilterSavitzkyGolay::dmatrix(long nrl, long nrh, long ncl, long nch) const
233 {
234  long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
235  pappso_double **m;
236  m = (pappso_double **)malloc((size_t)((nrow + 1) * sizeof(pappso_double *)));
237  if(!m)
238  {
239  qFatal("Error: Allocation failure.");
240  }
241  m += 1;
242  m -= nrl;
243  m[nrl] = (pappso_double *)malloc(
244  (size_t)((nrow * ncol + 1) * sizeof(pappso_double)));
245  if(!m[nrl])
246  {
247  qFatal("Error: Allocation failure.");
248  }
249  m[nrl] += 1;
250  m[nrl] -= ncl;
251  for(i = nrl + 1; i <= nrh; i++)
252  m[i] = m[i - 1] + ncol;
253  return m;
254 }
255 
256 
257 void
259  long nl,
260  [[maybe_unused]] long nh) const
261 {
262  free((char *)(v + nl - 1));
263 }
264 
265 
266 void
268  long nl,
269  [[maybe_unused]] long nh) const
270 {
271  free((char *)(v + nl - 1));
272 }
273 
274 
275 void
277  long nrl,
278  [[maybe_unused]] long nrh,
279  long ncl,
280  [[maybe_unused]] long nch) const
281 {
282  free((char *)(m[nrl] + ncl - 1));
283  free((char *)(m + nrl - 1));
284 }
285 
286 
287 void
289  int n,
290  int *indx,
291  pappso_double b[]) const
292 {
293  int i, ii = 0, ip, j;
295 
296  for(i = 1; i <= n; i++)
297  {
298  ip = indx[i];
299  sum = b[ip];
300  b[ip] = b[i];
301  if(ii)
302  for(j = ii; j <= i - 1; j++)
303  sum -= a[i][j] * b[j];
304  else if(sum)
305  ii = i;
306  b[i] = sum;
307  }
308  for(i = n; i >= 1; i--)
309  {
310  sum = b[i];
311  for(j = i + 1; j <= n; j++)
312  sum -= a[i][j] * b[j];
313  b[i] = sum / a[i][i];
314  }
315 }
316 
317 
318 void
320  int n,
321  int *indx,
322  pappso_double *d) const
323 {
324  int i, imax = 0, j, k;
325  pappso_double big, dum, sum, temp;
326  pappso_double *vv;
327 
328  vv = dvector(1, n);
329  *d = 1.0;
330  for(i = 1; i <= n; i++)
331  {
332  big = 0.0;
333  for(j = 1; j <= n; j++)
334  if((temp = fabs(a[i][j])) > big)
335  big = temp;
336  if(big == 0.0)
337  {
338  qFatal("Error: Singular matrix found in routine ludcmp().");
339  }
340  vv[i] = 1.0 / big;
341  }
342  for(j = 1; j <= n; j++)
343  {
344  for(i = 1; i < j; i++)
345  {
346  sum = a[i][j];
347  for(k = 1; k < i; k++)
348  sum -= a[i][k] * a[k][j];
349  a[i][j] = sum;
350  }
351  big = 0.0;
352  for(i = j; i <= n; i++)
353  {
354  sum = a[i][j];
355  for(k = 1; k < j; k++)
356  sum -= a[i][k] * a[k][j];
357  a[i][j] = sum;
358  if((dum = vv[i] * fabs(sum)) >= big)
359  {
360  big = dum;
361  imax = i;
362  }
363  }
364  if(j != imax)
365  {
366  for(k = 1; k <= n; k++)
367  {
368  dum = a[imax][k];
369  a[imax][k] = a[j][k];
370  a[j][k] = dum;
371  }
372  *d = -(*d);
373  vv[imax] = vv[j];
374  }
375  indx[j] = imax;
376  if(a[j][j] == 0.0)
377  a[j][j] = std::numeric_limits<pappso_double>::epsilon();
378  if(j != n)
379  {
380  dum = 1.0 / (a[j][j]);
381  for(i = j + 1; i <= n; i++)
382  a[i][j] *= dum;
383  }
384  }
385  free_dvector(vv, 1, n);
386 }
387 
388 
389 void
390 FilterSavitzkyGolay::four1(pappso_double data[], unsigned long nn, int isign)
391 {
392  unsigned long n, mmax, m, j, istep, i;
393  pappso_double wtemp, wr, wpr, wpi, wi, theta;
394  pappso_double tempr, tempi;
395 
396  n = nn << 1;
397  j = 1;
398  for(i = 1; i < n; i += 2)
399  {
400  if(j > i)
401  {
402  SWAP(data[j], data[i]);
403  SWAP(data[j + 1], data[i + 1]);
404  }
405  m = n >> 1;
406  while(m >= 2 && j > m)
407  {
408  j -= m;
409  m >>= 1;
410  }
411  j += m;
412  }
413  mmax = 2;
414  while(n > mmax)
415  {
416  istep = mmax << 1;
417  theta = isign * (6.28318530717959 / mmax);
418  wtemp = sin(0.5 * theta);
419  wpr = -2.0 * wtemp * wtemp;
420  wpi = sin(theta);
421  wr = 1.0;
422  wi = 0.0;
423  for(m = 1; m < mmax; m += 2)
424  {
425  for(i = m; i <= n; i += istep)
426  {
427  j = i + mmax;
428  tempr = wr * data[j] - wi * data[j + 1];
429  tempi = wr * data[j + 1] + wi * data[j];
430  data[j] = data[i] - tempr;
431  data[j + 1] = data[i + 1] - tempi;
432  data[i] += tempr;
433  data[i + 1] += tempi;
434  }
435  wr = (wtemp = wr) * wpr - wi * wpi + wr;
436  wi = wi * wpr + wtemp * wpi + wi;
437  }
438  mmax = istep;
439  }
440 }
441 
442 
443 void
445  pappso_double data2[],
446  pappso_double fft1[],
447  pappso_double fft2[],
448  unsigned long n)
449 {
450  unsigned long nn3, nn2, jj, j;
451  pappso_double rep, rem, aip, aim;
452 
453  nn3 = 1 + (nn2 = 2 + n + n);
454  for(j = 1, jj = 2; j <= n; j++, jj += 2)
455  {
456  fft1[jj - 1] = data1[j];
457  fft1[jj] = data2[j];
458  }
459  four1(fft1, n, 1);
460  fft2[1] = fft1[2];
461  fft1[2] = fft2[2] = 0.0;
462  for(j = 3; j <= n + 1; j += 2)
463  {
464  rep = 0.5 * (fft1[j] + fft1[nn2 - j]);
465  rem = 0.5 * (fft1[j] - fft1[nn2 - j]);
466  aip = 0.5 * (fft1[j + 1] + fft1[nn3 - j]);
467  aim = 0.5 * (fft1[j + 1] - fft1[nn3 - j]);
468  fft1[j] = rep;
469  fft1[j + 1] = aim;
470  fft1[nn2 - j] = rep;
471  fft1[nn3 - j] = -aim;
472  fft2[j] = aip;
473  fft2[j + 1] = -rem;
474  fft2[nn2 - j] = aip;
475  fft2[nn3 - j] = rem;
476  }
477 }
478 
479 
480 void
481 FilterSavitzkyGolay::realft(pappso_double data[], unsigned long n, int isign)
482 {
483  unsigned long i, i1, i2, i3, i4, np3;
484  pappso_double c1 = 0.5, c2, h1r, h1i, h2r, h2i;
485  pappso_double wr, wi, wpr, wpi, wtemp, theta;
486 
487  theta = 3.141592653589793 / (pappso_double)(n >> 1);
488  if(isign == 1)
489  {
490  c2 = -0.5;
491  four1(data, n >> 1, 1);
492  }
493  else
494  {
495  c2 = 0.5;
496  theta = -theta;
497  }
498  wtemp = sin(0.5 * theta);
499  wpr = -2.0 * wtemp * wtemp;
500  wpi = sin(theta);
501  wr = 1.0 + wpr;
502  wi = wpi;
503  np3 = n + 3;
504  for(i = 2; i <= (n >> 2); i++)
505  {
506  i4 = 1 + (i3 = np3 - (i2 = 1 + (i1 = i + i - 1)));
507  h1r = c1 * (data[i1] + data[i3]);
508  h1i = c1 * (data[i2] - data[i4]);
509  h2r = -c2 * (data[i2] + data[i4]);
510  h2i = c2 * (data[i1] - data[i3]);
511  data[i1] = h1r + wr * h2r - wi * h2i;
512  data[i2] = h1i + wr * h2i + wi * h2r;
513  data[i3] = h1r - wr * h2r + wi * h2i;
514  data[i4] = -h1i + wr * h2i + wi * h2r;
515  wr = (wtemp = wr) * wpr - wi * wpi + wr;
516  wi = wi * wpr + wtemp * wpi + wi;
517  }
518  if(isign == 1)
519  {
520  data[1] = (h1r = data[1]) + data[2];
521  data[2] = h1r - data[2];
522  }
523  else
524  {
525  data[1] = c1 * ((h1r = data[1]) + data[2]);
526  data[2] = c1 * (h1r - data[2]);
527  four1(data, n >> 1, -1);
528  }
529 }
530 
531 
532 char
534  unsigned long n,
535  pappso_double respns[],
536  unsigned long m,
537  int isign,
538  pappso_double ans[])
539 {
540  unsigned long i, no2;
541  pappso_double dum, mag2, *fft;
542 
543  fft = dvector(1, n << 1);
544  for(i = 1; i <= (m - 1) / 2; i++)
545  respns[n + 1 - i] = respns[m + 1 - i];
546  for(i = (m + 3) / 2; i <= n - (m - 1) / 2; i++)
547  respns[i] = 0.0;
548  twofft(data, respns, fft, ans, n);
549  no2 = n >> 1;
550  for(i = 2; i <= n + 2; i += 2)
551  {
552  if(isign == 1)
553  {
554  ans[i - 1] =
555  (fft[i - 1] * (dum = ans[i - 1]) - fft[i] * ans[i]) / no2;
556  ans[i] = (fft[i] * dum + fft[i - 1] * ans[i]) / no2;
557  }
558  else if(isign == -1)
559  {
560  if((mag2 = ans[i - 1] * ans[i - 1] + ans[i] * ans[i]) == 0.0)
561  {
562  qDebug("Attempt of deconvolving at zero response in convlv().");
563  return (1);
564  }
565  ans[i - 1] =
566  (fft[i - 1] * (dum = ans[i - 1]) + fft[i] * ans[i]) / mag2 / no2;
567  ans[i] = (fft[i] * dum - fft[i - 1] * ans[i]) / mag2 / no2;
568  }
569  else
570  {
571  qDebug("No meaning for isign in convlv().");
572  return (1);
573  }
574  }
575  ans[2] = ans[n + 1];
576  realft(ans, n, -1);
577  free_dvector(fft, 1, n << 1);
578  return (0);
579 }
580 
581 
582 char
584  pappso_double c[], int np, int nl, int nr, int ld, int m) const
585 {
586  int imj, ipj, j, k, kk, mm, *indx;
587  pappso_double d, fac, sum, **a, *b;
588 
589  if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m)
590  {
591  qDebug("Inconsistent arguments detected in routine sgcoeff.");
592  return (1);
593  }
594  indx = ivector(1, m + 1);
595  a = dmatrix(1, m + 1, 1, m + 1);
596  b = dvector(1, m + 1);
597  for(ipj = 0; ipj <= (m << 1); ipj++)
598  {
599  sum = (ipj ? 0.0 : 1.0);
600  for(k = 1; k <= nr; k++)
601  sum += pow((pappso_double)k, (pappso_double)ipj);
602  for(k = 1; k <= nl; k++)
603  sum += pow((pappso_double)-k, (pappso_double)ipj);
604  mm = (ipj < 2 * m - ipj ? ipj : 2 * m - ipj);
605  for(imj = -mm; imj <= mm; imj += 2)
606  a[1 + (ipj + imj) / 2][1 + (ipj - imj) / 2] = sum;
607  }
608  ludcmp(a, m + 1, indx, &d);
609  for(j = 1; j <= m + 1; j++)
610  b[j] = 0.0;
611  b[ld + 1] = 1.0;
612  lubksb(a, m + 1, indx, b);
613  for(kk = 1; kk <= np; kk++)
614  c[kk] = 0.0;
615  for(k = -nl; k <= nr; k++)
616  {
617  sum = b[1];
618  fac = 1.0;
619  for(mm = 1; mm <= m; mm++)
620  sum += b[mm + 1] * (fac *= k);
621  kk = ((np - k) % np) + 1;
622  c[kk] = sum;
623  }
624  free_dvector(b, 1, m + 1);
625  free_dmatrix(a, 1, m + 1, 1, m + 1);
626  free_ivector(indx, 1, m + 1);
627  return (0);
628 }
629 
630 
631 //! Perform the Savitzky-Golay filtering process.
632 /*
633  The results are in the \c y_filtered_data_p C array of pappso_double
634  values.
635  */
636 char
638  double *y_filtered_data_p,
639  int data_point_count) const
640 {
641  int np = m_params.nL + 1 + m_params.nR;
642  pappso_double *c;
643  char retval;
644 
645 #if CONVOLVE_WITH_NR_CONVLV
646  c = dvector(1, data_point_count);
647  retval = sgcoeff(c, np, m_nL, m_nR, m_lD, m_m);
648  if(retval == 0)
649  convlv(y_data_p, data_point_count, c, np, 1, y_filtered_data_p);
650  free_dvector(c, 1, data_point_count);
651 #else
652  int j;
653  long int k;
654  c = dvector(1, m_params.nL + m_params.nR + 1);
655  retval = sgcoeff(c, np, m_params.nL, m_params.nR, m_params.lD, m_params.m);
656  if(retval == 0)
657  {
658  qDebug() << __FILE__ << __LINE__ << __FUNCTION__ << "()"
659  << "retval is 0";
660 
661  for(k = 1; k <= m_params.nL; k++)
662  {
663  for(y_filtered_data_p[k] = 0.0, j = -m_params.nL; j <= m_params.nR;
664  j++)
665  {
666  if(k + j >= 1)
667  {
668  y_filtered_data_p[k] +=
669  c[(j >= 0 ? j + 1 : m_params.nR + m_params.nL + 2 + j)] *
670  y_data_p[k + j];
671  }
672  }
673  }
674  for(k = m_params.nL + 1; k <= data_point_count - m_params.nR; k++)
675  {
676  for(y_filtered_data_p[k] = 0.0, j = -m_params.nL; j <= m_params.nR;
677  j++)
678  {
679  y_filtered_data_p[k] +=
680  c[(j >= 0 ? j + 1 : m_params.nR + m_params.nL + 2 + j)] *
681  y_data_p[k + j];
682  }
683  }
684  for(k = data_point_count - m_params.nR + 1; k <= data_point_count; k++)
685  {
686  for(y_filtered_data_p[k] = 0.0, j = -m_params.nL; j <= m_params.nR;
687  j++)
688  {
689  if(k + j <= data_point_count)
690  {
691  y_filtered_data_p[k] +=
692  c[(j >= 0 ? j + 1 : m_params.nR + m_params.nL + 2 + j)] *
693  y_data_p[k + j];
694  }
695  }
696  }
697  }
698 
699  free_dvector(c, 1, m_params.nR + m_params.nL + 1);
700 #endif
701 
702  return (retval);
703 }
704 
705 
706 } // namespace pappso
excetion to use when an item type is not recognized
uses Savitsky-Golay filter on trace
Definition: savgolfilter.h:129
void free_dmatrix(pappso_double **m, long nrl, long nrh, long ncl, long nch) const
pappso_double * dvector(long nl, long nh) const
void buildFilterFromString(const QString &strBuildParams) override
build this filter using a string
void free_dvector(pappso_double *v, long nl, long nh) const
QString name() const override
pappso_double ** dmatrix(long nrl, long nrh, long ncl, long nch) const
QString toString() const override
Return a string with the textual representation of the configuration data.
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)
void free_ivector(int *v, long nl, long nh) const
char convlv(pappso_double data[], unsigned long n, pappso_double respns[], unsigned long m, int isign, pappso_double ans[])
void lubksb(pappso_double **a, int n, int *indx, pappso_double b[]) const
int * ivector(long nl, long nh) const
FilterSavitzkyGolay(int nL, int nR, int m, int lD, bool convolveWithNr=false)
void ludcmp(pappso_double **a, int n, int *indx, pappso_double *d) const
void four1(pappso_double data[], unsigned long nn, int isign)
FilterSavitzkyGolay & operator=(const FilterSavitzkyGolay &other)
Trace & filter(Trace &data_points) const override
char sgcoeff(pappso_double c[], int np, int nl, int nr, int ld, int m) const
char runFilter(double *y_data_p, double *y_filtered_data_p, int data_point_count) const
Perform the Savitzky-Golay filtering process.
SavGolParams getParameters() const
A simple container of DataPoint instances.
Definition: trace.h:148
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
double pappso_double
A type definition for doubles.
Definition: types.h:49
@ sum
sum of intensities
#define SWAP(a, b)
Parameters for the Savitzky-Golay filter.
Definition: savgolfilter.h:50
QString toString() const
Definition: savgolfilter.h:107
int nR
number of data points on the right of the filtered point
Definition: savgolfilter.h:53
int nL
number of data points on the left of the filtered point
Definition: savgolfilter.h:51
bool convolveWithNr
set to false for best results
Definition: savgolfilter.h:61