libdspl-2.0
Digital Signal Processing Algorithm Library
Convolution and digital filtration.

This group describes functions for linear and circular convolutions, FIR filtration and IIR filtration. More...

Functions

int DSPL_API conv (double *a, int na, double *b, int nb, double *c)
 Real vectors linear convolution. More...
 
int DSPL_API conv_cmplx (complex_t *a, int na, complex_t *b, int nb, complex_t *c)
 Complex vectors linear convolution. More...
 
int DSPL_API conv_fft (double *a, int na, double *b, int nb, fft_t *pfft, int nfft, double *c)
 Real vectors fast linear convolution by using fast Fourier transform algorithms. More...
 
int DSPL_API conv_fft_cmplx (complex_t *a, int na, complex_t *b, int nb, fft_t *pfft, int nfft, complex_t *c)
 Complex vectors fast linear convolution by using fast Fourier transform algorithms. More...
 
int DSPL_API filter_iir (double *b, double *a, int ord, double *x, int n, double *y)
 Real IIR filtration. More...
 

Detailed Description

This group describes functions for linear and circular convolutions, FIR filtration and IIR filtration.

Function Documentation

◆ conv()

int conv ( double *  a,
int  na,
double *  b,
int  nb,
double *  c 
)

Real vectors linear convolution.


Function convolves two real vectors $ c = a * b$ length na and nb. The output convolution is a vector c with length equal to na + nb - 1.

Parameters
[in]aPointer to the first vector a.
Vector size is [na x 1].

[in]naSize of the first vector a.

[in]bPointer to the second vector b.
Vector size is [nb x 1].

[in]nbSize of the second vector b.

[out]cPointer to the convolution output vector $ c = a * b$.
Vector size is [na + nb - 1 x 1].
Memory must be allocated.

Returns
RES_OK if convolution is calculated successfully.
Else code error.
Note
If vectors a and b are coefficients of two polynomials, then convolution of the vectors a and b returns polynomial product coefficients.

Example:

double ar[3] = {1.0, 2.0, 3.0};
double br[4] = {3.0, -1.0, 2.0, 4.0};
double cr[6];
int n;
conv(ar, 3, br, 4, cr);
for(n = 0; n < 6; n++)
printf("cr[%d] = %5.1f\n", n, cr[n]);


Output:

cr[0] =   3.0
cr[1] =   5.0
cr[2] =   9.0
cr[3] =   5.0
cr[4] =  14.0
cr[5] =  12.0
Author
Sergey Bakhurin www.dsplib.org

Definition at line 157 of file conv.c.

Referenced by ratcompos().

◆ conv_cmplx()

int conv_cmplx ( complex_t a,
int  na,
complex_t b,
int  nb,
complex_t c 
)

Complex vectors linear convolution.


Function convolves two complex vectors $ c = a * b$ length na and nb. The output convolution is a vector c with length equal to na + nb - 1.

Parameters
[in]aPointer to the first vector a.
Vector size is [na x 1].

[in]naSize of the first vector a.

[in]bPointer to the second vector b.
Vector size is [nb x 1].

[in]nbSize of the second vector b.

[out]cPointer to the convolution output vector $ c = a * b$.
Vector size is [na + nb - 1 x 1].
Memory must be allocated.

Returns
RES_OK if convolution is calculated successfully.
Else code error.
Note
If vectors a and b are coefficients of two polynomials, then convolution of the vectors a and b returns polynomial product coefficients.

Example:

complex_t ac[3] = {{0.0, 1.0}, {1.0, 1.0}, {2.0, 2.0}};
complex_t bc[4] = {{3.0, 3.0}, {4.0, 4.0}, {5.0, 5.0}, {6.0, 6.0}};
complex_t cc[6];
int n;
conv_cmplx(ac, 3, bc, 4, cc);
for(n = 0; n < 6; n++)
printf("cc[%d] = %5.1f%+5.1fj\n", n, RE(cc[n]),IM(cc[n]));


Output:

cc[0] =  -3.0 +3.0j
cc[1] =  -4.0+10.0j
cc[2] =  -5.0+25.0j
cc[3] =  -6.0+32.0j
cc[4] =   0.0+32.0j
cc[5] =   0.0+24.0j
Author
Sergey Bakhurin www.dsplib.org

Definition at line 327 of file conv.c.

◆ conv_fft()

int conv_fft ( double *  a,
int  na,
double *  b,
int  nb,
fft_t pfft,
int  nfft,
double *  c 
)

Real vectors fast linear convolution by using fast Fourier transform algorithms.


