60 int nL,
int nR,
int m,
int lD,
bool convolveWithNr)
124 if(parameters.startsWith(QString(
"%1|").arg(
name())))
126 QStringList params = parameters.split(
"|").back().split(
";");
137 QString(
"Building of FilterSavitzkyGolay from string %1 failed")
150 int data_point_count = data_points.size();
157 y_filtered_data_p =
dvector(1, 2 * data_point_count);
159 y_filtered_data_p =
dvector(1, data_point_count);
161 for(
int iter = 0; iter < data_point_count; ++iter)
163 x_data_p[iter] = data_points.at(iter).x;
164 y_initial_data_p[iter] = data_points.at(iter).y;
169 runFilter(y_initial_data_p, y_filtered_data_p, data_point_count);
172 auto iter_yf = y_filtered_data_p;
173 for(
auto &data_point : data_points)
175 data_point.y = *iter_yf;
202 return "Savitzky-Golay";
210 v = (
int *)malloc((
size_t)((nh - nl + 2) *
sizeof(
int)));
213 qFatal(
"Error: Allocation failure.");
227 qFatal(
"Error: Allocation failure.");
229 for(k = nl; k <= nh; k++)
238 long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
243 qFatal(
"Error: Allocation failure.");
251 qFatal(
"Error: Allocation failure.");
255 for(i = nrl + 1; i <= nrh; i++)
256 m[i] = m[i - 1] + ncol;
264 [[maybe_unused]]
long nh)
const
266 free((
char *)(v + nl - 1));
273 [[maybe_unused]]
long nh)
const
275 free((
char *)(v + nl - 1));
282 [[maybe_unused]]
long nrh,
284 [[maybe_unused]]
long nch)
const
286 free((
char *)(m[nrl] + ncl - 1));
287 free((
char *)(m + nrl - 1));
297 int i, ii = 0, ip, j;
300 for(i = 1; i <= n; i++)
306 for(j = ii; j <= i - 1; j++)
307 sum -=
a[i][j] *
b[j];
312 for(i = n; i >= 1; i--)
315 for(j = i + 1; j <= n; j++)
316 sum -=
a[i][j] *
b[j];
317 b[i] =
sum /
a[i][i];
328 int i, imax = 0, j, k;
334 for(i = 1; i <= n; i++)
337 for(j = 1; j <= n; j++)
338 if((temp = fabs(
a[i][j])) > big)
342 qFatal(
"Error: Singular matrix found in routine ludcmp().");
346 for(j = 1; j <= n; j++)
348 for(i = 1; i < j; i++)
351 for(k = 1; k < i; k++)
352 sum -=
a[i][k] *
a[k][j];
356 for(i = j; i <= n; i++)
359 for(k = 1; k < j; k++)
360 sum -=
a[i][k] *
a[k][j];
362 if((dum = vv[i] * fabs(
sum)) >= big)
370 for(k = 1; k <= n; k++)
373 a[imax][k] =
a[j][k];
381 a[j][j] = std::numeric_limits<pappso_double>::epsilon();
384 dum = 1.0 / (
a[j][j]);
385 for(i = j + 1; i <= n; i++)
396 unsigned long n, mmax, m, j, istep, i;
402 for(i = 1; i < n; i += 2)
406 SWAP(data[j], data[i]);
407 SWAP(data[j + 1], data[i + 1]);
410 while(m >= 2 && j > m)
421 theta = isign * (6.28318530717959 / mmax);
422 wtemp = sin(0.5 * theta);
423 wpr = -2.0 * wtemp * wtemp;
427 for(m = 1; m < mmax; m += 2)
429 for(i = m; i <= n; i += istep)
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;
437 data[i + 1] += tempi;
439 wr = (wtemp = wr) * wpr - wi * wpi + wr;
440 wi = wi * wpr + wtemp * wpi + wi;
454 unsigned long nn3, nn2, jj, j;
457 nn3 = 1 + (nn2 = 2 + n + n);
458 for(j = 1, jj = 2; j <= n; j++, jj += 2)
460 fft1[jj - 1] = data1[j];
465 fft1[2] = fft2[2] = 0.0;
466 for(j = 3; j <= n + 1; j += 2)
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]);
475 fft1[nn3 - j] = -aim;
487 unsigned long i, i1, i2, i3, i4, np3;
495 four1(data, n >> 1, 1);
502 wtemp = sin(0.5 * theta);
503 wpr = -2.0 * wtemp * wtemp;
508 for(i = 2; i <= (n >> 2); i++)
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;
524 data[1] = (h1r = data[1]) + data[2];
525 data[2] = h1r - data[2];
529 data[1] = c1 * ((h1r = data[1]) + data[2]);
530 data[2] = c1 * (h1r - data[2]);
531 four1(data, n >> 1, -1);
544 unsigned long i, no2;
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++)
552 twofft(data, respns, fft, ans, n);
554 for(i = 2; i <= n + 2; i += 2)
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;
564 if((mag2 = ans[i - 1] * ans[i - 1] + ans[i] * ans[i]) == 0.0)
566 qDebug(
"Attempt of deconvolving at zero response in convlv().");
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;
575 qDebug(
"No meaning for isign in convlv().");
590 int imj, ipj, j, k, kk, mm, *indx;
593 if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m)
595 qDebug(
"Inconsistent arguments detected in routine sgcoeff.");
601 for(ipj = 0; ipj <= (m << 1); ipj++)
603 sum = (ipj ? 0.0 : 1.0);
604 for(k = 1; k <= nr; k++)
606 for(k = 1; k <= nl; k++)
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;
613 for(j = 1; j <= m + 1; j++)
617 for(kk = 1; kk <= np; kk++)
619 for(k = -nl; k <= nr; k++)
623 for(mm = 1; mm <= m; mm++)
624 sum +=
b[mm + 1] * (fac *= k);
625 kk = ((np - k) % np) + 1;
642 double *y_filtered_data_p,
643 int data_point_count)
const
649#if CONVOLVE_WITH_NR_CONVLV
651 retval =
sgcoeff(
c, np, m_nL, m_nR, m_lD, m_m);
653 convlv(y_data_p, data_point_count,
c, np, 1, y_filtered_data_p);
662 qDebug() << __FILE__ << __LINE__ << __FUNCTION__ <<
"()"
672 y_filtered_data_p[k] +=
683 y_filtered_data_p[k] +=
688 for(k = data_point_count -
m_params.
nR + 1; k <= data_point_count; k++)
693 if(k + j <= data_point_count)
695 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
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.
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...
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