40 #include "../../exception/exceptionnotrecognized.h"
56 int nL,
int nR,
int m,
int lD,
bool convolveWithNr)
120 if(parameters.startsWith(QString(
"%1|").arg(
name())))
122 QStringList params = parameters.split(
"|").back().split(
";");
133 QString(
"Building of FilterSavitzkyGolay from string %1 failed")
146 int data_point_count = data_points.size();
153 y_filtered_data_p =
dvector(1, 2 * data_point_count);
155 y_filtered_data_p =
dvector(1, data_point_count);
157 for(
int iter = 0; iter < data_point_count; ++iter)
159 x_data_p[iter] = data_points.at(iter).x;
160 y_initial_data_p[iter] = data_points.at(iter).y;
165 runFilter(y_initial_data_p, y_filtered_data_p, data_point_count);
168 auto iter_yf = y_filtered_data_p;
169 for(
auto &data_point : data_points)
171 data_point.y = *iter_yf;
198 return "Savitzky-Golay";
206 v = (
int *)malloc((
size_t)((nh - nl + 2) *
sizeof(
int)));
209 qFatal(
"Error: Allocation failure.");
223 qFatal(
"Error: Allocation failure.");
225 for(k = nl; k <= nh; k++)
234 long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
239 qFatal(
"Error: Allocation failure.");
247 qFatal(
"Error: Allocation failure.");
251 for(i = nrl + 1; i <= nrh; i++)
252 m[i] = m[i - 1] + ncol;
260 [[maybe_unused]]
long nh)
const
262 free((
char *)(v + nl - 1));
269 [[maybe_unused]]
long nh)
const
271 free((
char *)(v + nl - 1));
278 [[maybe_unused]]
long nrh,
280 [[maybe_unused]]
long nch)
const
282 free((
char *)(m[nrl] + ncl - 1));
283 free((
char *)(m + nrl - 1));
293 int i, ii = 0, ip, j;
296 for(i = 1; i <= n; i++)
302 for(j = ii; j <= i - 1; j++)
303 sum -=
a[i][j] *
b[j];
308 for(i = n; i >= 1; i--)
311 for(j = i + 1; j <= n; j++)
312 sum -=
a[i][j] *
b[j];
313 b[i] =
sum /
a[i][i];
324 int i, imax = 0, j, k;
330 for(i = 1; i <= n; i++)
333 for(j = 1; j <= n; j++)
334 if((temp = fabs(
a[i][j])) > big)
338 qFatal(
"Error: Singular matrix found in routine ludcmp().");
342 for(j = 1; j <= n; j++)
344 for(i = 1; i < j; i++)
347 for(k = 1; k < i; k++)
348 sum -=
a[i][k] *
a[k][j];
352 for(i = j; i <= n; i++)
355 for(k = 1; k < j; k++)
356 sum -=
a[i][k] *
a[k][j];
358 if((dum = vv[i] * fabs(
sum)) >= big)
366 for(k = 1; k <= n; k++)
369 a[imax][k] =
a[j][k];
377 a[j][j] = std::numeric_limits<pappso_double>::epsilon();
380 dum = 1.0 / (
a[j][j]);
381 for(i = j + 1; i <= n; i++)
392 unsigned long n, mmax, m, j, istep, i;
398 for(i = 1; i < n; i += 2)
402 SWAP(data[j], data[i]);
403 SWAP(data[j + 1], data[i + 1]);
406 while(m >= 2 && j > m)
417 theta = isign * (6.28318530717959 / mmax);
418 wtemp = sin(0.5 * theta);
419 wpr = -2.0 * wtemp * wtemp;
423 for(m = 1; m < mmax; m += 2)
425 for(i = m; i <= n; i += istep)
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;
433 data[i + 1] += tempi;
435 wr = (wtemp = wr) * wpr - wi * wpi + wr;
436 wi = wi * wpr + wtemp * wpi + wi;
450 unsigned long nn3, nn2, jj, j;
453 nn3 = 1 + (nn2 = 2 + n + n);
454 for(j = 1, jj = 2; j <= n; j++, jj += 2)
456 fft1[jj - 1] = data1[j];
461 fft1[2] = fft2[2] = 0.0;
462 for(j = 3; j <= n + 1; j += 2)
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]);
471 fft1[nn3 - j] = -aim;
483 unsigned long i, i1, i2, i3, i4, np3;
491 four1(data, n >> 1, 1);
498 wtemp = sin(0.5 * theta);
499 wpr = -2.0 * wtemp * wtemp;
504 for(i = 2; i <= (n >> 2); i++)
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;
520 data[1] = (h1r = data[1]) + data[2];
521 data[2] = h1r - data[2];
525 data[1] = c1 * ((h1r = data[1]) + data[2]);
526 data[2] = c1 * (h1r - data[2]);
527 four1(data, n >> 1, -1);
540 unsigned long i, no2;
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++)
548 twofft(data, respns, fft, ans, n);
550 for(i = 2; i <= n + 2; i += 2)
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;
560 if((mag2 = ans[i - 1] * ans[i - 1] + ans[i] * ans[i]) == 0.0)
562 qDebug(
"Attempt of deconvolving at zero response in convlv().");
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;
571 qDebug(
"No meaning for isign in convlv().");
584 pappso_double c[],
int np,
int nl,
int nr,
int ld,
int m)
const
586 int imj, ipj, j, k, kk, mm, *indx;
589 if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m)
591 qDebug(
"Inconsistent arguments detected in routine sgcoeff.");
597 for(ipj = 0; ipj <= (m << 1); ipj++)
599 sum = (ipj ? 0.0 : 1.0);
600 for(k = 1; k <= nr; k++)
602 for(k = 1; k <= nl; k++)
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;
609 for(j = 1; j <= m + 1; j++)
613 for(kk = 1; kk <= np; kk++)
615 for(k = -nl; k <= nr; k++)
619 for(mm = 1; mm <= m; mm++)
620 sum +=
b[mm + 1] * (fac *= k);
621 kk = ((np - k) % np) + 1;
638 double *y_filtered_data_p,
639 int data_point_count)
const
645 #if CONVOLVE_WITH_NR_CONVLV
647 retval =
sgcoeff(
c, np, m_nL, m_nR, m_lD, m_m);
649 convlv(y_data_p, data_point_count,
c, np, 1, y_filtered_data_p);
658 qDebug() << __FILE__ << __LINE__ << __FUNCTION__ <<
"()"
668 y_filtered_data_p[k] +=
679 y_filtered_data_p[k] +=
684 for(k = data_point_count -
m_params.
nR + 1; k <= data_point_count; k++)
689 if(k + j <= data_point_count)
691 y_filtered_data_p[k] +=
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.
virtual ~FilterSavitzkyGolay()
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.
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
double pappso_double
A type definition for doubles.
Parameters for the Savitzky-Golay filter.
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