From 3a5395f1ac1359f10fa9d2afda9b81fad8beeca5 Mon Sep 17 00:00:00 2001 From: Sergey Bakhurin Date: Sat, 31 Aug 2019 16:40:54 +0300 Subject: [PATCH] added filter_ws1 function --- dspl/src/filter_ft.c | 66 ++++++++++++++- dspl/src/filter_iir.c | 7 +- examples/bin/gnuplot/iir_test.plt | 35 ++++++++ examples/src/iit_test.c | 128 ++++++++++++++++++++++++++++++ include/dspl.c | 2 + include/dspl.h | 5 ++ release/include/dspl.c | 2 + release/include/dspl.h | 5 ++ 8 files changed, 248 insertions(+), 2 deletions(-) create mode 100644 examples/bin/gnuplot/iir_test.plt create mode 100644 examples/src/iit_test.c diff --git a/dspl/src/filter_ft.c b/dspl/src/filter_ft.c index 4058e57..7f8fcd6 100644 --- a/dspl/src/filter_ft.c +++ b/dspl/src/filter_ft.c @@ -24,7 +24,71 @@ #include "dspl.h" - +/****************************************************************************** + * Recalculate ws frequency to 1 rad/s for HPF and BANDSTOP filters + ******************************************************************************/ +double DSPL_API filter_ws1(int ord, double rp, double rs, int type) +{ + double es2, ep2, gs2, x, ws; + + if(ord<1 || rp < 0.0 || rs < 0.0) + return -1.0; + + es2 = pow(10.0, rs*0.1) - 1.0; + ep2 = pow(10.0, rp*0.1) - 1.0; + gs2 = 1.0 / (1.0 + es2); + + x = (1.0 - gs2) / (gs2 * ep2); + + switch( type & DSPL_FILTER_APPROX_MASK) + { + case DSPL_FILTER_BUTTER: + ws = pow(x, 0.5 / (double)ord); + break; + case DSPL_FILTER_CHEBY1: + case DSPL_FILTER_CHEBY2: + x = sqrt(x) + sqrt(x - 1.0); + x = log(x) / (double)ord; + ws = 0.5 * (exp(-x) + exp(x)); + break; + case DSPL_FILTER_ELLIP: + { + double k, k1; + complex_t y, z; + int res; + k = sqrt(ep2 / es2); + res = ellip_modulareq(rp, rs, ord, &k1); + if(res != RES_OK) + { + ws = -1.0; + break; + } + RE(z) = sqrt(x); + IM(z) = 0.0; + + res = ellip_acd_cmplx(&z, 1, k, &y); + if(res != RES_OK) + { + ws = -1.0; + break; + } + RE(y) /= (double)ord; + IM(y) /= (double)ord; + res = ellip_cd_cmplx(&y, 1, k1, &z); + if(res != RES_OK) + { + ws = -1.0; + break; + } + ws = RE(z); + break; + } + default: + ws = -1.0; + break; + } + return ws; +} /****************************************************************************** * low 2 bandpass transformation diff --git a/dspl/src/filter_iir.c b/dspl/src/filter_iir.c index 74f1460..5d154a6 100644 --- a/dspl/src/filter_iir.c +++ b/dspl/src/filter_iir.c @@ -97,7 +97,7 @@ int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, double *as = NULL; double *bt = NULL; double *at = NULL; - double wa0, wa1; + double wa0, wa1, ws; int err, ord_ap = ord; if(((type & DSPL_FILTER_TYPE_MASK) == DSPL_FILTER_LPF) || @@ -140,6 +140,8 @@ int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, break; case DSPL_FILTER_HPF: + ws = filter_ws1(ord_ap, rp, rs, type); + err = low2low( bs, as, ord_ap, 1.0, 1.0 / ws, bs, as); err = low2high(bs, as, ord_ap, 1.0, wa0, bt, at); break; @@ -149,6 +151,9 @@ int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, case DSPL_FILTER_BSTOP: // need frequency transform ws -> 1 rad/s + + ws = filter_ws1(ord_ap, rp, rs, type); + err = low2low( bs, as, ord_ap, 1.0, 1.0 / ws, bs, as); err = low2bs(bs, as, ord_ap, 1.0, wa0, wa1, bt, at); break; diff --git a/examples/bin/gnuplot/iir_test.plt b/examples/bin/gnuplot/iir_test.plt new file mode 100644 index 0000000..45da514 --- /dev/null +++ b/examples/bin/gnuplot/iir_test.plt @@ -0,0 +1,35 @@ +unset key +set grid +set xlabel " normalized frequency" + +set terminal pngcairo size 920, 840 enhanced font 'Verdana,8' +set output 'img/iir_test.png' +set ylabel "Magnitude, dB" +set yrange [-100:5] + +set multiplot layout 4,4 rowsfirst + +plot 'dat/iir_butter_lpf.txt' with lines +plot 'dat/iir_butter_hpf.txt' with lines +plot 'dat/iir_butter_bpf.txt' with lines +plot 'dat/iir_butter_bsf.txt' with lines + + +plot 'dat/iir_cheby1_lpf.txt' with lines +plot 'dat/iir_cheby1_hpf.txt' with lines +plot 'dat/iir_cheby1_bpf.txt' with lines +plot 'dat/iir_cheby1_bsf.txt' with lines + + +plot 'dat/iir_cheby2_lpf.txt' with lines +plot 'dat/iir_cheby2_hpf.txt' with lines +plot 'dat/iir_cheby2_bpf.txt' with lines +plot 'dat/iir_cheby2_bsf.txt' with lines + + +plot 'dat/iir_ellip_lpf.txt' with lines +plot 'dat/iir_ellip_hpf.txt' with lines +plot 'dat/iir_ellip_bpf.txt' with lines +plot 'dat/iir_ellip_bsf.txt' with lines + +unset multiplot \ No newline at end of file diff --git a/examples/src/iit_test.c b/examples/src/iit_test.c new file mode 100644 index 0000000..ec047f4 --- /dev/null +++ b/examples/src/iit_test.c @@ -0,0 +1,128 @@ +#include +#include +#include +#include "dspl.h" + +// Порядок фильтра +#define LPF_ORD 6 +#define HPF_ORD 6 +#define BPF_ORD 12 +#define BSF_ORD 12 +#define MAX_ORD BPF_ORD +#define RS 60.0 +#define RP 1.0 + +// размер векторов частотной характеристики фильтра +#define N 2048 + +void freq_resp_write2txt(double* b, double* a, int ord, int n, char* fn) +{ + double w[N], mag[N]; + int k; + + // вектор нормированной частоты от 0 до pi + linspace(0, M_PI, N , DSPL_PERIODIC, w); + + //частотные характеристика фильтра + filter_freq_resp(b, a, ord, w, n, DSPL_FLAG_LOGMAG, mag, NULL, NULL); + + for(k = 0; k < N; k++) + w[k] /= M_PI; + + // Сохранить характеристики для построения графиков + writetxt(w, mag, n, fn); +} + + +int main() +{ + void* handle; // DSPL handle + handle = dspl_load(); // Load DSPL function + + double a[MAX_ORD+1], b[MAX_ORD+1]; // коэффициенты H(s) + int err; + + // рассчитываем цифровой ФНЧ Баттерворта + iir(RP, RS, LPF_ORD, 0.3, 0.0, DSPL_FILTER_BUTTER | DSPL_FILTER_LPF, b, a); + freq_resp_write2txt(b, a, LPF_ORD, N, "dat/iir_butter_lpf.txt"); + + // рассчитываем цифровой ФВЧ Баттерворта + iir(RP, RS, HPF_ORD, 0.3, 0.0, DSPL_FILTER_BUTTER | DSPL_FILTER_HPF, b, a); + freq_resp_write2txt(b, a, HPF_ORD, N, "dat/iir_butter_hpf.txt"); + + // рассчитываем цифровой полосовой фильтр Баттерворта + iir(RP, RS, BPF_ORD, 0.3, 0.7, DSPL_FILTER_BUTTER | DSPL_FILTER_BPASS, b, a); + freq_resp_write2txt(b, a, BPF_ORD, N, "dat/iir_butter_bpf.txt"); + + // рассчитываем цифровой режекторный фильтр Баттерворта + iir(RP, RS, BSF_ORD, 0.3, 0.7, DSPL_FILTER_BUTTER | DSPL_FILTER_BSTOP, b, a); + freq_resp_write2txt(b, a, BSF_ORD, N, "dat/iir_butter_bsf.txt"); + + + + + // рассчитываем цифровой ФНЧ Чебышева 1 рода + iir(RP, RS, LPF_ORD, 0.3, 0.0, DSPL_FILTER_CHEBY1 | DSPL_FILTER_LPF, b, a); + freq_resp_write2txt(b, a, LPF_ORD, N, "dat/iir_cheby1_lpf.txt"); + + // рассчитываем цифровой ФВЧ Чебышева 1 рода + iir(RP, RS, HPF_ORD, 0.3, 0.0, DSPL_FILTER_CHEBY1 | DSPL_FILTER_HPF, b, a); + freq_resp_write2txt(b, a, HPF_ORD, N, "dat/iir_cheby1_hpf.txt"); + + // рассчитываем цифровой полосовой фильтр Чебышева 1 рода + iir(RP, RS, BPF_ORD, 0.3, 0.7, DSPL_FILTER_CHEBY1 | DSPL_FILTER_BPASS, b, a); + freq_resp_write2txt(b, a, BPF_ORD, N, "dat/iir_cheby1_bpf.txt"); + + // рассчитываем цифровой режекторный фильтр Чебышева 1 рода + iir(RP, RS, BSF_ORD, 0.3, 0.7, DSPL_FILTER_CHEBY1 | DSPL_FILTER_BSTOP, b, a); + freq_resp_write2txt(b, a, BSF_ORD, N, "dat/iir_cheby1_bsf.txt"); + + + + + // рассчитываем цифровой ФНЧ Чебышева 2 рода + iir(RP, RS, LPF_ORD, 0.3, 0.0, DSPL_FILTER_CHEBY2 | DSPL_FILTER_LPF, b, a); + freq_resp_write2txt(b, a, LPF_ORD, N, "dat/iir_cheby2_lpf.txt"); + + // рассчитываем цифровой ФВЧ Чебышева 2 рода + iir(RP, RS, HPF_ORD, 0.3, 0.0, DSPL_FILTER_CHEBY2 | DSPL_FILTER_HPF, b, a); + freq_resp_write2txt(b, a, HPF_ORD, N, "dat/iir_cheby2_hpf.txt"); + + // рассчитываем цифровой полосовой фильтр Чебышева 2 рода + iir(RP, RS, BPF_ORD, 0.3, 0.7, DSPL_FILTER_CHEBY2 | DSPL_FILTER_BPASS, b, a); + freq_resp_write2txt(b, a, BPF_ORD, N, "dat/iir_cheby2_bpf.txt"); + + // рассчитываем цифровой режекторный фильтр Чебышева 2 рода + iir(RP, RS, BSF_ORD, 0.3, 0.7, DSPL_FILTER_CHEBY2 | DSPL_FILTER_BSTOP, b, a); + freq_resp_write2txt(b, a, BSF_ORD, N, "dat/iir_cheby2_bsf.txt"); + + + + + // рассчитываем цифровой эллиптический ФНЧ + iir(RP, RS, LPF_ORD, 0.3, 0.0, DSPL_FILTER_ELLIP | DSPL_FILTER_LPF, b, a); + freq_resp_write2txt(b, a, LPF_ORD, N, "dat/iir_ellip_lpf.txt"); + + // рассчитываем цифровой эллиптический ФВЧ + iir(RP, RS, HPF_ORD, 0.3, 0.0, DSPL_FILTER_ELLIP | DSPL_FILTER_HPF, b, a); + freq_resp_write2txt(b, a, HPF_ORD, N, "dat/iir_ellip_hpf.txt"); + + // рассчитываем цифровой полосовой эллиптический фильтр + iir(RP, RS, BPF_ORD, 0.3, 0.7, DSPL_FILTER_ELLIP | DSPL_FILTER_BPASS, b, a); + freq_resp_write2txt(b, a, BPF_ORD, N, "dat/iir_ellip_bpf.txt"); + + // рассчитываем цифровой режекторный эллиптический фильтр + iir(RP, RS, BSF_ORD, 0.3, 0.7, DSPL_FILTER_ELLIP | DSPL_FILTER_BSTOP, b, a); + freq_resp_write2txt(b, a, BSF_ORD, N, "dat/iir_ellip_bsf.txt"); + + + + + dspl_free(handle); // free dspl handle + + // выполнить скрипт GNUPLOT для построения графиков + // по рассчитанным данным + return system("gnuplot -p gnuplot/iir_test.plt");; +} + + diff --git a/include/dspl.c b/include/dspl.c index cd1337c..056b886 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -84,6 +84,7 @@ p_fft_shift fft_shift ; p_fft_shift_cmplx fft_shift_cmplx ; p_filter_freq_resp filter_freq_resp ; p_filter_iir filter_iir ; +p_filter_ws1 filter_ws1 ; p_filter_zp2ab filter_zp2ab ; p_find_max_abs find_max_abs ; p_fir_linphase fir_linphase ; @@ -248,6 +249,7 @@ void* dspl_load() LOAD_FUNC(fft_shift_cmplx); LOAD_FUNC(filter_freq_resp); LOAD_FUNC(filter_iir); + LOAD_FUNC(filter_ws1); LOAD_FUNC(filter_zp2ab); LOAD_FUNC(find_max_abs); LOAD_FUNC(fir_linphase); diff --git a/include/dspl.h b/include/dspl.h index c6e4b58..b029ac4 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -514,6 +514,11 @@ DECLARE_FUNC(int, filter_iir, double* COMMA int COMMA double*); //------------------------------------------------------------------------------ +DECLARE_FUNC(double, filter_ws1, int ord + COMMA double rp + COMMA double rs + COMMA int type); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, filter_zp2ab, complex_t* COMMA int COMMA complex_t* diff --git a/release/include/dspl.c b/release/include/dspl.c index cd1337c..056b886 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -84,6 +84,7 @@ p_fft_shift fft_shift ; p_fft_shift_cmplx fft_shift_cmplx ; p_filter_freq_resp filter_freq_resp ; p_filter_iir filter_iir ; +p_filter_ws1 filter_ws1 ; p_filter_zp2ab filter_zp2ab ; p_find_max_abs find_max_abs ; p_fir_linphase fir_linphase ; @@ -248,6 +249,7 @@ void* dspl_load() LOAD_FUNC(fft_shift_cmplx); LOAD_FUNC(filter_freq_resp); LOAD_FUNC(filter_iir); + LOAD_FUNC(filter_ws1); LOAD_FUNC(filter_zp2ab); LOAD_FUNC(find_max_abs); LOAD_FUNC(fir_linphase); diff --git a/release/include/dspl.h b/release/include/dspl.h index c6e4b58..b029ac4 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -514,6 +514,11 @@ DECLARE_FUNC(int, filter_iir, double* COMMA int COMMA double*); //------------------------------------------------------------------------------ +DECLARE_FUNC(double, filter_ws1, int ord + COMMA double rp + COMMA double rs + COMMA int type); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, filter_zp2ab, complex_t* COMMA int COMMA complex_t*