fixed complex sqrt problem and changed iir example

pull/6/merge
Dsplib 2019-10-06 23:57:15 +03:00
rodzic 6de764bd94
commit 1f55f0acfa
8 zmienionych plików z 282 dodań i 55 usunięć

Wyświetl plik

@ -28,9 +28,15 @@ long_line_behaviour=1
long_line_column=72
[files]
current_page=1
FILE_NAME_0=15567;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Csrc%5Ccomplex.c;0;2
FILE_NAME_1=251;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Ciir_test.c;0;2
current_page=7
FILE_NAME_0=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Csrc%5Cfilter_ap.c;0;2
FILE_NAME_1=360;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Ccheby_poly1_test.c;0;2
FILE_NAME_2=3040;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Csrc%5Ccheby.c;0;2
FILE_NAME_3=507;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Ciir_test.c;0;2
FILE_NAME_4=890;None;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Cbin%5Cgnuplot%5Ciir_test.plt;0;2
FILE_NAME_5=1097;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_series%5Cc%5Cfourier_series_dirichlet_ex.c;0;2
FILE_NAME_6=2672;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_series%5Cc%5Cfourier_series_pimp_spectrum.c;0;2
FILE_NAME_7=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_series%5Cc%5Cfourier_series_rec.c;0;2
[build-menu]
NF_00_LB=_Собрать
@ -45,3 +51,11 @@ NF_02_WD=%p
EX_00_LB=_Execute
EX_00_CM=%e.exe
EX_00_WD=%p/examples/bin
[prjorg]
source_patterns=*.c;*.C;*.cpp;*.cxx;*.c++;*.cc;*.m;
header_patterns=*.h;*.H;*.hpp;*.hxx;*.h++;*.hh;
ignored_dirs_patterns=.*;CVS;
ignored_file_patterns=*.o;*.obj;*.a;*.lib;*.so;*.dll;*.lo;*.la;*.class;*.jar;*.pyc;*.mo;*.gmo;
generate_tag_prefs=0
external_dirs=

Wyświetl plik

