A Survey of Gaussian Convolution Algorithms
|
Filtering utility functions. More...
Filtering utility functions.
Copyright (c) 2012-2013, Pascal Getreuer All rights reserved.
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... | |
|
static |
Half-sample symmetric boundary extension.
N | signal length |
n | requested sample, possibly outside {0,...,N -1} |
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.
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.
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,
with boundary extension by approximating the infinite sum
Definition at line 86 of file filter_util.c.
Compute taps of the impulse response of a causal recursive filter.
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 of the impulse response of the recursive filter
or equivalently
for n = 0, ..., N-1 where is the unit impulse. In the denominator coefficient array, element a[0] is not used.
Definition at line 47 of file filter_util.c.