libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
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
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
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
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:50
@ sum
sum of intensities
#define SWAP(a, b)
Parameters for the Savitzky-Golay filter.
QString toString() const
int nR
number of data points on the right of the filtered point
int nL
number of data points on the left of the filtered point
bool convolveWithNr
set to false for best results