Filtering utility functions. More...

Detailed Description

Filtering utility functions.

Author: Pascal Getreuer getre.nosp@m.uer@.nosp@m.cmla..nosp@m.ens-.nosp@m.cacha.nosp@m.n.fr

This program is free software: you can redistribute it and/or modify it under, at your option, the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version, or the terms of the simplified BSD license.

You should have received a copy of these licenses along with this program. If not, see http://www.gnu.org/licenses/ and http://www.opensource.org/licenses/bsd-license.html.

Definition in file filter_util.h.

#include "num.h"

Go to the source code of this file.

Functions
static long	extension (long N, long n)
	Half-sample symmetric boundary extension. More...

void	recursive_filter_impulse (num h, long N, const num b, int p, const num *a, int q)
	Compute taps of the impulse response of a causal recursive filter. More...

void	init_recursive_filter (num dest, const num src, long N, long stride, const num b, int p, const num a, int q, num sum, num tol, long max_iter)
	Initialize a causal recursive filter with boundary extension. More...

Function Documentation

static long extension	(	long	N,
		long	n
	)

static

Half-sample symmetric boundary extension.

Parameters

N	signal length
n	requested sample, possibly outside {0,...,`N`-1}

Returns: reflected sample in {0,...,N-1}

This function is used for boundary handling. Suppose that src is an array of length N, then src[extension(N, n)] evaluates the symmetric extension of src at location n.

Half-sample symmetric extension is implemented by the pseudocode

    repeat
        if n < 0, reflect n over -1/2
        if n >= N, reflect n over N - 1/2
    until 0 <= n < N

The loop is necessary as some n require multiple reflections to bring them into the domain {0,...,N-1}.

This function is used by all of the Gaussian convolution algorithms included in this work except for DCT-based convolution (where symmetric boundary handling is performed implicitly by the transform). For FIR, box, extended box, SII, and Deriche filtering, this function could be replaced to apply some other boundary extension (e.g., periodic or constant extrapolation) without any further changes. However, Alvarez-Mazorra and Vliet-Young-Verbeek are hard-coded for symmetric extension on the right boundary, and would require specific modification to change the handling on the right boundary.

A note on efficiency: This function is a computational bottleneck, as it is used within core filtering loops. As a small optimization, we encourage inlining by defining the function as static. We refrain from further optimization since this is a pedagogical implementation, and code readability is more important. Ideally, filtering routines should take advantage of algorithm-specific properties such as exploiting sequential sample locations (to update the extension cheaply) and samples that are provably in the interior (where boundary checks may omitted be entirely).

Definition at line 68 of file filter_util.h.

void init_recursive_filter	(	num *	dest,
		const num *	src,
		long	N,
		long	stride,
		const num *	b,
		int	p,
		const num *	a,
		int	q,
		num	sum,
		num	tol,
		long	max_iter
	)

Initialize a causal recursive filter with boundary extension.

Parameters

dest	destination array with size of at least q
src	input signal of length N
N	length of src
stride	the stride between successive samples of src
b	numerator coefficients
p	largest delay of the numerator
a	denominator coefficients
q	largest delay of the denominator
sum	the L^1 norm of the impulse response
tol	accuracy tolerance
max_iter	maximum number of samples to use for approximation

This routine initializes a recursive filter,

$\begin{aligned} u_n &= b_0 f_n + b_1 f_{n-1} + \cdots + b_p f_{n-p} &\quad - a_1 u_{n-1} - a_2 u_{n-2} - \cdots - a_q u_{n-q}, \end{aligned}$

with boundary extension by approximating the infinite sum

$u_m=\sum_{n=-m}^\infty h_{n+m}\Tilde{f}_{-n} \approx \sum_{n=-m}^{k-1} h_{n+m} \Tilde{f}_{-n}, \quad m = 0, \ldots, q - 1.$

Definition at line 86 of file filter_util.c.

void recursive_filter_impulse	(	num *	h,
		long	N,
		const num *	b,
		int	p,
		const num *	a,
		int	q
	)

Compute taps of the impulse response of a causal recursive filter.

Parameters

h	destination array of size N
N	number of taps to compute (beginning from n = 0)
b	numerator coefficients
p	largest delay of the numerator
a	denominator coefficients
q	largest delay of the denominator

Computes taps $h_0, \ldots, h_{N-1}$ of the impulse response of the recursive filter

$H(z) = \frac{b[0] + b[1]z^{-1} + \cdots + b[p]z^{-p}} {1 + a[1]z^{-1} + \cdots + a[q]z^{-q}},$

or equivalently

$\begin{aligned} h[n] = & b[0]\delta_n + b[1]\delta_{n-1}+\cdots + b[p]\delta_{n-p}\\ & -a[1]h[n-1]-\cdots -a[q]h[n-q], \end{aligned}$

for n = 0, ..., N-1 where $\delta$ is the unit impulse. In the denominator coefficient array, element a[0] is not used.

Definition at line 47 of file filter_util.c.

Detailed Description

Functions

Function Documentation