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
45namespace 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
116void
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
139Trace &
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.
188QString
190{
191 return QString("%1|%2").arg(name()).arg(m_params.toString());
192}
193
194
195QString
197{
198 return "Savitzky-Golay";
199}
200
201
202int *
203FilterSavitzkyGolay::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
216FilterSavitzkyGolay::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
232FilterSavitzkyGolay::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
257void
259 long nl,
260 [[maybe_unused]] long nh) const
261{
262 free((char *)(v + nl - 1));
263}
264
265
266void
268 long nl,
269 [[maybe_unused]] long nh) const
270{
271 free((char *)(v + nl - 1));
272}
273
274
275void
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
287void
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
318void
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
389void
390FilterSavitzkyGolay::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
443void
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
480void
481FilterSavitzkyGolay::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
532char
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
582char
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 */
636char
638 double *y_filtered_data_p,
639 int data_point_count) const
640{
641 int np = m_params.nL + 1 + m_params.nR;
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