@ -308,6 +308,128 @@ void DSPL_API fft_free(fft_t *pfft)
/*******************************************************************************
FFT magnitude for the real signal
*******************************************************************************/
int DSPL_API fft_mag(double* x, int n, fft_t* pfft,
double fs, int flag,
double* mag, double* freq)
{
int k, err = RES_OK;
complex_t *X = NULL;
if(!x || !pfft)
return ERROR_PTR;
if(n<1)
return ERROR_SIZE;
if(mag)
{
X = (complex_t*)malloc(n*sizeof(complex_t));
err = fft(x, n, pfft, X);
if(err!=RES_OK)
goto error_proc;
if(flag & DSPL_FLAG_LOGMAG)
for(k = 0; k < n; k++)
mag[k] = 10.0*log10(ABSSQR(X[k]));
else
for(k = 0; k < n; k++)
mag[k] = ABS(X[k]);
if(flag & DSPL_FLAG_FFT_SHIFT)
{
err = fft_shift(mag, n, mag);
if(err!=RES_OK)
goto error_proc;
}
}
if(freq)
{
if(flag & DSPL_FLAG_FFT_SHIFT)
if(n%2)
err = linspace(-fs*0.5 + fs*0.5/(double)n,
fs*0.5 - fs*0.5/(double)n,
n, DSPL_SYMMETRIC, freq);
else
err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
else
err = linspace(0, fs, n, DSPL_PERIODIC, freq);
}
error_proc:
if(X)
free(X);
return err;
}
/*******************************************************************************
FFT magnitude for the real signal
*******************************************************************************/
int DSPL_API fft_mag_cmplx(complex_t* x, int n, fft_t* pfft,
double fs, int flag,
double* mag, double* freq)
{
int k, err = RES_OK;
complex_t *X = NULL;
if(!x || !pfft)
return ERROR_PTR;
if(n<1)
return ERROR_SIZE;
if(mag)
{
X = (complex_t*)malloc(n*sizeof(complex_t));
err = fft_cmplx(x, n, pfft, X);
if(err!=RES_OK)
goto error_proc;
if(flag & DSPL_FLAG_LOGMAG)
for(k = 0; k < n; k++)
mag[k] = 10.0*log10(ABSSQR(X[k]));
else
for(k = 0; k < n; k++)
mag[k] = ABS(X[k]);
if(flag & DSPL_FLAG_FFT_SHIFT)
{
err = fft_shift(mag, n, mag);
if(err!=RES_OK)
goto error_proc;
}
}
if(freq)
{
if(flag & DSPL_FLAG_FFT_SHIFT)
if(n%2)
err = linspace(-fs*0.5 + fs*0.5/(double)n,
fs*0.5 - fs*0.5/(double)n,
n, DSPL_SYMMETRIC, freq);
else
err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, freq);
else
err = linspace(0, fs, n, DSPL_PERIODIC, freq);
}
error_proc:
if(X)
free(X);
return err;
}
/*******************************************************************************
FFT shifting
@ -338,11 +460,11 @@ int DSPL_API fft_shift(double* x, int n, double* y)
}
else
{
n2 = (n-1) >> 1;
n2 = (n+1) >> 1;
buf = (double*) malloc(n2*sizeof(double));
memcpy(buf, x, n2*sizeof(double));
memcpy(y, x+n2, (n2+1)*sizeof(double));
memcpy(y+n2+1, buf, n2*sizeof(double));
memcpy(y, x+n2, (n2-1)*sizeof(double));
memcpy(y+n2-1, buf, n2*sizeof(double));
free(buf);
}
return RES_OK;
@ -384,11 +506,11 @@ int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y)
}
else
{
n2 = (n-1) >> 1;
n2 = (n+1) >> 1;
buf = (complex_t*) malloc(n2*sizeof(complex_t));
memcpy(buf, x, n2*sizeof(complex_t));
memcpy(y, x+n2, (n2+1)*sizeof(complex_t));
memcpy(y+n2+1, buf, n2*sizeof(complex_t));
memcpy(y, x+n2, (n2-1)*sizeof(complex_t));
memcpy(y+n2-1, buf, n2*sizeof(complex_t));
free(buf);
}
return RES_OK;

Wyświetl plik

@ -2,31 +2,40 @@ unset key
set grid
set xlabel " normalized frequency"
set terminal pngcairo size 920, 840 enhanced font 'Verdana,8'
set terminal wxt size 920, 840 enhanced font 'Verdana,8'
set output 'img/iir_test.png'
set ylabel "Magnitude, dB"
set yrange [-100:5]
set xtics 0,1
set xtics add ("0.3" 0.3)
set xtics add ("0.7" 0.7)
set xtics add ("1" 1)
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
unset multiplot

Wyświetl plik

@ -3,124 +3,164 @@
#include <string.h>
#include "dspl.h"
// Порядок фильтра
/* Low-pass filter order */
#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
// размер векторов частотной характеристики фильтра
/* High-pass filter order */
#define HPF_ORD 6
/* band-pass filter order */
#define BPF_ORD 12
/* Band-stop filter order */
#define BSF_ORD 12
/* Maximum filter order */
#define MAX_ORD BPF_ORD
/* Stopband suppression (dB) */
#define RS 60.0
/* Pass-band maximum distortion (dB) */
#define RP 2.0
/* Frequency response vector size */
#define N 1024
/*******************************************************************************
* function calculates filter frequency response and save magnitude to
* the text file.
* params: b - pointer to the transfer fuction H(z) numerator vector
* a - pointer to the transfer fuction H(z) denominator vector
* ord - filter order
* n - number of magnitude vector size
* fn - file name
******************************************************************************/
void freq_resp_write2txt(double* b, double* a, int ord, int n, char* fn)
{
double w[N], mag[N];
double *w = NULL, *mag = NULL;
int k;
// вектор нормированной частоты от 0 до pi
linspace(0, M_PI, N , DSPL_PERIODIC, w);
w = (double*)malloc(n*sizeof(double));
mag = (double*)malloc(n*sizeof(double));
//частотные характеристика фильтра
/* Normalized frequency from 0 to pi */
linspace(0, M_PI, n , DSPL_PERIODIC, w);
/* Magnitude (dB) calculation */
filter_freq_resp(b, a, ord, w, n, DSPL_FLAG_LOGMAG, mag, NULL, NULL);
/* Frequency normaliztion from 0 to 1 */
for(k = 0; k < N; k++)
w[k] /= M_PI;
// Сохранить характеристики для построения графиков
/* Save magnitude to the txt file */
writetxt(w, mag, n, fn);
free(w);
free(mag);
}
int main()
/*******************************************************************************
* Main program
******************************************************************************/
int main(int argc, char* argv[])
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
void* handle; /* DSPL handle */
handle = dspl_load(); /* Load DSPL function */
double a[MAX_ORD+1], b[MAX_ORD+1]; // коэффициенты H(s)
int err;
/* Transfer function H(z) coeff. vectors */
double a[MAX_ORD+1], b[MAX_ORD+1];
/*--------------------------------------------------------------------------*/
// рассчитываем цифровой ФНЧ Баттерворта
/* LPF Batterworth */
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");
// рассчитываем цифровой ФВЧ Баттерворта
/* HPF Batterworth */
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");
// рассчитываем цифровой полосовой фильтр Баттерворта
/* Band-pass Batterworth */
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");
// рассчитываем цифровой режекторный фильтр Баттерворта
/* Band-stop Batterworth */
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 рода
/* LPF Chebyshev type 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 рода
/* HPF Chebyshev type 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 рода
/* Bnad-pass Chebyshev type 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 рода
/* Bnad-stop Chebyshev type 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 рода
/* LPF Chebyshev type 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 рода
/* HPF Chebyshev type 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 рода
/* Band-pass Chebyshev type 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 рода
/* Band-stop Chebyshev type 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");
/*--------------------------------------------------------------------------*/
// рассчитываем цифровой эллиптический ФНЧ
/* LPF Elliptic */
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");
// рассчитываем цифровой эллиптический ФВЧ
/* HPF Elliptic */
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");
// рассчитываем цифровой полосовой эллиптический фильтр
/* Band-pass Elliptic */
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");
// рассчитываем цифровой режекторный эллиптический фильтр
/* Band-stop Elliptic */
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
/* free dspl handle */
dspl_free(handle);
// выполнить скрипт GNUPLOT для построения графиков
// по рассчитанным данным
if(argc>1)
{
if(!strcmp(argv[1], "--noplot"))
return 0;
}
/* Run GNUPLOT for magnitudes plotting */
return system("gnuplot -p gnuplot/iir_test.plt");
}

