# 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) = fh(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 = fh 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. ## 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  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  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  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  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` 2) `demo2.m` 3) `demo3.m` 4) `demo4.m` 5) `demo5.m` ## Reference

 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
}

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.

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