Function convolves two real vectors $ c = a * b$ length na and nb in the frequency domain by using FFT algorithms. This approach provide high-performance convolution which increases with na and nb increasing. The output convolution is a vector c with length equal to na + nb - 1.

Parameters
[in]aPointer to the first vector a.
Vector size is [na x 1].

[in]naSize of the first vector a.

[in]bPointer to the second vector b.
Vector size is [nb x 1].

[in]nbSize of the second vector b.

[in]pfftPointer to the structure fft_t.
Function changes fft_t structure fields so fft_t must be clear before program returns.

[in]nfftFFT size.
This parameter set which FFT size will be used for overlapped frequency domain convolution.
FFT size must be more of minimal na and nb value. For example if na = 10, nb = 4 then nfft parameter must be more than 4.
[out]cPointer to the convolution output vector $ c = a * b$.
Vector size is [na + nb - 1 x 1].
Memory must be allocated.

Returns
RES_OK if convolution is calculated successfully.
Else code error.

Example:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 13
#define M 7
int main()
{
void* handle; /* DSPL handle */
handle = dspl_load(); /* Load DSPL function */
double a[N], b[M], c[N+M-1], d[N+M-1];
fft_t pfft;
int n, err;
linspace(0, N, N, DSPL_PERIODIC, a);
linspace(0, M, M, DSPL_PERIODIC, b);
memset(&pfft, 0, sizeof(fft_t));
err = conv_fft(a, N, b, M, &pfft, 16, c);
printf("conv_fft error: 0x%.8x\n", err);
err = conv(a, N, b, M, d);
printf("conv error: 0x%.8x\n", err);
/* print result */
for(n = 0; n < N+M-1; n++)
printf("c[%3d] = %9.2f d[%3d] = %9.2f\n", n, c[n], n, d[n]);
fft_free(&pfft); /* free fft structure memory */
dspl_free(handle); /* free dspl handle */
return 0;
}

Program output:

conv_fft error: 0x00000000
conv error:     0x00000000
c[  0] =     -0.00    d[  0] =      0.00
c[  1] =     -0.00    d[  1] =      0.00
c[  2] =      1.00    d[  2] =      1.00
c[  3] =      4.00    d[  3] =      4.00
c[  4] =     10.00    d[  4] =     10.00
c[  5] =     20.00    d[  5] =     20.00
c[  6] =     35.00    d[  6] =     35.00
c[  7] =     56.00    d[  7] =     56.00
c[  8] =     77.00    d[  8] =     77.00
c[  9] =     98.00    d[  9] =     98.00
c[ 10] =    119.00    d[ 10] =    119.00
c[ 11] =    140.00    d[ 11] =    140.00
c[ 12] =    161.00    d[ 12] =    161.00
c[ 13] =    182.00    d[ 13] =    182.00
c[ 14] =    190.00    d[ 14] =    190.00
c[ 15] =    184.00    d[ 15] =    184.00
c[ 16] =    163.00    d[ 16] =    163.00
c[ 17] =    126.00    d[ 17] =    126.00
c[ 18] =     72.00    d[ 18] =     72.00
Author
Sergey Bakhurin www.dsplib.org

Definition at line 543 of file conv.c.

◆ conv_fft_cmplx()

int conv_fft_cmplx ( complex_t a,
int  na,
complex_t b,
int  nb,
fft_t pfft,
int  nfft,
complex_t c 
)

Complex vectors fast linear convolution by using fast Fourier transform algorithms.


Function convolves two complex vectors $ c = a * b$ length na and nb in the frequency domain by using FFT algorithms. This approach provide high-performance convolution which increases with na and nb increasing. The output convolution is a vector c with length equal to na + nb - 1.

Parameters
[in]aPointer to the first vector a.
Vector size is [na x 1].

[in]naSize of the first vector a.

[in]bPointer to the second vector b.
Vector size is [nb x 1].

[in]nbSize of the second vector b.

[in]pfftPointer to the structure fft_t.
Function changes fft_t structure fields so fft_t must be clear before program returns.

[in]nfftFFT size.
This parameter set which FFT size will be used for overlapped frequency domain convolution.
FFT size must be more of minimal na and nb value. For example if na = 10, nb = 4 then nfft parameter must be more than 4.
[out]cPointer to the convolution output vector $ c = a * b$.
Vector size is [na + nb - 1 x 1].
Memory must be allocated.

Returns
RES_OK if convolution is calculated successfully.
Else code error.

