55 int nL,
int nR,
int m,
int lD,
bool convolveWithNr)
106 int data_point_count = data_points.size();
113 y_filtered_data_p =
dvector(1, 2 * data_point_count);
115 y_filtered_data_p =
dvector(1, data_point_count);
117 for(
int iter = 0; iter < data_point_count; ++iter)
119 x_data_p[iter] = data_points.at(iter).x;
120 y_initial_data_p[iter] = data_points.at(iter).y;
125 runFilter(y_initial_data_p, y_filtered_data_p, data_point_count);
128 auto iter_yf = y_filtered_data_p;
129 for(
auto &data_point : data_points)
131 data_point.y = *iter_yf;
150 return QString::asprintf(
151 "Savitzy-Golay filter parameters:\n"
152 "nL: %d ; nR: %d ; m: %d ; lD: %d ; convolveWithNr : %s",
165 v = (
int *)malloc((
size_t)((nh - nl + 2) *
sizeof(
int)));
168 qFatal(
"Error: Allocation failure.");
182 qFatal(
"Error: Allocation failure.");
184 for(k = nl; k <= nh; k++)
193 long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
198 qFatal(
"Error: Allocation failure.");
206 qFatal(
"Error: Allocation failure.");
210 for(i = nrl + 1; i <= nrh; i++)
211 m[i] = m[i - 1] + ncol;
219 [[maybe_unused]]
long nh)
const
221 free((
char *)(v + nl - 1));
228 [[maybe_unused]]
long nh)
const
230 free((
char *)(v + nl - 1));
237 [[maybe_unused]]
long nrh,
239 [[maybe_unused]]
long nch)
const
241 free((
char *)(m[nrl] + ncl - 1));
242 free((
char *)(m + nrl - 1));
252 int i, ii = 0, ip, j;
255 for(i = 1; i <= n; i++)
261 for(j = ii; j <= i - 1; j++)
262 sum -=
a[i][j] *
b[j];
267 for(i = n; i >= 1; i--)
270 for(j = i + 1; j <= n; j++)
271 sum -=
a[i][j] *
b[j];
272 b[i] =
sum /
a[i][i];
283 int i, imax = 0, j, k;
289 for(i = 1; i <= n; i++)
292 for(j = 1; j <= n; j++)
293 if((temp = fabs(
a[i][j])) > big)
297 qFatal(
"Error: Singular matrix found in routine ludcmp().");
301 for(j = 1; j <= n; j++)
303 for(i = 1; i < j; i++)
306 for(k = 1; k < i; k++)
307 sum -=
a[i][k] *
a[k][j];
311 for(i = j; i <= n; i++)
314 for(k = 1; k < j; k++)
315 sum -=
a[i][k] *
a[k][j];
317 if((dum = vv[i] * fabs(
sum)) >= big)
325 for(k = 1; k <= n; k++)
328 a[imax][k] =
a[j][k];
336 a[j][j] = std::numeric_limits<pappso_double>::epsilon();
339 dum = 1.0 / (
a[j][j]);
340 for(i = j + 1; i <= n; i++)
351 unsigned long n, mmax, m, j, istep, i;
357 for(i = 1; i < n; i += 2)
361 SWAP(data[j], data[i]);
362 SWAP(data[j + 1], data[i + 1]);
365 while(m >= 2 && j > m)
376 theta = isign * (6.28318530717959 / mmax);
377 wtemp = sin(0.5 * theta);
378 wpr = -2.0 * wtemp * wtemp;
382 for(m = 1; m < mmax; m += 2)
384 for(i = m; i <= n; i += istep)
387 tempr = wr * data[j] - wi * data[j + 1];
388 tempi = wr * data[j + 1] + wi * data[j];
389 data[j] = data[i] - tempr;
390 data[j + 1] = data[i + 1] - tempi;
392 data[i + 1] += tempi;
394 wr = (wtemp = wr) * wpr - wi * wpi + wr;
395 wi = wi * wpr + wtemp * wpi + wi;
409 unsigned long nn3, nn2, jj, j;
412 nn3 = 1 + (nn2 = 2 + n + n);
413 for(j = 1, jj = 2; j <= n; j++, jj += 2)
415 fft1[jj - 1] = data1[j];
420 fft1[2] = fft2[2] = 0.0;
421 for(j = 3; j <= n + 1; j += 2)
423 rep = 0.5 * (fft1[j] + fft1[nn2 - j]);
424 rem = 0.5 * (fft1[j] - fft1[nn2 - j]);
425 aip = 0.5 * (fft1[j + 1] + fft1[nn3 - j]);
426 aim = 0.5 * (fft1[j + 1] - fft1[nn3 - j]);
430 fft1[nn3 - j] = -aim;
442 unsigned long i, i1, i2, i3, i4, np3;
450 four1(data, n >> 1, 1);
457 wtemp = sin(0.5 * theta);
458 wpr = -2.0 * wtemp * wtemp;
463 for(i = 2; i <= (n >> 2); i++)
465 i4 = 1 + (i3 = np3 - (i2 = 1 + (i1 = i + i - 1)));
466 h1r = c1 * (data[i1] + data[i3]);
467 h1i = c1 * (data[i2] - data[i4]);
468 h2r = -c2 * (data[i2] + data[i4]);
469 h2i = c2 * (data[i1] - data[i3]);
470 data[i1] = h1r + wr * h2r - wi * h2i;
471 data[i2] = h1i + wr * h2i + wi * h2r;
472 data[i3] = h1r - wr * h2r + wi * h2i;
473 data[i4] = -h1i + wr * h2i + wi * h2r;
474 wr = (wtemp = wr) * wpr - wi * wpi + wr;
475 wi = wi * wpr + wtemp * wpi + wi;
479 data[1] = (h1r = data[1]) + data[2];
480 data[2] = h1r - data[2];
484 data[1] = c1 * ((h1r = data[1]) + data[2]);
485 data[2] = c1 * (h1r - data[2]);
486 four1(data, n >> 1, -1);
499 unsigned long i, no2;
503 for(i = 1; i <= (m - 1) / 2; i++)
504 respns[n + 1 - i] = respns[m + 1 - i];
505 for(i = (m + 3) / 2; i <= n - (m - 1) / 2; i++)
507 twofft(data, respns, fft, ans, n);
509 for(i = 2; i <= n + 2; i += 2)
514 (fft[i - 1] * (dum = ans[i - 1]) - fft[i] * ans[i]) / no2;
515 ans[i] = (fft[i] * dum + fft[i - 1] * ans[i]) / no2;
519 if((mag2 = ans[i - 1] * ans[i - 1] + ans[i] * ans[i]) == 0.0)
521 qDebug(
"Attempt of deconvolving at zero response in convlv().");
525 (fft[i - 1] * (dum = ans[i - 1]) + fft[i] * ans[i]) / mag2 / no2;
526 ans[i] = (fft[i] * dum - fft[i - 1] * ans[i]) / mag2 / no2;
530 qDebug(
"No meaning for isign in convlv().");
543 pappso_double c[],
int np,
int nl,
int nr,
int ld,
int m)
const
545 int imj, ipj, j, k, kk, mm, *indx;
548 if(np < nl + nr + 1 || nl < 0 || nr < 0 || ld > m || nl + nr < m)
550 qDebug(
"Inconsistent arguments detected in routine sgcoeff.");
556 for(ipj = 0; ipj <= (m << 1); ipj++)
558 sum = (ipj ? 0.0 : 1.0);
559 for(k = 1; k <= nr; k++)
561 for(k = 1; k <= nl; k++)
563 mm = (ipj < 2 * m - ipj ? ipj : 2 * m - ipj);
564 for(imj = -mm; imj <= mm; imj += 2)
565 a[1 + (ipj + imj) / 2][1 + (ipj - imj) / 2] =
sum;
568 for(j = 1; j <= m + 1; j++)
572 for(kk = 1; kk <= np; kk++)
574 for(k = -nl; k <= nr; k++)
578 for(mm = 1; mm <= m; mm++)
579 sum +=
b[mm + 1] * (fac *= k);
580 kk = ((np - k) % np) + 1;
596 double *y_filtered_data_p,
597 int data_point_count)
const
603 #if CONVOLVE_WITH_NR_CONVLV
607 convlv(y_data_p, data_point_count,
c, np, 1, y_filtered_data_p);
616 qDebug() << __FILE__ << __LINE__ << __FUNCTION__ <<
"()"
619 for(k = 1; k <=
m_nL; k++)
621 for(y_filtered_data_p[k] = 0.0, j = -
m_nL; j <=
m_nR; j++)
625 y_filtered_data_p[k] +=
626 c[(j >= 0 ? j + 1 :
m_nR +
m_nL + 2 + j)] * y_data_p[k + j];
630 for(k =
m_nL + 1; k <= data_point_count -
m_nR; k++)
632 for(y_filtered_data_p[k] = 0.0, j = -
m_nL; j <=
m_nR; j++)
634 y_filtered_data_p[k] +=
635 c[(j >= 0 ? j + 1 :
m_nR +
m_nL + 2 + j)] * y_data_p[k + j];
638 for(k = data_point_count -
m_nR + 1; k <= data_point_count; k++)
640 for(y_filtered_data_p[k] = 0.0, j = -
m_nL; j <=
m_nR; j++)
642 if(k + j <= data_point_count)
644 y_filtered_data_p[k] +=
645 c[(j >= 0 ? j + 1 :
m_nR +
m_nL + 2 + j)] * y_data_p[k + j];