GCV-FFT: FFT-based computation of continuous convolution operations using generalized sampling theory

Introduction

The GCV-FFT package is a set of Matlab functions for the computation of analog convolutions of the form

g(x) = f ∗ h(x),

that is, convolutions between an analog signal f(x) (defind by a set of discrete samples c[k]) and an analog convolution kernel h(x), defined by its analog frequency response ĥ(ν). Functions in the package return discrete samples d[k] that approximate the analog result g(x).

While for bandlimited signals f(x) and bandlimited convolution kernels h(x), the samples of g = f ∗ h are simply given by

g(k) = d[k] = c[k] ∗ h[k],

that is, the discrete convolution between the discrete samples c[k] = f(k) of f and the discrete samples h[k]=h(k) of h, this relation no longer holds when the signal f or the convolution kernel h are not bandlimited. Implementations that ignore this point result in aliasing. Furthermore, artifically forcing the signals or operator to be band-limited (using a band-pass pre-filter) can result in strong artifacts (such as Gibbs oscillations or ringing).

The GCV-FFT package (the GCV-FFT “acronym” loosely stands for Generalized ConVolution via Fast Fourier Ttransforms) generalizes the simple discrete convolution result from classical sampling theory by relaxing the requiment for the input signal and convolution kernel to be bandlimited. For example, this allows handling finite-duration signals (which are never bandlimited). Also, the package handles input-output sampling step ratios that are rational (i.e. if the input is sampled with sampling step Δx the output can be sampled with sampling step p/qΔx, with p and q positive integers). The calculated coefficients correspond to least-squares approximations of the analog convolutions in shift-invariant function spaces (e.g. spline spaces) and includes the classical result as a special case (simply set the approximation spaces to shifted sinc functions).

The package also provides functions for inverse-convolution, that is, to optimally (in the least-squares sense) recover the input samples given the result of the forward convolution.

Fast computations are ensured for periodic and mirror periodic signals. Aperiodic (finite-length) signals require approriate padding.

For a complete description of the algorithms, refer to:

N. Chacko, M. Liebling, T. Blu, "Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography", J. Opt. Soc. Am. A, vol. 30, no10, pp. 2012-2020, 2013.

Discretization of continuous convolution operations based on generalized sampling theory; (a) A continuous convolution operation can be approximated using two suitable shift invariant (SI) spaces and can be numerically computed by a discrete convolution without aliasing, even when the signals are not band-limited; (b-c) When the ratio of the second sampling step to the first is p/q, (p and q being natural numbers), the input and output expansion coefficients are related by digital filters for both (b) the forward and (c) the inverse convolution operation; (d-e) Equivalent filters to (b-c) when the signals are defined by their discrete samples rather than expansion coefficients.

Functions Descriptions

gcvfft is a MATLAB function designed to enable the discretization of continuous convolution operations using FFTs, without aliasing, even when the signals are not band-limited. The algorithm [1] exploits an underlying discrete relation between two sets of expansion coefficients (or discrete samples) that uniquely characterize the continuous input and convolved output functions, respectively.

gcvfft_mirror is the more efficient variant of gcvfft when the input sequence has half-sample (HS) symmetry and the effective discrete convolution kernel has whole-sample (WS) symmetry.

igcvfft is the inverse counterpart to gcvfft and is designed to enable the reconstruction of the input expansion coefficients (or discrete samples) from the output by a filtering operation equivalent to least-squares computation.

igcvfft_mirror is the similar inverse counterpart to gcvfft_mirror.

Syntax

d = gcvfft(c,s)
[d,s] = gcvfft(c,s)

d = gcvfft_mirror(c,s)
[d,s] = gcvfft_mirror(c,s)

c = igcvfft(d,s)
[c,s] = igcvfft(d,s)

c = igcvfft_mirror(d,s)
[c,s] = igcvfft_mirror(d,s)

Parameter Description

d = gcvfft(c,s) computes the output coefficients using the input coefficients c and a structure array, s, that holds the operator parameters, as tabulated below.

Field	Usage	Meaning	Default Value
s.h	'FrT'	Unitary Fresnel transform	'FrT'
	'IFrT'	Inverse unitary Fresnel transform
	'RS'	Rayleigh-Sommerfield diffraction integral
	'bspline-derivative'	Derivative of a B-spline function
	'Identity'	No convolution
	'Shift'	Spatial shift
s.h_arg	tau or [tau,~]	Parameter for s.h='FrT','IFrT'. Any purely real or purely imaginary value.	[1,0]
	[lambda,z]	Arguments for s.h='RS'. [wavelength (positive and real), axial distance (real)]
	[~,~]	No required arguments for s.h='Identity' and 'bspline-derivative'.
	shiftX or [shiftX,shiftY]	Arguments for s.h='Shift'. Required shifts along X and Y (only real values). Positive entries for shiftX and shiftY shift values to right and down, respectively.
s.phi1	'bspline'	Zero-centered B-spline basis	'cardinal-bspline'
	'causal-bspline'	Causal B-spline basis
	'cardinal-bspline'	Cardinal (interpolating) B-spline basis
	'sinc'	Zero-centered sinc basis (band-limited)
	'Delta'	Dirac-delta
	'CCD'	Zero-centered box function with finite fill-factor (analogous to CCD sensors)
s.deg1	deg1	Degree of basis for s.phi1='bspline','causal-bspline','cardinal-bspline'. Any integer between (including) 0 and 7.	3
s.phi2	same as s.phi1	Basis that spans second SI space.	'Delta'
s.deg2	deg2	Degree of basis for s.phi2='bspline','causal-bspline','cardinal-bspline'. Any integer between (including) 0 and 7.	1
s.T1	(1-D case) DeltaX1	First sampling step: Any positive real number.	[1,1]
	(2-D case) [DeltaX1,DeltaY1]	First sampling step along x,y: Any positive real numbers.