Example:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 15
#define M 5
int main()
{
void* handle; /* DSPL handle */
handle = dspl_load(); /* Load DSPL function */
complex_t a[N], b[M], c[N+M-1], d[N+M-1];
fft_t pfft;
int n;
linspace(0, 2*N, 2*N, DSPL_PERIODIC, (double*)a);
linspace(0, 2*M, 2*M, DSPL_PERIODIC, (double*)b);
memset(&pfft, 0, sizeof(fft_t));
conv_fft_cmplx(a, N, b, M, &pfft, 8, c);
conv_cmplx(a, N, b, M, d);
/* print result */
for(n = 0; n < N+M-1; n++)
{
printf("c[%3d] = %9.2f%+9.2fj ", n, RE(c[n]), IM(c[n]));
printf("d[%3d] = %9.2f%+9.2fj \n", n, RE(d[n]), IM(d[n]));
}
fft_free(&pfft); /* free fft structure memory */
dspl_free(handle); /* free dspl handle */
return 0;
}

Program output:

c[  0] =     -1.00    -0.00j    d[  0] =     -1.00    +0.00j
c[  1] =     -6.00    +4.00j    d[  1] =     -6.00    +4.00j
c[  2] =    -15.00   +20.00j    d[  2] =    -15.00   +20.00j
c[  3] =    -28.00   +56.00j    d[  3] =    -28.00   +56.00j
c[  4] =    -45.00  +120.00j    d[  4] =    -45.00  +120.00j
c[  5] =    -55.00  +210.00j    d[  5] =    -55.00  +210.00j
c[  6] =    -65.00  +300.00j    d[  6] =    -65.00  +300.00j
c[  7] =    -75.00  +390.00j    d[  7] =    -75.00  +390.00j
c[  8] =    -85.00  +480.00j    d[  8] =    -85.00  +480.00j
c[  9] =    -95.00  +570.00j    d[  9] =    -95.00  +570.00j
c[ 10] =   -105.00  +660.00j    d[ 10] =   -105.00  +660.00j
c[ 11] =   -115.00  +750.00j    d[ 11] =   -115.00  +750.00j
c[ 12] =   -125.00  +840.00j    d[ 12] =   -125.00  +840.00j
c[ 13] =   -135.00  +930.00j    d[ 13] =   -135.00  +930.00j
c[ 14] =   -145.00 +1020.00j    d[ 14] =   -145.00 +1020.00j
c[ 15] =   -124.00 +1080.00j    d[ 15] =   -124.00 +1080.00j
c[ 16] =    -99.00 +1016.00j    d[ 16] =    -99.00 +1016.00j
c[ 17] =    -70.00  +820.00j    d[ 17] =    -70.00  +820.00j
c[ 18] =    -37.00  +484.00j    d[ 18] =    -37.00  +484.00j
Author
Sergey Bakhurin www.dsplib.org

Definition at line 745 of file conv.c.

Referenced by conv_fft().

◆ filter_iir()

int filter_iir ( double *  b,
double *  a,
int  ord,
double *  x,
int  n,
double *  y 
)

Real IIR filtration.


Function calculates real IIR filter output for real signal. The real filter contains real coefficients of the transfer function $H(z)$ numerator and denominator:

\[ H(z) = \frac{\sum_{n = 0}^{N} b_n z^{-n}} {1+{\frac{1}{a_0}}\sum_{m = 1}^{M} a_m z^{-n}}, \]

here $a_0$ cannot be equals zeros, $N=M=$ord.

Parameters
[in]bPointer to the vector $b$ of IIR filter transfer function numerator coefficients.
Vector size is [ord + 1 x 1].

[in]aPointer to the vector $a$ of IIR filter transfer function denominator coefficients.
Vector size is [ord + 1 x 1].
This pointer can be NULL if filter is FIR.

[in]ordFilter order. Number of the transfer function numerator and denominator coefficients (length of vectors b and a) is ord + 1.

[in]xPointer to the input signal vector.
Vector size is [n x 1].

[in]nSize of the input signal vector x.

[out]yPointer to the IIR filter output vector.
Vector size is [n x 1].
Memory must be allocated.

