diff --git a/dspl/src/filter_ft.c b/dspl/src/filter_ft.c index 01d888a..4058e57 100644 --- a/dspl/src/filter_ft.c +++ b/dspl/src/filter_ft.c @@ -53,6 +53,34 @@ int DSPL_API low2bp(double* b, double* a, int ord, + +/****************************************************************************** + * low 2 bandstop transformation + ******************************************************************************/ +int DSPL_API low2bs(double* b, double* a, int ord, + double w0, double wsl, double wsh, + double* beta, double* alpha) +{ + + double den[3] = {0.0, 0.0, 1.0}; + double num[3] = {0.0, 0.0, 0.0}; + + if(!b || !a || !beta || !alpha) + return ERROR_PTR; + if(ord < 1) + return ERROR_FILTER_ORD; + if(w0 <= 0.0 || wsl <= 0.0 || wsh <= 0.0 || wsh <= wsl) + return ERROR_FILTER_FT; + + den[0] = (wsh * wsl) / (w0 * w0); + num[1] = (wsh - wsl) / w0; + + return ratcompos(b, a, ord, num, den, 2, beta, alpha); +} + + + + /****************************************************************************** * low 2 high transformation ******************************************************************************/ diff --git a/dspl/src/filter_iir.c b/dspl/src/filter_iir.c index df2413b..0aaae0f 100644 --- a/dspl/src/filter_iir.c +++ b/dspl/src/filter_iir.c @@ -147,8 +147,10 @@ int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, err = low2bp(bs, as, ord_ap, 1.0, wa0, wa1, bt, at); break; - // case DSPL_FILTER_BSTOP: - // break; + case DSPL_FILTER_BSTOP: + // need frequency transform ws -> 1 rad/s + err = low2bs(bs, as, ord_ap, 1.0, wa0, wa1, bt, at); + break; default: err = ERROR_FILTER_TYPE; diff --git a/examples/bin/gnuplot/iir_bstop.plt b/examples/bin/gnuplot/iir_bstop.plt new file mode 100644 index 0000000..bc79888 --- /dev/null +++ b/examples/bin/gnuplot/iir_bstop.plt @@ -0,0 +1,21 @@ +unset key +set grid +set xlabel " normalized frequency" + +set terminal pngcairo size 920, 260 enhanced font 'Verdana,8' +set output 'img/iir_bstop.png' + +set multiplot layout 1,3 rowsfirst +set ylabel "Magnitude, dB" +set yrange [-100:5] +plot 'dat/iir_bstop_mag.txt' with lines + +set ylabel "Phase response, rad" +unset yrange +plot 'dat/iir_bstop_phi.txt' with lines + +set ylabel "Groupdelay, samples" +unset yrange +plot 'dat/iir_bstop_tau.txt' with lines + +unset multiplot \ No newline at end of file diff --git a/examples/src/iir_bstop.c b/examples/src/iir_bstop.c new file mode 100644 index 0000000..3732352 --- /dev/null +++ b/examples/src/iir_bstop.c @@ -0,0 +1,61 @@ +#include +#include +#include +#include "dspl.h" + +// Filter order (must be even for bandstop IIR) +#define ORD 14 + +// Frequency response vector size +#define N 1001 + + +int main() +{ + void* handle; // DSPL handle + handle = dspl_load(); // Load DSPL function + + double a[ORD+1], b[ORD+1]; // H(s) coefficients + double rs = 60.0; // Bandstop suppression equals 60 dB + double rp = 1.0; // Bandpass ripple equals 1 dB + + // Frequency (w), magnitude (mag), phase response (phi) + // and group delay (tau) + double w[N], mag[N], phi[N], tau[N]; + int k; + + // Calculate Chebyshev type 2 digital bandstop filter + int res = iir(rp, rs, ORD, 0.3, 0.7, + DSPL_FILTER_CHEBY2 | DSPL_FILTER_BSTOP, b, a); + if(res != RES_OK) + printf("error code = 0x%8x\n", res); + + // Print coefficients + for(k = 0; k < ORD+1; k++) + printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]); + + // Normalized circular frequency from 0 to pi + linspace(0, M_PI, N , DSPL_PERIODIC, w); + + // Filter magnitude, phase response and group delay + filter_freq_resp(b, a, ORD, w, N, DSPL_FLAG_LOGMAG|DSPL_FLAG_UNWRAP, + mag, phi, tau); + + // Normalized frequency from 0 to 1. + // w = 1 corresponds to Fs/2 + for(k = 0; k < N; k++) + w[k] /= M_PI; + + // Save filter frequency response to the txt-files + // for plotting by GNUPLOT + writetxt(w, mag, N, "dat/iir_bstop_mag.txt"); + writetxt(w, phi, N, "dat/iir_bstop_phi.txt"); + writetxt(w, tau, N, "dat/iir_bstop_tau.txt"); + + dspl_free(handle); // free dspl handle + + // run GNUPLOT script + return system("gnuplot -p gnuplot/iir_bstop.plt");; +} + + diff --git a/include/dspl.c b/include/dspl.c index 8c990c5..fa6fed5 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -107,6 +107,7 @@ p_linspace linspace ; p_log_cmplx log_cmplx ; p_logspace logspace ; p_low2bp low2bp ; +p_low2bs low2bs ; p_low2high low2high ; p_low2low low2low ; p_matrix_create matrix_create ; @@ -269,6 +270,7 @@ void* dspl_load() LOAD_FUNC(log_cmplx); LOAD_FUNC(logspace); LOAD_FUNC(low2bp); + LOAD_FUNC(low2bs); LOAD_FUNC(low2high); LOAD_FUNC(low2low); LOAD_FUNC(matrix_create); diff --git a/include/dspl.h b/include/dspl.h index 4e4f814..190eaf0 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -660,6 +660,15 @@ DECLARE_FUNC(int, low2bp, double* b COMMA double* beta COMMA double* alpha); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, low2bs, double* b + COMMA double* a + COMMA int ord + COMMA double w0 + COMMA double wsl + COMMA double wsh + COMMA double* beta + COMMA double* alpha); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, low2high, double* b COMMA double* a COMMA int ord diff --git a/release/include/dspl.c b/release/include/dspl.c index 8c990c5..fa6fed5 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -107,6 +107,7 @@ p_linspace linspace ; p_log_cmplx log_cmplx ; p_logspace logspace ; p_low2bp low2bp ; +p_low2bs low2bs ; p_low2high low2high ; p_low2low low2low ; p_matrix_create matrix_create ; @@ -269,6 +270,7 @@ void* dspl_load() LOAD_FUNC(log_cmplx); LOAD_FUNC(logspace); LOAD_FUNC(low2bp); + LOAD_FUNC(low2bs); LOAD_FUNC(low2high); LOAD_FUNC(low2low); LOAD_FUNC(matrix_create); diff --git a/release/include/dspl.h b/release/include/dspl.h index 4e4f814..190eaf0 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -660,6 +660,15 @@ DECLARE_FUNC(int, low2bp, double* b COMMA double* beta COMMA double* alpha); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, low2bs, double* b + COMMA double* a + COMMA int ord + COMMA double w0 + COMMA double wsl + COMMA double wsh + COMMA double* beta + COMMA double* alpha); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, low2high, double* b COMMA double* a COMMA int ord