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

## 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.

## ◆ 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] a Pointer to the first vector a. Vector size is [na x 1]. [in] na Size of the first vector a. [in] b Pointer to the second vector b. Vector size is [nb x 1]. [in] nb Size of the second vector b. [out] c Pointer 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 = {1.0, 2.0, 3.0};
double br = {3.0, -1.0, 2.0, 4.0};
double cr;
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 =   3.0
cr =   5.0
cr =   9.0
cr =   5.0
cr =  14.0
cr =  12.0


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] a Pointer to the first vector a. Vector size is [na x 1]. [in] na Size of the first vector a. [in] b Pointer to the second vector b. Vector size is [nb x 1]. [in] nb Size of the second vector b. [out] c Pointer 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 = {{0.0, 1.0}, {1.0, 1.0}, {2.0, 2.0}};
complex_t bc = {{3.0, 3.0}, {4.0, 4.0}, {5.0, 5.0}, {6.0, 6.0}};
complex_t cc;
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 =  -3.0 +3.0j
cc =  -4.0+10.0j
cc =  -5.0+25.0j
cc =  -6.0+32.0j
cc =   0.0+32.0j
cc =   0.0+24.0j


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] a Pointer to the first vector a. Vector size is [na x 1]. [in] na Size of the first vector a. [in] b Pointer to the second vector b. Vector size is [nb x 1]. [in] nb Size of the second vector b. [in] pfft Pointer to the structure fft_t. Function changes fft_t structure fields so fft_t must be clear before program returns. [in] nfft FFT 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] c Pointer 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 */
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


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] a Pointer to the first vector a. Vector size is [na x 1]. [in] na Size of the first vector a. [in] b Pointer to the second vector b. Vector size is [nb x 1]. [in] nb Size of the second vector b. [in] pfft Pointer to the structure fft_t. Function changes fft_t structure fields so fft_t must be clear before program returns. [in] nfft FFT 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] c Pointer 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 */
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


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] b Pointer to the vector $$b$$ of IIR filter transfer function numerator coefficients. Vector size is [ord + 1 x 1]. [in] a Pointer 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] ord Filter order. Number of the transfer function numerator and denominator coefficients (length of vectors b and a) is ord + 1. [in] x Pointer to the input signal vector. Vector size is [n x 1]. [in] n Size of the input signal vector x. [out] y Pointer 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;
/* 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:

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:350
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:420
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:1177
Fast Fourier Transform Object Data Structure.
Definition: dspl.h:278
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
Complex data type.
Definition: dspl.h:86
void dspl_free(void *handle)
Cleans up the previously linked DSPL-2.0 dynamic library.
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:1041
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:478
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