Fractional Fourier shift

Back to Blog

03.03.2016 (20.03.2016)

Pick a sequence of points , and interpret them as points in . We would like to make sense of a fractional periodic shift for the sequence , which looks like this:

The sequence is denoted by red points, while its arc-length reparametrization is denoted by green points. The black curve is the continuous curve reconstructed from the green points. The moving blue circles visualize the parametrization of the black curve.

Motivation: Fourier descriptors

A Fourier descriptor is computed from the Fourier-series of an interpolating curve for by normalizing the desired geometric properties. A shape, given by its boundary, can then be searched among other shapes based on its Fourier descriptor, no matter what translation, scale, rotation, or fractional shift.

Property Normalized Fourier-series coefficient
Translation invariance
Rotation invariance
Fractional-shift invariance
Scale invariance ,

where is the index of Fourier-series coefficient with largest absolute value.

We demonstrate here that it is essential to take care of both noise-reduction and arc-length reparametrization before normalizing the fractional-shift.

Fractional shift

Shift of a periodic function

A shift of a periodic function is a translation with :

A shift is fractional if .

We compute the shift by taking the Fourier-series transform (FST) of , multiplying the Fourier-series coefficients by certain linear-phase exponentials, and then taking the inverse Fourier-series transform (IFST).

Shift of a finite sequence

We are studying a sequence for which the intermediate points are not explicitly defined.

To rotate the index of the sequence , we take the discrete Fourier transform (DFT) of , multiply the Fourier-transform coefficients by certain linear-phase exponentials, and take the inverse discrete Fourier transform (IDFT).

While a priori this process is defined only for integers, the mathematics allows us to specify a fractional shift. If we blindly attempt this, the drawn result seems not to generalize the shift geometrically.

However, a small change on the linear-phase exponentials fixes the problem. Why is this? How can a discrete sequence encode information about a continuous curve? We answer these questions by showing the following:

Reconstruction

We construct a smooth, periodic, band-limited , which interpolates the -sequence on the interval by

where is such that

and

Orthogonal basis

We define an inner-product on the interval by

which induces the -norm on by

The form an orthogonal basis for their span under this inner-product.

Relation to sinc-reconstruction

We rewrite the as

taking limits where appropriate. This reveals the interpolatory behaviour

The somewhat resembles the classical -reconstruction function for non-periodic functions:

However, -reconstruction is only applicable to functions of finite -norm on , which excludes non-zero periodic functions.

Instead, the is the analogous interpolation function to reconstruct an -bandlimited finite-L2-norm periodic function from uniformly-spaced samples; see S. Goldman, Information Theory, Prentice Hall, New York, 1953.

Finite Fourier series

We rewrite the expression for as a finite Fourier series

where

This shows that

where is the sequence zero-extended to a sequence in . This immediately provides the Fourier series of , and shows that is -bandlimited.

Fourier series coefficients by the discrete Fourier transform

We rewrite the Fourier-series coefficients as

where is such that

This shows that

Sampling by the inverse discrete Fourier transform

Let

where , , and . By definition, we have that

Depending on the value of , we now divide into two cases: oversampling , and undersampling . The distinction exists only because we want perform the sampling by IDFT for efficiency.

Oversampling

For the oversampling case ,

where

We rewrite the as

This shows that

That is, adding zero-coefficients can be used to increase the number of samples.

Undersampling

For the undersampling case , the extra coefficients sum into the lower coefficients:

where

This shows that

Fractional shift

To summarize, the previous shows that the fractional shift of by , resampled to , where , is achieved by

where is such that

Ringing and Gibb’s phenomenon

Using fractional shifting, we may now draw the reconstruction curve efficiently and in a numerically stable manner.

We then observe that exhibits ringing; is wavy in both transversal and longitudinal directions. In addition, when the sequence contains a jump, the point moves fast towards the next point, overshoots it, oscillates around it, and then dampens back to smoother behavior. This latter behavior is called Gibbs phenomenon.

We would like to eliminate the ringing. We do this by resampling the polygonal boundary defined by the sequence at unit-length intervals in time (green points). We call this (approximate) arc-length reparametrization.

Noise

When the points are corrupted by noise, the reconstruction curve starts to curl.

This is a consequence of requiring the curve to pass through all the points in order. We remove the curling by filtering the high-frequency components before arc-length reparametrization or approximation. This recovers the underlying curve, assuming the curve itself does not contain high-frequency components, such as sharp turns.

General filtering

We choose a filter function to smooth out the curve. The filtered Fourier series coefficients are then given by

where is how many times the filtering is repeated. Denote the inverse Fourier-series transform of the coefficients by . Then the filtering is equivalent to replacing the curve with

where is convolved with itself times. By eliminating the noise we lose the interpolation property — interpolates the local average of instead. This is exactly what we wanted.

Lanczos filtering

One of the most popular filters is given by the Lanczos filter

Then is given by

Lanczos filtering is equivalent to replacing with the local average

The :th order Lanczos filtering is equivalent to repeating this local averaging times. The is a non-negative -degree bump polynomial with support . Here is an example with :

Kaiser filter

Another popular choice is the Kaiser filter, defined by

where is the zeroth-order Bessel function, and is a parameter which controls the smoothness of the filter. The corresponds to the brick-wall filter (zeroing terms).

Approximation

We produce an approximation to the reconstructed curve by zeroing out high-frequency Fourier coefficients in the Fourier series. Suppose we use a three-term approximation, with terms -1, 0, +1. The reconstructed curves are then ellipses, and we have something like this:

The way Fourier series is constructed, the approximation with fewer terms is optimal with respect to the parametrization. However, we expect geometric optimality instead; the curve should pass close to the points. We recover geometric approximation by using an arc-length reparametrization:

For fractional shifting, both arc-length reparametrization and noise removal are important for eliminating ringing and retaining a good geometric approximation.

Software

You can find a real-time Python 2.7 program to visualize the fractional Fourier shift from here. The videos on this page have been captured from the program. The program requires NumPy, SciPy, and PyQt to run. Left mouse-button inserts new points. Right mouse-button clears all points.