Returns
RES_OK if filter output is calculated successfully.
Else code error.
Example:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define ORD 6
#define N 2000
int main(int argc, char* argv[])
{
void* hdspl; /* DSPL handle */
void* hplot; /* GNUPLOT handle */
double b[ORD+1], a[ORD+1];
double t[N], s[N], n[N], sf[N];
random_t rnd;
int k;
int err;
/* Load DSPL function */
hdspl = dspl_load();
/* random generator init */
random_init(&rnd, RAND_TYPE_MT19937, NULL);
/* fill time vector */
linspace(0, N, N, DSPL_PERIODIC, t);
/* generate noise */
randn(n, N, 0, 1.0, &rnd);
/* input signal s = sin(2*pi*t) + n(t) */
for(k = 0; k < N; k++)
s[k] = sin(M_2PI*0.02*t[k]) + n[k];
/* IIR filter coefficients calculation */
iir(1.0, 70.0, ORD, 0.06, 0.0, DSPL_FILTER_ELLIP | DSPL_FILTER_LPF, b, a);
/* input signal filtration */
filter_iir(b, a, ORD, s, N, sf);
/* save input signal and filter output to the txt-files */
writetxt(t,s, N, "dat/s.txt");
writetxt(t,sf,N, "dat/sf.txt");
/* plotting by GNUPLOT */
gnuplot_create(argc, argv, 820, 340, "img/filter_iir_test.png", &hplot);
gnuplot_cmd(hplot, "unset key");
gnuplot_cmd(hplot, "set grid");
gnuplot_cmd(hplot, "set xlabel 'n'");
gnuplot_cmd(hplot, "set ylabel 's(n)'");
gnuplot_cmd(hplot, "set yrange [-3:3]");
gnuplot_cmd(hplot, "set multiplot layout 2,1 rowsfirst");
gnuplot_cmd(hplot, "plot 'dat/s.txt' with lines");
gnuplot_cmd(hplot, "set ylabel 's_f(n)'");
gnuplot_cmd(hplot, "plot 'dat/sf.txt' with lines");
gnuplot_cmd(hplot, "unset multiplot");
gnuplot_close(hplot);
/* free DSPL handle */
dspl_free(hdspl);
return err;
}

Input signal is $s(t) = \sin(2\pi \cdot 0.05 t) + n(t)$, here $n(t)$ white Gaussian noise with zero mean value and unit standard deviation.
Input signal is filtered by elliptic LPF order 6 and output signal and data saves in the txt-files

dat/s.txt  - input signal + noise
dat/sf.txt - filter output.

Plots:

GNUPLOT script for make plots is:

Author
Sergey Bakhurin www.dsplib.org

Definition at line 993 of file conv.c.

int DSPL_API conv_fft(double *a, int na, double *b, int nb, fft_t *pfft, int nfft, double *c)
Real vectors fast linear convolution by using fast Fourier transform algorithms.
Definition: conv.c:543
int DSPL_API filter_iir(double *b, double *a, int ord, double *x, int n, double *y)
Real IIR filtration.
Definition: conv.c:993
Definition: dspl.h:289
int DSPL_API conv(double *a, int na, double *b, int nb, double *c)
Real vectors linear convolution.
Definition: conv.c:157
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:359
void DSPL_API gnuplot_close(void *h)
Close GNUPLOT handle.
Definition: gnuplot.c:325
int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, int type, double *b, double *a)
Digital IIR filter design.
Definition: filter_iir.c:404
int DSPL_API conv_fft_cmplx(complex_t *a, int na, complex_t *b, int nb, fft_t *pfft, int nfft, complex_t *c)
Complex vectors fast linear convolution by using fast Fourier transform algorithms.
Definition: conv.c:745
void DSPL_API fft_free(fft_t *pfft)
Free fft_t structure.
Definition: fft.c:892
Fast Fourier Transform Object Data Structure.
Definition: dspl.h:227
void DSPL_API gnuplot_cmd(void *h, char *cmd)
Function sends cmd command to GNUPLOT corresponds to h handle.
Definition: gnuplot.c:390
int DSPL_API writetxt(double *x, double *y, int n, char *fn)
Save real data to the text file fn. .
Definition: inout.c:491
int DSPL_API random_init(random_t *prnd, int type, void *seed)
Pseudorandom numbers generators initialization.
Definition: randgen.c:121
int DSPL_API conv_cmplx(complex_t *a, int na, complex_t *b, int nb, complex_t *c)
Complex vectors linear convolution.
Definition: conv.c:327
double complex_t[2]
Complex data type.
Definition: dspl.h:86
int DSPL_API linspace(double x0, double x1, int n, int type, double *x)
Function fills a vector with n linearly spaced elements between x0 and x1.
Definition: array.c:1009
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:417
int DSPL_API gnuplot_create(int argc, char *argv[], int w, int h, char *fn_png, void **hplot)
Create GNUPLOT chart.
Definition: gnuplot.c:200