File indexing completed on 2024-05-12 15:27:07
0001 /*************************************************************************** 0002 File : nsl_filter.c 0003 Project : LabPlot 0004 Description : NSL Fourier filter functions 0005 -------------------------------------------------------------------- 0006 Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) 0007 0008 ***************************************************************************/ 0009 0010 /*************************************************************************** 0011 * * 0012 * This program is free software; you can redistribute it and/or modify * 0013 * it under the terms of the GNU General Public License as published by * 0014 * the Free Software Foundation; either version 2 of the License, or * 0015 * (at your option) any later version. * 0016 * * 0017 * This program is distributed in the hope that it will be useful, * 0018 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 0019 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 0020 * GNU General Public License for more details. * 0021 * * 0022 * You should have received a copy of the GNU General Public License * 0023 * along with this program; if not, write to the Free Software * 0024 * Foundation, Inc., 51 Franklin Street, Fifth Floor, * 0025 * Boston, MA 02110-1301 USA * 0026 * * 0027 ***************************************************************************/ 0028 0029 #include "nsl_filter.h" 0030 #include "nsl_common.h" 0031 #include "nsl_sf_poly.h" 0032 #include <gsl/gsl_sf_pow_int.h> 0033 #include <gsl/gsl_fft_real.h> 0034 #include <gsl/gsl_fft_halfcomplex.h> 0035 #ifdef HAVE_FFTW3 0036 #include <fftw3.h> 0037 #endif 0038 0039 const char* nsl_filter_type_name[] = { i18n("Low pass"), i18n("High pass"), i18n("Band pass"), i18n("Band reject") }; 0040 const char* nsl_filter_form_name[] = { i18n("Ideal"), i18n("Butterworth"), i18n("Chebyshev type I"), i18n("Chebyshev type II"), i18n("Legendre (Optimum L)"), i18n("Bessel (Thomson)") }; 0041 const char* nsl_filter_cutoff_unit_name[] = { i18n("Frequency"), i18n("Fraction"), i18n("Index") }; 0042 0043 /* n - order, x = w/w0 */ 0044 double nsl_filter_gain_bessel(int n, double x) { 0045 gsl_complex z0 = gsl_complex_rect(0.0, 0.0); 0046 gsl_complex z = gsl_complex_rect(0.0, x); 0047 double norm = gsl_complex_abs(nsl_sf_poly_reversed_bessel_theta(n, z)); 0048 double value = GSL_REAL(nsl_sf_poly_reversed_bessel_theta(n, z0)); 0049 return value/norm; 0050 } 0051 0052 /* size of data should be n+2 */ 0053 int nsl_filter_apply(double data[], size_t n, nsl_filter_type type, nsl_filter_form form, int order, double cutindex, double bandwidth) { 0054 if (cutindex < 0) { 0055 printf("index for cutoff must be >= 0\n"); 0056 return -1; 0057 } 0058 if (cutindex > n/2+1) { 0059 printf("index for cutoff must be <= n/2+1\n"); 0060 return -1; 0061 } 0062 0063 size_t i; 0064 double factor, centerindex = cutindex + bandwidth/2.; 0065 switch (type) { 0066 case nsl_filter_type_low_pass: 0067 switch (form) { 0068 case nsl_filter_form_ideal: 0069 for (i = (size_t)cutindex; i < n/2+1; i++) 0070 data[2*i] = data[2*i+1] = 0; 0071 break; 0072 case nsl_filter_form_butterworth: 0073 for (i = 0; i < n/2+1; i++) { 0074 factor = 1./sqrt(1. + gsl_sf_pow_int(i/cutindex, 2*order)); 0075 data[2*i] *= factor; 0076 data[2*i+1] *= factor; 0077 } 0078 break; 0079 case nsl_filter_form_chebyshev_i: 0080 for (i = 0; i < n/2+1; i++) { 0081 factor = 1./sqrt(1. + gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i/cutindex), 2)); 0082 data[2*i] *= factor; 0083 data[2*i+1] *= factor; 0084 } 0085 break; 0086 case nsl_filter_form_chebyshev_ii: 0087 for (i = 1; i < n/2+1; i++) { /* i==0: factor=1 */ 0088 factor = 1./sqrt(1. + 1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, cutindex/i), 2)); 0089 data[2*i] *= factor; 0090 data[2*i+1] *= factor; 0091 } 0092 break; 0093 case nsl_filter_form_legendre: 0094 for (i = 0; i < n/2+1; i++) { 0095 factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, i*i/(cutindex*cutindex) )); 0096 data[2*i] *= factor; 0097 data[2*i+1] *= factor; 0098 } 0099 break; 0100 case nsl_filter_form_bessel: 0101 for (i = 0; i < n/2+1; i++) { 0102 factor = nsl_filter_gain_bessel(order, i/cutindex); 0103 data[2*i] *= factor; 0104 data[2*i+1] *= factor; 0105 } 0106 break; 0107 } 0108 break; 0109 case nsl_filter_type_high_pass: 0110 switch (form) { 0111 case nsl_filter_form_ideal: 0112 for (i = 0; i < cutindex; i++) 0113 data[2*i] = data[2*i+1] = 0; 0114 break; 0115 case nsl_filter_form_butterworth: 0116 data[0]=data[1]=0; 0117 for (i = 1; i < n/2+1; i++) { 0118 factor = 1./sqrt(1.+gsl_sf_pow_int(cutindex/i, 2*order)); 0119 data[2*i] *= factor; 0120 data[2*i+1] *= factor; 0121 } 0122 break; 0123 case nsl_filter_form_chebyshev_i: 0124 data[0]=data[1]=0; 0125 for (i = 1; i < n/2+1; i++) { 0126 factor = 1./sqrt(1.+gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, cutindex/i), 2)); 0127 data[2*i] *= factor; 0128 data[2*i+1] *= factor; 0129 } 0130 break; 0131 case nsl_filter_form_chebyshev_ii: 0132 for (i = 0; i < n/2+1; i++) { 0133 factor = 1./sqrt(1.+1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i/cutindex), 2)); 0134 data[2*i] *= factor; 0135 data[2*i+1] *= factor; 0136 } 0137 break; 0138 case nsl_filter_form_legendre: 0139 data[0]=data[1]=0; 0140 for (i = 1; i < n/2+1; i++) { 0141 factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, cutindex*cutindex/(i*i) )); 0142 data[2*i] *= factor; 0143 data[2*i+1] *= factor; 0144 } 0145 break; 0146 case nsl_filter_form_bessel: 0147 data[0]=data[1]=0; 0148 for (i = 1; i < n/2+1; i++) { 0149 factor = nsl_filter_gain_bessel(order, cutindex/i); 0150 data[2*i] *= factor; 0151 data[2*i+1] *= factor; 0152 } 0153 break; 0154 } 0155 break; 0156 case nsl_filter_type_band_pass: 0157 switch (form) { 0158 case nsl_filter_form_ideal: 0159 for (i = 0; i < (size_t)cutindex; i++) 0160 data[2*i] = data[2*i+1] = 0; 0161 for (i = (size_t)(cutindex + bandwidth); i < n/2+1; i++) 0162 data[2*i] = data[2*i+1] = 0; 0163 break; 0164 case nsl_filter_form_butterworth: 0165 data[0]=data[1]=0; 0166 for (i = 1; i < n/2+1; i++) { 0167 factor = 1./sqrt(1.+gsl_sf_pow_int((i*i - centerindex*centerindex)/i/bandwidth, 2*order)); 0168 data[2*i] *= factor; 0169 data[2*i+1] *= factor; 0170 } 0171 break; 0172 case nsl_filter_form_chebyshev_i: 0173 data[0]=data[1]=0; 0174 for (i = 1; i < n/2+1; i++) { 0175 factor = 1./sqrt(1.+gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, (i*i - centerindex*centerindex)/i/bandwidth), 2)); 0176 data[2*i] *= factor; 0177 data[2*i+1] *= factor; 0178 } 0179 break; 0180 case nsl_filter_form_chebyshev_ii: 0181 for (i = 0; i < n/2+1; i++) { 0182 factor = 1./sqrt(1.+1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i*bandwidth/(i*i - centerindex*centerindex)), 2)); 0183 data[2*i] *= factor; 0184 data[2*i+1] *= factor; 0185 } 0186 break; 0187 case nsl_filter_form_legendre: 0188 data[0]=data[1]=0; 0189 for (i = 1; i < n/2+1; i++) { 0190 factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, 0191 (i*i-2.*centerindex*centerindex+gsl_pow_4(centerindex)/(i*i))/gsl_pow_2(bandwidth) )); 0192 data[2*i] *= factor; 0193 data[2*i+1] *= factor; 0194 } 0195 break; 0196 case nsl_filter_form_bessel: 0197 data[0]=data[1]=0; 0198 for (i = 1; i < n/2+1; i++) { 0199 factor = nsl_filter_gain_bessel(order, (i*i - centerindex*centerindex)/i/bandwidth); 0200 data[2*i] *= factor; 0201 data[2*i+1] *= factor; 0202 } 0203 break; 0204 } 0205 break; 0206 case nsl_filter_type_band_reject: 0207 switch (form) { 0208 case nsl_filter_form_ideal: 0209 for (i = (size_t)cutindex; i < (size_t)(cutindex + bandwidth); i++) 0210 data[2*i] = data[2*i+1] = 0; 0211 break; 0212 case nsl_filter_form_butterworth: 0213 for (i = 0; i < n/2+1; i++) { 0214 factor = 1./sqrt(1.+gsl_sf_pow_int(i*bandwidth/(i*i - centerindex*centerindex), 2*order)); 0215 data[2*i] *= factor; 0216 data[2*i+1] *= factor; 0217 } 0218 break; 0219 case nsl_filter_form_chebyshev_i: 0220 for (i = 0; i < n/2+1; i++) { 0221 factor = 1./sqrt(1.+gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, i*bandwidth/(i*i - centerindex*centerindex)), 2)); 0222 data[2*i] *= factor; 0223 data[2*i+1] *= factor; 0224 } 0225 break; 0226 case nsl_filter_form_chebyshev_ii: 0227 for (i = 1; i < n/2+1; i++) { /* i==0: factor=1 */ 0228 factor = 1./sqrt(1.+1./gsl_sf_pow_int(nsl_sf_poly_chebyshev_T(order, (i*i - centerindex*centerindex)/i/bandwidth), 2)); 0229 data[2*i] *= factor; 0230 data[2*i+1] *= factor; 0231 } 0232 break; 0233 case nsl_filter_form_legendre: 0234 for (i = 0; i < n/2+1; i++) { 0235 factor = 1./sqrt(1. + nsl_sf_poly_optimal_legendre_L(order, 0236 gsl_pow_2(i*bandwidth)/gsl_pow_2(i*i-centerindex*centerindex) )); 0237 data[2*i] *= factor; 0238 data[2*i+1] *= factor; 0239 } 0240 break; 0241 case nsl_filter_form_bessel: 0242 for (i = 0; i < n/2+1; i++) { 0243 if (i == centerindex) 0244 factor = 0; 0245 else 0246 factor = nsl_filter_gain_bessel(order, i*bandwidth/(i*i - centerindex*centerindex)); 0247 data[2*i] *= factor; 0248 data[2*i+1] *= factor; 0249 } 0250 break; 0251 } 0252 break; 0253 } 0254 0255 return 0; 0256 } 0257 0258 void print_fdata(double data[], size_t n) { 0259 size_t i; 0260 for(i=0; i < 2*(n/2+1); i++) 0261 printf("%g ", data[i]); 0262 printf("\nreal: "); 0263 for(i=0; i < n/2+1; i++) 0264 printf("%g ", data[2*i]); 0265 printf("\nimag: "); 0266 for(i=0; i < n/2+1; i++) 0267 printf("%g ", data[2*i+1]); 0268 puts(""); 0269 } 0270 0271 int nsl_filter_fourier(double data[], size_t n, nsl_filter_type type, nsl_filter_form form, int order, int cutindex, int bandwidth) { 0272 /* 1. transform */ 0273 double* fdata = (double*)malloc(2*n*sizeof(double)); /* contains re0,im0,re1,im1,re2,im2,... */ 0274 #ifdef HAVE_FFTW3 0275 fftw_plan plan = fftw_plan_dft_r2c_1d((int)n, data, (fftw_complex *) fdata, FFTW_ESTIMATE); 0276 fftw_execute(plan); 0277 fftw_destroy_plan(plan); 0278 #else 0279 gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(n); 0280 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(n); 0281 0282 gsl_fft_real_transform(data, 1, n, real, work); 0283 gsl_fft_real_wavetable_free(real); 0284 0285 gsl_fft_halfcomplex_unpack(data, fdata, 1, n); 0286 #endif 0287 0288 /* 2. apply filter */ 0289 /*print_fdata(fdata, n);*/ 0290 int status = nsl_filter_apply(fdata, n, type, form, order, cutindex, bandwidth); 0291 /*print_fdata(fdata, n);*/ 0292 0293 /* 3. back transform */ 0294 #ifdef HAVE_FFTW3 0295 plan = fftw_plan_dft_c2r_1d((int)n, (fftw_complex *) fdata, data, FFTW_ESTIMATE); 0296 fftw_execute(plan); 0297 fftw_destroy_plan(plan); 0298 /* normalize*/ 0299 size_t i; 0300 for (i=0; i < n; i++) 0301 data[i] /= n; 0302 #else 0303 gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc(n); 0304 gsl_fft_halfcomplex_inverse(data, 1, n, hc, work); 0305 gsl_fft_halfcomplex_wavetable_free(hc); 0306 gsl_fft_real_workspace_free (work); 0307 #endif 0308 free(fdata); 0309 0310 return status; 0311 }