kopia lustrzana https://github.com/Dsplib/libdspl-2.0
added filter_ws1 function
Sergey Bakhurin
rodzic
52a6a38c03
commit
3a5395f1ac
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,128 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#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");;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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*
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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*
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue