The full Fast Sinc Magnification (FSM) [1] procedure consists of modified forward and inverse FFT algorithms. The forward part uses decimation-in-frequency (DIF) type FFT without the bit-reversal reordering. No other alterations are needed for the forward FFT and thus the C code listing is not provided here. The reader can take any such implementation as long as the reordering section is removed.

The inverse algorithm uses modified decimation-in-time (DIT) type FFT and incorporates the array extension, zero padding, multiplication with a transfer function, and short-cuts to avoid unnecessary multiplications and additions by zero. This is the novel part and the full C code listing is provided here. The code is based on [2], but the arrays start from index zero as is common to C language style and the sign of the exponent in Fourier formulas is opposite using the definitions of the paper [1].

For 2D or 3D data interpolation, the one-dimensional FSM is repeated many times for the same array size and interpolation factor. Thus, it is best to have a data structure for keeping together some frequently used information. Tables of pre-calculated cosines and sines are also used.

#define REAL float typedef struct _fsm { unsigned int n; /* Original array size (must be power of two) */ unsigned int interp; /* Interpolation factor (must be power of two) */ double *cosines; /* Table of cosines */ double *sines; /* Table of sines */ REAL *imuls_r; /* Real multipliers, n floats OR */ REAL *imuls_c; /* Complex multipliers, 2*n floats */ } FSMRec, *FSM;

The structure above is used for keeping together the information
needed for the sinc interpolation. Member *n* is the size of the
original array to be interpolated and member *interp* is the
interpolation factor. Both of these values must be powers of
two. Arrays *cosines* and *sines* hold the
pre-calculated cosines and sines for angles 2π/(*j*+1), where
*j* is the index to the arrays. The index starts from value 0
(C language) and the size of these arrays must be
log_{2}(*n*)+log_{2}(*interp*)+2 entries both.

Either one of the last two arrays can be used for specifying a desired
transfer function on top of the sinc interpolation. One of these must
be NULL and the other one non-NULL. They contain either the real
valued multipliers (*imuls_r*) or complex multipliers
(*imuls_c*) for the Fourier transform of the transfer
function. These arrays hold *n* multipliers both in bit-reversal
order as they are used in the scrambled frequency domain. Values for
the complex multipliers are encoded in the C language array
*imuls_c* so that entries 0 and 1 contain the real and
imaginary parts for multiplier 0, respectively. Entries 2 and 3
contain the real and imaginary parts for the next multiplier (array
position 1 after bit reversal, position *n*/2 in the
original order), and so on. The total lengths of these C language
arrays are *n* for *imuls_r* and 2**n* for
*imuls_c*. The entries of the arrays must be scaled down by a
factor of 1/*n* to incorporate the corresponding scaling
for the inverse Fourier transform.

Typically, array *imuls_r* is used and set to value 1/*n* for
all entries. This means simple sinc interpolation without any extra
transfer function applied. If the interpolated values need to be
translated, array *imuls_c* is used with appropriate entries.

void fsmi(const FSMRec *fsm, REAL *data) ...

[1] | M. Seppä, "High-quality two-stage resampling for 3-D volumes in medical imaging", Medical Image Analysis, vol. 11, no. 4, pp. 346-360, 2007. |

[2] | W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in C. Cambridge: Cambridge University Press, 1988. |