Wyświetl plik

@ -80,6 +80,8 @@ p_fft fft ;
p_fft_cmplx fft_cmplx ;
p_fft_create fft_create ;
p_fft_free fft_free ;
p_fft_mag fft_mag ;
p_fft_mag_cmplx fft_mag_cmplx ;
p_fft_shift fft_shift ;
p_fft_shift_cmplx fft_shift_cmplx ;
p_filter_freq_resp filter_freq_resp ;
@ -245,6 +247,8 @@ void* dspl_load()
LOAD_FUNC(fft_cmplx);
LOAD_FUNC(fft_create);
LOAD_FUNC(fft_free);
LOAD_FUNC(fft_mag);
LOAD_FUNC(fft_mag_cmplx);
LOAD_FUNC(fft_shift);
LOAD_FUNC(fft_shift_cmplx);
LOAD_FUNC(filter_freq_resp);

Wyświetl plik

@ -170,6 +170,7 @@ typedef struct
#define DSPL_FLAG_ANALOG 0x00000001
#define DSPL_FLAG_LOGMAG 0x00000002
#define DSPL_FLAG_UNWRAP 0x00000004
#define DSPL_FLAG_FFT_SHIFT 0x00000008
@ -489,6 +490,22 @@ DECLARE_FUNC(int, fft_create, fft_t*
//------------------------------------------------------------------------------
DECLARE_FUNC(void, fft_free, fft_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft_mag, double* x
COMMA int n
COMMA fft_t* pfft
COMMA double fs
COMMA int flag
COMMA double* mag
COMMA double* freq);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft_mag_cmplx, complex_t* x
COMMA int n
COMMA fft_t* pfft
COMMA double fs
COMMA int flag
COMMA double* mag
COMMA double* freq);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft_shift, double*
COMMA int n
COMMA double*);

Wyświetl plik

@ -80,6 +80,8 @@ p_fft fft ;
p_fft_cmplx fft_cmplx ;
p_fft_create fft_create ;
p_fft_free fft_free ;
p_fft_mag fft_mag ;
p_fft_mag_cmplx fft_mag_cmplx ;
p_fft_shift fft_shift ;
p_fft_shift_cmplx fft_shift_cmplx ;
p_filter_freq_resp filter_freq_resp ;
@ -245,6 +247,8 @@ void* dspl_load()
LOAD_FUNC(fft_cmplx);
LOAD_FUNC(fft_create);
LOAD_FUNC(fft_free);
LOAD_FUNC(fft_mag);
LOAD_FUNC(fft_mag_cmplx);
LOAD_FUNC(fft_shift);
LOAD_FUNC(fft_shift_cmplx);
LOAD_FUNC(filter_freq_resp);

Wyświetl plik

@ -170,6 +170,7 @@ typedef struct
#define DSPL_FLAG_ANALOG 0x00000001
#define DSPL_FLAG_LOGMAG 0x00000002
#define DSPL_FLAG_UNWRAP 0x00000004
#define DSPL_FLAG_FFT_SHIFT 0x00000008
@ -489,6 +490,22 @@ DECLARE_FUNC(int, fft_create, fft_t*
//------------------------------------------------------------------------------
DECLARE_FUNC(void, fft_free, fft_t*);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft_mag, double* x
COMMA int n
COMMA fft_t* pfft
COMMA double fs
COMMA int flag
COMMA double* mag
COMMA double* freq);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft_mag_cmplx, complex_t* x
COMMA int n
COMMA fft_t* pfft
COMMA double fs
COMMA int flag
COMMA double* mag
COMMA double* freq);
//------------------------------------------------------------------------------
DECLARE_FUNC(int, fft_shift, double*
COMMA int n
COMMA double*);