s.p	p	Numerator of T2/T1 = p/q. Any natural number.	1
s.q	q	Denominator of T2/T1 = p/q. Any natural number.	1
s.gamma	gamma	Finite fill-factor for s.phi1 or s.phi2 = 'CCD'. 0 < gamma <= 1.	1
s.K	K	The no. of terms to be considered in the infinite sum for the computation of the FFT convolution kernel. See [1] for mathematical details. Any natural number.	500

[d,s] = gcvfft(c,s) returns the output coefficients d , as well as the structure array s containing the operator parameters.

[d,s] = gcvfft_mirror(c,s) accepts the non-redundant half (quadrant in 2D case) of the half-sample (HS) symmetric input sequence c and returns the non-redundant half (or quadrant) output coefficients d, as well as the structure array s containing the operator parameters. The structure array s holds exactly the same operator parameters as in gcvfft, except for the following restrictions.

Field	Invalid entries	Reason
s.h	'Shift'	The effective discrete convolution kernel should have whole-sample (WS) symmetry about the origin to exploit gcvfft_mirror. See [1] for mathematical details.
s.h_arg	shiftX or [shiftX,shiftY]
s.phi1	'causal-bspline'
s.phi2	'causal-bspline'
s.p	Values other than 1	Output approximation coefficients have HS symmetry similar to the input only when T2/T1 = 1. See [1] for mathematical details.
s.q	Values other than 1

[c,s] = igcvfft(d,s) computes the least-squares solution for the input coefficients c that would give rise to the given coefficients d, if the structure s were to be used. The input parameters in s should be similar to that in the forward model, as if d = gcvfft(c,s) were to be computed.

[c,s] = igcvfft_mirror(d,s) computes the least-squares solution for the HS symmetric input coefficients c that would give rise to the given HS symmetric coefficients d, if the structure s were to be used. The structure array s has the same restrictions as in gcvfft_mirror.

Demonstrations

1) demo1.m

(a) A box signal formed with N=4096 samples where the sampling step=10e-6 m and aperture width = 5.15 mm; (b) FrT computed using Fresnel integrals with 'lambda' = 632 nm and 'z' = 5 mm (only real values shown); (c) Inverse FrT of (b) computed using CV-FFT (conventional band-limited case) and (d) using GCV-FFT.

2) demo2.m

(a) Periodic input signal; (b) the continuous (exact) signal obtained by the convolution and its best approximation in the second SI space. Note that s.h = 'FrT', 'tau' = 2.5 in this illustration and only the corresponding real values have been plotted in (b).

3) demo3.m

(a) The input signal composed of three types of basis functions (box, linear and cubic B-splines) and its samples; (b) The FrT with 'tau' = 1 (only real values shown) and its samples subsequently used for reconstruction of the input signal; (c) The reconstructed samples and signal in the band-limited space obtained with the conventional CV-FFT method; (d)-(f) The reconstructed samples and signal in the shift-invariant spaces spanned by integer shifts of box, linear and cubic B-splines respectively. Clover leaves indicate reconstruction artifacts (eg. Gibbs oscillation) and asterisks denote perfect reconstruction.

4) demo4.m

(a) A typical CCD with finite-size detector elements, and (b) its corresponding family of 1-D basis functions. (c)FrT coefficients ('lambda' = 632e-9 m, 'z' = 1e-2 m, 'gamma' = 0.7) computed; (d) Reconstruction computed using CV-FFT (band-limited case) and (e) its associated SSIM map showing the presence of spatially aliased components; (f) perfect reconstructed result from (c) obtained using igcvfft yielding (g) the associated plain white SSIM map.

5) demo5.m

Comparison of methods to estimate the input coefficients from output, using discretized-continuous-inverse (gcvfft, s.h = 'IFrT' ) and alternatively, using discrete-inverse (igcvfft, s.h = 'FrT').

Reference

[1] N. Chacko, M. Liebling, T. Blu, "Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography", J. Opt. Soc. Am. A, vol. 30, no10, pp. 2012-2020, 2013.

BibTeX
@ARTICLE{gcvfft:2013,
author={N. Chacko and M. Liebling and T. Blu},
title={Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography},
journal={J. Opt. Soc. Am. A},
volume = {30},
issue = {10},
pages = {2012-2020},
year={2013},
month=oct
}

Download

By downloading this software from this site, you agree to the following terms and conditions.

Copyright © 2013. The Regents of the University of California (Regents). All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted for non-commercial educational, non-commercial research, and not-for-profit purposes only, provided that the following conditions are met:

• Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

• Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

•Neither the name of The Regents or the University of California nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 

If you publish results that are based on any of the functions included in the GCV-FFT distribution, we expect you to acknowledge our work by citing:

N. Chacko, M. Liebling, T. Blu, "Discretization of continuous convolution operators for accurate modeling of wave propagation in digital holography", J. Opt. Soc. Am. A, vol. 30, no10, pp. 2012-2020, 2013.

In that paper, you will also find the mathematical details concerning this algorithm.

Download: gcvfft_distribution_5930.zip.

 

Feedback

Please send bug reports and comments about GCV-FFT and this page to N. Chacko or M. Liebling, but please understand that this software is provided "AS IS" and that we cannot provide any support.

 

Copyright

Copyright © 2013 The Regents of the University of California
All Rights Reserved

This MATLAB package is maintained by N. Chacko and M. Liebling.