kopia lustrzana https://github.com/Dsplib/libdspl-2.0
low2bs and its example
Sergey Bakhurin
rodzic
77368c98b1
commit
73ccf3adbf
|
|
@ -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
|
||||
******************************************************************************/
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,61 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#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");;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue