GCVFFT: FFTbased computation of continuous convolution operations using generalized sampling theory
Introduction
The GCVFFT 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 bandlimited (using a bandpass prefilter) can result in strong artifacts (such as Gibbs oscillations or ringing).
The GCVFFT package (the GCVFFT “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 finiteduration signals (which are never bandlimited). Also, the package handles inputoutput 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 leastsquares approximations of the analog convolutions in shiftinvariant 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 inverseconvolution, that is, to optimally (in the leastsquares sense) recover the input samples given the result of the forward convolution.
Fast computations are ensured for periodic and mirror periodic signals. Aperiodic (finitelength) 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. 20122020, 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 bandlimited. 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 halfsample (HS) symmetry and the effective discrete convolution kernel has wholesample (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 leastsquares 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' 
RayleighSommerfield diffraction integral  
'bsplinederivative' 
Derivative of a Bspline 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 'bsplinederivative' . 

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' 
Zerocentered Bspline basis  'cardinalbspline' 
'causalbspline' 
Causal Bspline basis  
'cardinalbspline' 
Cardinal (interpolating) Bspline basis  
'sinc' 
Zerocentered sinc basis (bandlimited)  
'Delta' 
Diracdelta  
'CCD' 
Zerocentered box function with finite fillfactor (analogous to CCD sensors) 

s.deg1 
deg1 
Degree of basis for s.phi1='bspline','causalbspline','cardinalbspline' . 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','causalbspline','cardinalbspline' . Any integer between (including) 0 and 7. 
1 
s.T1 
(1D case) DeltaX1 
First sampling step: Any positive real number. 
[1,1] 
(2D 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 fillfactor 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 nonredundant half (quadrant in 2D case) of the halfsample (HS) symmetric input sequence c
and returns the nonredundant 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 wholesample (WS) symmetry about the origin to exploit gcvfft_mirror. See [1] for mathematical details. 
s.h_arg 
shiftX or [shiftX,shiftY] 

s.phi1 
'causalbspline' 

s.phi2 
'causalbspline' 

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 leastsquares 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 leastsquares 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
[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. 20122020, 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 = {20122020},
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 noncommercial educational, noncommercial research, and notforprofit 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 GCVFFT 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. 20122020, 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 GCVFFT 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.