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.

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]	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[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 159 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[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 160 of file conv_cmplx.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        */
    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 202 of file conv_fft.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        */
    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 192 of file conv_fft_cmplx.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;
 
    /* 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 167 of file filter_iir.c.