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; _______________________________________________ kde-doc-english mailing list [email protected] https://mail.kde.org/mailman/listinfo/kde-doc-english
