FSM in C

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.

Data structure and initialization

#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₂(n)+log₂(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.

Inverse transform

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