[kde-doc-english] [labplot] /: added Bessel form for Fourier filter
Stefan Gerlach
stefan.gerlach at uni-konstanz.de
Fri Jul 15 21:38:14 UTC 2016
Git commit 711dd489d0e2406433ccc093f7c447decf240789 by Stefan Gerlach.
Committed on 15/07/2016 at 21:38.
Pushed by sgerlach into branch 'master'.
added Bessel form for Fourier filter
M +1 -1 ChangeLog
M +4 -3 doc/index.docbook
M +43 -5 src/backend/nsl/nsl_filter.c
M +5 -2 src/backend/nsl/nsl_filter.h
M +8 -1 src/backend/nsl/nsl_filter_test.c
M +43 -19 src/backend/nsl/nsl_sf_poly.c
M +7 -1 src/backend/nsl/nsl_sf_poly.h
M +1 -0 src/kdefrontend/dockwidgets/XYFourierFilterCurveDock.cpp
http://commits.kde.org/labplot/711dd489d0e2406433ccc093f7c447decf240789
diff --git a/ChangeLog b/ChangeLog
index a59a3a7..0c42a5f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -5,7 +5,7 @@ New features:
* Export of spreadsheets and matrices to LaTeX tables
* Interpolation of data including different splines, cosine, exponential, cubic Hermite (Catmull-Rom, cardinal, Kochanek-Bartels) and rational functions
* Data smoothing using moving average (centered or lagged), percentile filter or Savitzky-Golay algorithm
- * Fourier filter (low pass, high pass, band pass, band reject) with ideal, Butterworth, Chebychev I+II or Legendre filter
+ * Fourier filter (low pass, high pass, band pass, band reject) with ideal, Butterworth, Chebychev I+II, Legendre or Bessel-Thomson filter
* Fourier transform with many window functions (Welch, Hann, Hamming, etc.) calculating magnitude, amplitude, power, phase, dB, etc. and supporting
one/two sided spectrum with or without shift and x scaling to frequency, index or period
* Filter and search capabilities in the drop down box for the selection of data sources
diff --git a/doc/index.docbook b/doc/index.docbook
index bef5754..5543b66 100644
--- a/doc/index.docbook
+++ b/doc/index.docbook
@@ -51,8 +51,8 @@
</copyright>
<legalnotice>&FDLNotice;</legalnotice>
-<date>2016-07-02</date>
-<releaseinfo>3.2.2</releaseinfo>
+<date>2016-07-15</date>
+<releaseinfo>3.2.3</releaseinfo>
<abstract>
<para>
@@ -869,7 +869,8 @@ The menu is only available when a datapicker object is selected on the <guilabel
<listitem><para>Ideal</para></listitem>
<listitem><para>Butterworth (order 1 to 10)</para></listitem>
<listitem><para>Chebyshev type I or II (order 1 to 10)</para></listitem>
- <listitem><para>Optimal Legendre (order 1 to 10)</para></listitem>
+ <listitem><para>Optimal "L"egendre (order 1 to 10)</para></listitem>
+ <listitem><para>Bessel-Thomson (any order)</para></listitem>
</itemizedlist>
<para>
The cutoff value(s) can be specified in the units frequency (Hertz), fraction (0.0 to 1.0) or index
diff --git a/src/backend/nsl/nsl_filter.c b/src/backend/nsl/nsl_filter.c
index 58d6c2d..ae5c8a9 100644
--- a/src/backend/nsl/nsl_filter.c
+++ b/src/backend/nsl/nsl_filter.c
@@ -38,9 +38,14 @@
#endif
const char* nsl_filter_type_name[] = { "Low pass", "High pass", "Band pass", "Band reject" };
-const char* nsl_filter_form_name[] = { "Ideal", "Butterworth", "Chebyshev type I", "Chebyshev type II", "Legendre (Optimum L)" };
+const char* nsl_filter_form_name[] = { "Ideal", "Butterworth", "Chebyshev type I", "Chebyshev type II", "Legendre (Optimum L)", "Bessel (Thomson)" };
const char* nsl_filter_cutoff_unit_name[] = { "Frequency", "Fraction", "Index" };
+/* n - order, x = w/w0 */
+double nsl_filter_gain_bessel(int n, double x) {
+ return nsl_sf_poly_reversed_bessel_theta(n, 0)/cabs(nsl_sf_poly_reversed_bessel_theta(n, I*x));
+}
+
int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_form form, int order, double cutindex, double bandwidth) {
if (cutindex < 0) {
printf("index for cutoff must be >= 0\n");
@@ -83,7 +88,14 @@ int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_f
break;
case nsl_filter_form_legendre:
for (i = 0; i < n/2+1; i++) {
- factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre(order, i*i/(cutindex*cutindex) ));
+ factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, i*i/(cutindex*cutindex) ));
+ data[2*i] *= factor;
+ data[2*i+1] *= factor;
+ }
+ break;
+ case nsl_filter_form_bessel:
+ for (i = 0; i < n/2+1; i++) {
+ factor = nsl_filter_gain_bessel(order, i/cutindex);
data[2*i] *= factor;
data[2*i+1] *= factor;
}
@@ -122,7 +134,15 @@ int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_f
case nsl_filter_form_legendre:
data[0]=data[1]=0;
for (i = 1; i < n/2+1; i++) {
- factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre(order, cutindex*cutindex/(i*i) ));
+ factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, cutindex*cutindex/(i*i) ));
+ data[2*i] *= factor;
+ data[2*i+1] *= factor;
+ }
+ break;
+ case nsl_filter_form_bessel:
+ data[0]=data[1]=0;
+ for (i = 1; i < n/2+1; i++) {
+ factor = nsl_filter_gain_bessel(order, cutindex/i);
data[2*i] *= factor;
data[2*i+1] *= factor;
}
@@ -163,12 +183,20 @@ int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_f
case nsl_filter_form_legendre:
data[0]=data[1]=0;
for (i = 1; i < n/2+1; i++) {
- factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre(order,
+ factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order,
(i*i-2.*centerindex*centerindex+gsl_pow_4(centerindex)/(i*i))/gsl_pow_2(bandwidth) ));
data[2*i] *= factor;
data[2*i+1] *= factor;
}
break;
+ case nsl_filter_form_bessel:
+ data[0]=data[1]=0;
+ for (i = 1; i < n/2+1; i++) {
+ factor = nsl_filter_gain_bessel(order, (i*i - centerindex*centerindex)/i/bandwidth);
+ data[2*i] *= factor;
+ data[2*i+1] *= factor;
+ }
+ break;
}
break;
case nsl_filter_type_band_reject:
@@ -200,12 +228,22 @@ int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_f
break;
case nsl_filter_form_legendre:
for (i = 0; i < n/2+1; i++) {
- factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre(order,
+ factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order,
gsl_pow_2(i*bandwidth)/gsl_pow_2(i*i-centerindex*centerindex) ));
data[2*i] *= factor;
data[2*i+1] *= factor;
}
break;
+ case nsl_filter_form_bessel:
+ for (i = 0; i < n/2+1; i++) {
+ if (i == centerindex)
+ factor = 0;
+ else
+ factor = nsl_filter_gain_bessel(order, i*bandwidth/(i*i - centerindex*centerindex));
+ data[2*i] *= factor;
+ data[2*i+1] *= factor;
+ }
+ break;
}
break;
}
diff --git a/src/backend/nsl/nsl_filter.h b/src/backend/nsl/nsl_filter.h
index 2732093..effdbbd 100644
--- a/src/backend/nsl/nsl_filter.h
+++ b/src/backend/nsl/nsl_filter.h
@@ -35,9 +35,9 @@
typedef enum {nsl_filter_type_low_pass, nsl_filter_type_high_pass, nsl_filter_type_band_pass,
nsl_filter_type_band_reject} nsl_filter_type; /*TODO: Threshold */
extern const char* nsl_filter_type_name[];
-#define NSL_FILTER_FORM_COUNT 5
+#define NSL_FILTER_FORM_COUNT 6
typedef enum {nsl_filter_form_ideal, nsl_filter_form_butterworth, nsl_filter_form_chebyshev_i,
- nsl_filter_form_chebyshev_ii, nsl_filter_form_legendre} nsl_filter_form; /*TODO: Gaussian, Bessel, ... */
+ nsl_filter_form_chebyshev_ii, nsl_filter_form_legendre, nsl_filter_form_bessel} nsl_filter_form;
extern const char* nsl_filter_form_name[];
/* unit for cutoff
Frequency=0..f_max, Fraction=0..1, Index=0..N-1
@@ -47,6 +47,9 @@ typedef enum {nsl_filter_cutoff_unit_frequency, nsl_filter_cutoff_unit_fraction,
nsl_filter_cutoff_unit_index} nsl_filter_cutoff_unit;
extern const char* nsl_filter_cutoff_unit_name[];
+/* Gain G(x) for Bessel filter */
+double nsl_filter_gain_bessel(int n, double x);
+
int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_form form, int order, double cutindex, double bandwidth);
int nsl_filter_fourier(double data[], size_t n, nsl_filter_type type, nsl_filter_form form, int order, int cutindex, int bandwidth);
diff --git a/src/backend/nsl/nsl_filter_test.c b/src/backend/nsl/nsl_filter_test.c
index f891488..3379c27 100644
--- a/src/backend/nsl/nsl_filter_test.c
+++ b/src/backend/nsl/nsl_filter_test.c
@@ -48,11 +48,18 @@ int main() {
for(i=0;i<N;i++)
data[i]=1.0;
+ /*Bessel gain*/
+/* printf("G = %g\n",nsl_filter_gain_bessel(3,1)); */
+
/* filter form */
/*nsl_filter_apply(data, N, nsl_filter_type_low_pass, nsl_filter_form_legendre, 2, 50, 2);*/
/*nsl_filter_apply(data, N, nsl_filter_type_high_pass, nsl_filter_form_legendre, 2, 50, 2);*/
/*nsl_filter_apply(data, N, nsl_filter_type_band_pass, nsl_filter_form_legendre, 2, 100, 50);*/
- nsl_filter_apply(data, N, nsl_filter_type_band_reject, nsl_filter_form_legendre, 2, 100, 100);
+ /*nsl_filter_apply(data, N, nsl_filter_type_band_reject, nsl_filter_form_legendre, 2, 100, 100);*/
+ /*nsl_filter_apply(data, N, nsl_filter_type_low_pass, nsl_filter_form_bessel, 2, 100, 100);*/
+ /*nsl_filter_apply(data, N, nsl_filter_type_high_pass, nsl_filter_form_bessel, 2, 100, 100);*/
+ /*nsl_filter_apply(data, N, nsl_filter_type_band_pass, nsl_filter_form_bessel, 2, 100, 100);*/
+ nsl_filter_apply(data, N, nsl_filter_type_band_reject, nsl_filter_form_bessel, 2, 100, 100);
for(i=0; i < N/2; i++)
printf("%d %g\n", i, data[2*i]);
diff --git a/src/backend/nsl/nsl_sf_poly.c b/src/backend/nsl/nsl_sf_poly.c
index 414b962..31bccf5 100644
--- a/src/backend/nsl/nsl_sf_poly.c
+++ b/src/backend/nsl/nsl_sf_poly.c
@@ -33,48 +33,72 @@
/* see https://en.wikipedia.org/wiki/Chebyshev_polynomials */
double nsl_sf_poly_chebyshev_T(int n, double x) {
- if(fabs(x) <= 1)
- return cos(n*acos(x));
+ if (fabs(x) <= 1)
+ return cos(n * acos(x));
else if (x > 1)
- return cosh(n*gsl_acosh(x));
+ return cosh(n * gsl_acosh(x));
else
- return pow(-1.,n)*cosh(n*gsl_acosh(-x));
+ return pow(-1., n)*cosh(n * gsl_acosh(-x));
}
double nsl_sf_poly_chebyshev_U(int n, double x) {
- double sq=sqrt(x*x-1);
- return (pow(x+sq,n+1)-pow(x-sq,n+1))/2./sq;
+ double sq = sqrt(x*x - 1);
+ return (pow(x + sq, n + 1) - pow(x - sq, n + 1))/2./sq;
}
/* from http://www.crbond.com/papers/lopt.pdf */
-double nsl_sf_poly_optimal_legendre(int n, double x) {
- if (n<1 || n>10)
+double nsl_sf_poly_optimal_legendre_L(int n, double x) {
+ if (n < 1 || n > 10)
return 0.0;
switch (n) {
case 1:
return x;
case 2:
- return gsl_pow_2(x);
+ return x*x;
case 3:
- return x - 3.*gsl_pow_2(x) + 3.*gsl_pow_3(x);
+ return (1. + (-3. + 3.*x)*x)*x;
case 4:
- return 3.*gsl_pow_2(x) - 8.*gsl_pow_3(x) + 6.*gsl_pow_4(x);
+ return (3. + (-8. + 6*x)*x)*x*x;
case 5:
- return x - 8.*gsl_pow_2(x) + 28.*gsl_pow_3(x) - 40.*gsl_pow_4(x) + 20.*gsl_pow_5(x);
+ return (1. + (-8. + (28. + (-40. + 20*x)*x)*x)*x)*x;
case 6:
- return 6.*gsl_pow_2(x) - 40.*gsl_pow_3(x) + 105.*gsl_pow_4(x) - 120.*gsl_pow_5(x) + 50.*gsl_pow_6(x);
+ return (6. + (-40. + (105. + (-120. + 50.*x)*x)*x)*x)*x*x;
case 7:
- return x - 15.*gsl_pow_2(x) + 105.*gsl_pow_3(x) - 355.*gsl_pow_4(x) + 615.*gsl_pow_5(x) - 525.*gsl_pow_6(x) + 175.*gsl_pow_7(x);
+ return (1. + (-15. + (105. + (-355. + (615. + (-525. + 175.*x)*x)*x)*x)*x)*x)*x;
case 8:
- return 10.*gsl_pow_2(x) - 120.*gsl_pow_3(x) + 615.*gsl_pow_4(x) - 1624.*gsl_pow_5(x) + 2310.*gsl_pow_6(x) - 1680.*gsl_pow_7(x) + 490.*gsl_pow_8(x);
+ return (10. + (-120. + (615. + (-1624. + (2310. + (-1680. + 490.*x)*x)*x)*x)*x)*x)*x*x;
case 9:
- return x - 24.*gsl_pow_2(x) + 276.*gsl_pow_3(x) - 1624.*gsl_pow_4(x) + 5376.*gsl_pow_5(x) - 10416.*gsl_pow_6(x) + 11704.*gsl_pow_7(x)
- - 7056.*gsl_pow_8(x) + 1764.*gsl_pow_9(x);
+ return (1. + (-24. + (276. + (-1624. + (5376. + (-10416. + (11704. + (-7056. + 1764*x)*x)*x)*x)*x)*x)*x)*x)*x;
case 10:
- return 15.*gsl_pow_2(x) - 280.*gsl_pow_3(x) + 2310.*gsl_pow_4(x) - 10416.*gsl_pow_5(x) + 27860.*gsl_pow_6(x) - 45360.*gsl_pow_7(x)
- + 44100.*gsl_pow_8(x) - 23520.*gsl_pow_9(x) + 5292.*gsl_pow_5(gsl_pow_2(x));
+ return (15. + (-280. + (2310. + (-10416. + (27860. + (-45360. + (44100. + (23520. + 5292.*x)*x)*x)*x)*x)*x)*x)*x)*x*x;
}
return 0.0;
}
+
+/*
+ * https://en.wikipedia.org/wiki/Bessel_polynomials
+ * using recursion
+*/
+double complex nsl_sf_poly_bessel_y(int n, double complex x) {
+ if (n == 0)
+ return 1.0;
+ else if (n == 1)
+ return 1.0 + x;
+
+ return (2*n - 1)*x*nsl_sf_poly_bessel_y(n - 1, x) + nsl_sf_poly_bessel_y(n - 2, x);
+}
+
+/*
+ * https://en.wikipedia.org/wiki/Bessel_polynomials
+ * using recursion
+*/
+double complex nsl_sf_poly_reversed_bessel_theta(int n, double complex x) {
+ if (n == 0)
+ return 1.0;
+ else if (n == 1)
+ return 1.0 + x;
+
+ return (2*n - 1)*nsl_sf_poly_reversed_bessel_theta(n - 1, x) + x*x*nsl_sf_poly_reversed_bessel_theta(n - 2, x);
+}
diff --git a/src/backend/nsl/nsl_sf_poly.h b/src/backend/nsl/nsl_sf_poly.h
index c5ead2f..5fa2caa 100644
--- a/src/backend/nsl/nsl_sf_poly.h
+++ b/src/backend/nsl/nsl_sf_poly.h
@@ -30,6 +30,7 @@
#define NSL_SF_POLY_H
#include <stdlib.h>
+#include <complex.h>
/* Chebychev T_n(x) */
double nsl_sf_poly_chebyshev_T(int n, double x);
@@ -37,6 +38,11 @@ double nsl_sf_poly_chebyshev_T(int n, double x);
double nsl_sf_poly_chebyshev_U(int n, double x);
/* Optimal "L"egendre Polynomials */
-double nsl_sf_poly_optimal_legendre(int n, double x);
+double nsl_sf_poly_optimal_legendre_L(int n, double x);
+
+/* Bessel polynomials y_n(x) */
+double complex nsl_sf_poly_bessel_y(int n, double complex x);
+/* reversed Bessel polynomials \theta_n(x) */
+double complex nsl_sf_poly_reversed_bessel_theta(int n, double complex x);
#endif /* NSL_SF_POLY_H */
diff --git a/src/kdefrontend/dockwidgets/XYFourierFilterCurveDock.cpp b/src/kdefrontend/dockwidgets/XYFourierFilterCurveDock.cpp
index b6f8ef9..275e0ce 100644
--- a/src/kdefrontend/dockwidgets/XYFourierFilterCurveDock.cpp
+++ b/src/kdefrontend/dockwidgets/XYFourierFilterCurveDock.cpp
@@ -300,6 +300,7 @@ void XYFourierFilterCurveDock::formChanged() {
case nsl_filter_form_chebyshev_i:
case nsl_filter_form_chebyshev_ii:
case nsl_filter_form_legendre:
+ case nsl_filter_form_bessel:
uiGeneralTab.sbOrder->setVisible(true);
uiGeneralTab.lOrder->setVisible(true);
break;
More information about the kde-doc-english
mailing list