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