We’ve had a number of inquires about the CenterSpace FFT benchmarks, so I thought I would code up a few tests and run them on my machine. I’ve included our FFT performance numbers and the code that generated those numbers so you can try them on your machine. (If you don’t have NMath, you’ll need to download the eval version). I also did a comparison of 1 dimensional real DFTs, with FFTW, one of the fastest desktop FFT implementations available.

Benchmarks

These benchmarks were run on a 2.80 Ghz, Intel Core i7 CPU, with 4Gb of memory installed.

The clock resolution is 0.003 ns
1024 point, forward, real FFT required 4361.364 ns, Mflops 4069
1000 point, forward, real FFT required 5338.785 ns, Mflops 3235
4096 point, forward, real FFT required 21708.565 ns, Mflops 3924
4095 point, forward, real FFT required 43012.010 ns, Mflops 1980
1024 * 1024 point, forward, real FFT required 15.635 ms, Mflops 2324

I’m estimating the megaflop performance during the FFT using:

This is the asymptotic number of floating point operations for the radix-2 Cooley-Tukey FFT algorithm. This FFT MFlop estimate is used in a number of FFT benchmark reports and serves as a good basis for comparing algorithm efficiency.

As expected we take a performance hit for non-power of 2 lengths, but due to various optimizations for processing prime length FFT kernels (3, 5, 7 & 11), the performance hit is minimal in many cases. The 1000-point FFT has prime factors (2)(2)(2)(5)(5)(5), and the 4095-point FFT has prime factors (3)(3)(5)(7)(13), so those larger prime factors in the 4095-point FFT cost us some performance. Typically, user’s zero pad their data vectors to a power-of-two length to get optimal performance.

Side by side comparison with FFTW

FFTW claims to be the “Fastest Fourier Transform in the West”, and is a clever, high performance implementation of the discrete Fourier transform. This algorithm is shipped with all copies of MATLAB. FFTW is implemented in C and has the reputation as being one of the fastest desktop FFT algorithm.

Both the NMath FFT and the FFTW have a pre-computation setup that establishes the best algorithmic approach for the DFT at hand, before computing any FFT’s. This pre-computational phase is not included in the times below. In the case of the NMath FFT classes, this pre-computational phase in done in the class constructor; Therefore users must avoid constructing NMath FFT classes in tight loops for best performance (as shown in the benchmark code below). Below is a small side-by-side comparison between FFTW and NMath’s FFT (using the numbers from above).

 Comparison of a forward, real, out-of-place FFT. 


 FFT length
 FFTW 
 NMATH FFT 


 1024 
 4.14 μs
 4.36 μs 


 1000
 5.98 μs 
 5.33 μs 


 4096
 20.31 μs 
 21.71 μs 


 4095
 49.90 μs 
 43.01 μs 


 1024^2 
 17.16 ms 
 15.63 ms 

Clearly NMATH is very competitive with, and at times out-performs FFTW for real FFT’s of both power-of-2 length signals and otherwise. I chose 1D real signals as a test case because this is one of the most frequent use cases of our NMATH FFT library.

On a subjective scale, running a 1024-point FFT on a desktop commodity machine at around (an algorithm normalized) 4 GFlops is amazing. That means that in a real time measurement situation, users can compute 1024-point FFT’s at around 220kHz – all with just a couple of lines of code.

Happy Computing,
Paul

Benchmark Code

 public void BenchMarks()
    {
      Double numberTrials = 10000;
      Double flops;
 
      Stopwatch timer = new System.Diagnostics.Stopwatch();
      Console.WriteLine( String.Format("The clock resolution is {0:0.000} ns", Stopwatch.Frequency / 1000000000.0 ) );
 
      // Snip one - power of two
      RandGenUniform rand = new RandGenUniform();
      DoubleForward1DFFT fft = new DoubleForward1DFFT( 1024 );
      DoubleVector realsignal = new DoubleVector( 1024, rand );
 
      DoubleVector result = new DoubleVector( 1024 * 1024 );
 
      timer.Reset();
      for( int i = 0; i < numberTrials; i++ )
      {
        timer.Start();
        fft.FFT( realsignal, ref result );
        timer.Stop();
      }
      flops = (2.5 * 1024 * NMathFunctions.Log(1024)) / (((timer.ElapsedTicks / numberTrials) / Stopwatch.Frequency) * 1000000.0 );
      Console.WriteLine( String.Format( "1024 point, forward, real FFT required {0:0.000} ns, Mflops {1:0}", ( ( timer.ElapsedTicks / numberTrials ) / Stopwatch.Frequency ) * 1000000000.0, flops ) );
 
      // length 1000
      fft = new DoubleForward1DFFT( 1000 );
      realsignal = new DoubleVector( 1000, rand );
 
      timer.Reset();
      for( int i = 0; i < numberTrials; i++ )
      {
        timer.Start();
        fft.FFT( realsignal, ref result );
        timer.Stop();
      }
      flops = ( 2.5 * 1000 * NMathFunctions.Log( 1000 ) ) / ( ( ( timer.ElapsedTicks / numberTrials ) / Stopwatch.Frequency ) * 1000000.0 );
      Console.WriteLine( String.Format( "1000 point, forward, real FFT required {0:0.000} ns, Mflops {1:0}", ( ( timer.ElapsedTicks / numberTrials ) / Stopwatch.Frequency ) * 1000000000.0, flops ) );
 
      // length 4096
      fft = new DoubleForward1DFFT( 4096 );
      realsignal = new DoubleVector( 4096, rand );
 
      timer.Reset();
      for( int i = 0; i < numberTrials; i++ )
      {
        timer.Start();
        fft.FFT( realsignal, ref result );
        timer.Stop();
      }
      flops = ( 2.5 * 4096 * NMathFunctions.Log( 4096 ) ) / ( ( ( timer.ElapsedTicks / numberTrials ) / Stopwatch.Frequency ) * 1000000.0 );
      Console.WriteLine( String.Format( "4096 point, forward, real FFT required {0:0.000} ns, Mflops {1:0}", ( ( timer.ElapsedTicks / numberTrials ) / Stopwatch.Frequency ) * 1000000000.0, flops ) );
 
      // length 4095
      fft = new DoubleForward1DFFT( 4095 );
      realsignal = new DoubleVector( 4095, rand );
 
      timer.Reset();
      for( int i = 0; i < numberTrials; i++ )
      {
        timer.Start();
        fft.FFT( realsignal, ref result );
        timer.Stop();
      }
      flops = ( 2.5 * 4095 * NMathFunctions.Log( 4095 ) ) / ( ( ( timer.ElapsedTicks / numberTrials ) / Stopwatch.Frequency ) * 1000000.0 );
      Console.WriteLine( String.Format( "4095 point, forward, real FFT required {0:0.000} ns, Mflops {1:0}", ( ( timer.ElapsedTicks / numberTrials ) / Stopwatch.Frequency ) * 1000000000.0, flops ) );
 
 
      // length 1M
      fft = new DoubleForward1DFFT( 1024 * 1024 );
      realsignal = new DoubleVector( 1024 * 1024, rand );
 
      timer.Reset();
      for( int i = 0; i < 100; i++ )
      {
        timer.Start();
        fft.FFT( realsignal, ref result );
        timer.Stop();
      }
      flops = ( 2.5 * 1024 * 1024 * NMathFunctions.Log( 1024 * 1024 ) ) / ( ( ( timer.ElapsedTicks / 100.0 ) / Stopwatch.Frequency ) * 1000000.0 );
      Console.WriteLine( String.Format( "Million point (1024 * 1024), forward, real point FFT required {0:0.000} ms, Mflops {1:0}", ( ( timer.ElapsedTicks / 100.0 ) / Stopwatch.Frequency ) * 1000.0, flops ) );
 
    }

Most scientists and programmers understand the basic implementation details of their chosen math library. However, when algorithms are ported from one library to another, problems are hard to avoid. This seems to be particularly so when dealing with convolutions, correlations and the FFT – fundamental building blocks in many areas of computation. Frequently the theoretical concepts are clear, but when the bits hit the silicon, the confusion (at least for myself) starts.

To start eliminating some of this confusion, it’s important to understand two fundamental relationships between these three transforms.

The star is the correlation operator. Note that the first Fourier transform is conjugated, and that this breaks some basic symmetries in correlation that are found in convolution.

Convolution in R

If you happen to be porting code from R, the R language (v. 2.10.0) is distributed with a convolve function which actually, by default, returns the correlation. This is unfortunate and is a source of confusion for anybody porting a prototype from R to CenterSpace’s NMath, or to other math library for that matter.

Briefly, in R

kernel <- c(1, 2, 3, 1, 0, 0)
data <- c(1, 2, 3, 4, 5, 6)
convolve(kernel, data)

will result in:

[18 17 22 33 32 25]

Now this is a strange result, on a couple of fronts. First, this isn’t the convolution, and second, it isn’t the correlation either! Supposing that we want to compute the correlation between this kernel and signal, we would be expecting precisely,

[0 0 1 5 11 18 25 32 32 17 6 ].

Conflating correlation with convolution in one function is certain to cause confusion, because among many other reasons, convolution is communitive, and correlation is not. In general, for correlation,

yet for convolution,

If g or h satisfy certain symmetry properties, correlation can gain back the communitive property.

Back to our correlation example above, if we exchange the arguments kernel and data in R’s convolve() function we get a different answer.

convolve(data,kernel)

gives us,

[18 25 32 33 22 17 ].

There is part of a correlation swimming around in that vector, but the last three numbers given by R are not part of a linear correlation. Many users naturally take those 6 numbers incorrectly as the linear correlation (or worse convolution) of the kernel and data . This brings us to our next topic.

Fast Convolution

The fast Fourier transform is used to compute the convolution or correlation for performance reasons. This FFT based algorithm is often referred to as ‘fast convolution’, and is given by,

In the discrete case, when the two sequences are the same length, N , the FFT based method requires O(N log N) time, where a direct summation would require O(N*N) time.

This asymptotic runtime performance makes the FFT method the defacto standard for computing convolution. However, this is unfortunate because if the kernel is much smaller than the data, the direct summation is actually faster than using the FFT. This is not a rare special case, and is actually very common in signal filtering, wavelet transforms, and image processing applications. This also brings to light that many libraries (including R) require both inputs to be zero padded to the same length (typically a power of 2) – immediately eliminating this optimization and always forcing the use of the FFT technique.

Returning to our example above. If we remove the unnecessary padding from the kernel, and recompute the correlation, we arrive at,

Now since both the correlation (and the convolution) spread the signal data by kernel.Length() - 1 = 3 elements, most (engineering) users are interested in the correlation exclusively where the kernel fully overlaps the signal data. This windowing would then give us,

which are the first three numbers provided by R’s convolve function. The latter three numbers are the results of a circular , not a linear, correlation. This is probably not the result most engineers are looking for unless they are filtering a periodic signal or an image wrapped on a cylinder. Circular correlation wraps the data end-to-end in a continuous loop when summing, by effectively joining the first and last elements of the data array.

The circular correlation for this running example would look like the following table.

[1 2 3 4 5 6]
[1 2 3 1 - -] = 18
[- 1 2 3 1 -] = 25
[- - 1 2 3 1] = 32
[1 - - 1 2 3] = 33 (circular)
[3 1 - - 1 2] = 22 (circular)
[2 3 1 - - 1] = 17 (circular)

The top array is the data, and the arrays below represent the kernel sweeping across the data step by step.

The difference between the circular and linear correlation is restricted to the edges of the correlation where the (unpadded) kernel does not fully overlap the data. The circular and linear correlations are identical in the areas where the kernel fully overlaps the data – which in many applications is the area of interest.

Convolution & Correlation Classes in the NMath library

CenterSpace’s convolution and correlation classes rigorously and efficiently compute their respective transformation correctly, regardless of the computational technique used. This means that zero padding by the application programmer is no longer necessary, and in fact is discouraged. As is reflexively using the ‘fast convolution’ technique when direct summation is actually faster.

When a NMath convolution or correlation class is constructed, it estimates the number of MFlops needed by all competing techniques and chooses the fastest computational method. Zero padding will introduce errors into this MFlops estimation process.

Classes

The CenterSpace NMath library offers the following eight classes.

• {Double | Float}1DConvolution
• {DoubleComplex | FloatComplex}1DConvolution
• {Double | Float}1DCorrelation
• {DoubleComplex | FloatComplex}1DCorrelation

The two sets of correlation and convolution classes have completely symmetric interfaces.

Code Examples
If you are currently porting code from a system that uses the FFT ‘fast correlation’ technique, I will now outline how you would port that code to NMath.

Porting our running R-example from above to NMath, and assuming that what you need is linear correlation, the NMath code would look like:

DoubleVector kernel = new DoubleVector(1, 2, 3, 1);
DoubleVector data = new DoubleVector(1, 2, 3, 4, 5, 6);
 
Double1DCorrelation corr = new Double1DCorrelation(kernel, data.Length);
DoubleVector correlation = corr.Correlate(data);
DoubleVector corr_full_kernel_overlap = 
  corr.TrimConvolution(correlation, CorrelationBase.Windowing.FullKernelOverlap);
 
DoubleVector corr_centered = 
  corr.TrimConvolution(correlation, CorrelationBase.Windowing.CenterWindow);
 
// correlation =            [1 5 11 18 25 32 32 17 6]
// corr_centered =              [11 18 25 32 32 17]
// corr_full_kernel_overlap =      [18 25 32]

Note that the windowing method, TrimConvolution() does not copy any data. It just creates a windowed view (reference) into the underlying convolution vector. Windowing of native arrays are not supported because a copy would be required.

The CenterSpace NMath libraries currently do not support circular convolution, so if that is required due to the circular symmetry / periodicity of the data, the circular convolution or correlation must be computed using our FFT classes directly.

// Compute circular correlation via FFT's.
// Zero-padding is required here.
// Typically pad to the nearest power of 2.
double[] nhkernel = { 1, 2, 3, 1, 0, 0};      
double[] data = { 1, 2, 3, 4, 5, 6 };
 
// Build FFT classes
 // and setup the correct scaling.
 fft = new DoubleComplexForward1DFFT(nhkernel.Length);
 ifft = new DoubleComplexBackward1DFFT(nhkernel.Length);
 ifft.SetScaleFactorByLength();
 
// Build the complex vectors of the real data
DoubleComplexVector kernelz = 
new DoubleComplexVector(new DoubleVector(nhkernel), new DoubleVector(nhkernel.Length));
 
DoubleComplexVector dataz = 
  new DoubleComplexVector(new DoubleVector(data), new DoubleVector(nhkernel.Length));
 
// Compute.  The next five lines 
// implement the fast correlation algorithm.
fft.FFTInPlace(kernelz);
fft.FFTInPlace(dataz);
dataz = NMathFunctions.Conj(dataz)
DoubleComplexVector prodz = kernelz * dataz;
ifft.FFTInPlace(prod);
r = new DoubleVector(NMathFunctions.Real(prod));
// r = [18 17 22 33 32 25]

-Paul

See our FFT landing page for complete documentation and code examples.

Convolution is a fundamental operation in data smoothing and filtering, and is used in many other applications ranging from discrete wavelet transform’s to LTI system theory. NMath supports a high performance, forward scaling set of convolution classes that support both complex and real data. These classes will scale in performance in proportion to the number of processing cores – eliminating code rewrites to take advantage of new multi-core hardware upgrades.

The following four convolution classes are available in the NMath 4.0 library.

Additionally a symmetric set of correlation classes will be available.

Example

Computing a convolution is as simple as defining the convolution kernel, creating the right class object for the data type, and running the convolution.

// Create some random signal data using the
// Mersenne Twist random number generator.
RandomNumberGenerator rand = new RandGenMTwist(4230987);
DoubleVector data = new DoubleVector(500, rand);

// Create a simple averaging kernel.
DoubleVector kernel =  new DoubleVector("[ .25 .25 .25 .25 ]");

// Create the real number domain convolution class.
Double1DConvolution conv =
    new Double1DConvolution(kernel, data.Length);

// Compute the convolution.
DoubleVector smoothed_data = conv.Convolve(data); 

Optimal performance for all convolution problems

Exploiting the fundamental duality between convolution and the Fourier transform, and the O(n ln n) FFT algorithm, convolutions can be computed in O(n ln n) time between two sequences g and h.

For convolutions on very large data sets this is clearly the most time efficient algorithm, even though two forward and one backward FFT are required. However for shorter data sequences, this is slower in practice that just directly summing the convolution sum – particularly so on modern multi-core processors with large on-chip caches. Also, direct summation is often faster when the kernel is much shorter that the data, which is frequently the case in signal processing applications.

The decision machinery for choosing which technique to use for the problem at hand is done automatically when the class is constructed. This way the user is always getting the best available convolution performance without worrying about which technique to use.

-Paul

All variants of the original Cooley-Tukey O(n log n) fast Fourier transform fundamentally exploit different ways to factor the discrete Fourier summation of length N.

In modern implementations there really isn’t a single static FFT algorithm, but more a dynamic collection of FFT algorithms and tools that are cleverly collated for the Fourier transform type at hand. Major algorithmic changes occur in the underlying implementation as the length and forward domain (real or complex) of the problem vary. Sophisticated FFT implementations insulate the end-user programmer from all of this background machinery.

DFT length is fundamental to performance

The days of power-of-2-only FFT algorithms are dead. Users of modern FFT libraries should not need to worry about the large complexities involved in finding the optimal algorithm for the FFT computation at hand; the library should look at the FFT length, problem domain (real or complex), number of machine cores, and machine architecture, and find and compute with the best hybridized FFT algorithm available. However, it is still helpful to understand that your realized performance will depend fundamentally on the various factorization of the length of your FFT. Most know that the best FFT performance will be had when N is a power of 2. If this stringent length requirement cannot be met, then it is best to use a length that be factored into small primes. CenterSpace’s FFT algorithms contain optimized kernels for prime factor lengths of 2, 3, 5, 7 and 11. The table below demonstrates the FFT performance sensitivity to FFT length.

Clearly the fastest FFT’s are for lengths that can be factored into small primes (512, 510, 507, 506, 504), and especially small primes that have optimized kernels (512 and 504). The more kernel optimized primes your FFT length contains the faster it will run. This is a universal fact that all FFT implementations confront and holds true for higher dimension FFT’s as well. Slight changes in length can have a profound impact on FFT performance.

You can factor your FFT length using an online service to assess how your FFT will perform.

Multi-core Scalability

The ability to factor a particular FFT into a set independent computations makes it fundamentally suitable for parallelization. All modern desktop and many laptop computers today contain at least two processor cores and any modern math library should be exploiting this fact where possible. CenterSpace’s complex domain FFT’s (and related convolutions) are multi-core aware, and automatically expand to fully utilize the available processor cores. Small problems are run on a single core, but once the computational advantages of algorithm parallelization overcome the overhead costs of multi-core parallelization, the computation is spread across all available cores. This automatic parallelization is gained simply by using CenterSpace’s NMath class libraries. No end-user programming effort is involved.

The power of two FFT’s are so computationally efficient on modern processors that the gain between one and eight cores is only about 3 seconds on a 2^20-point FFT. However, for the non-power-of-two case we get a 4.5 times speed improvement going from one core to eight. Looked at another way, with multi-core scalability of the FFT, we suffered only a 2X loss in performance going from a 2^20 length FFT to a 2^20+1 length FFT, instead of a 10X loss in performance. In other words, the multi-core scalability of CenterSpace’s NMath FFT algorithms mitigate the performance loss in using non-power-of-2 lengths, and this simplifies the end-user programmer’s job.

-Paul

See our FFT landing page for complete documentation and code examples.

The next release of Center Space’s NMATH .NET libraries will contain high performance, multi-core aware, fast fourier transform classes. This set of classes will elegantly support all common 1D and 2D FFT computations in a robust easy to use object-oriented interface.

The following FFT classes will be available.

All classes efficiently support FFT’s of arbitrary length, with a simple interface for both in-place and out-of-place computations. Additionally, there is a parallel set of classes for single precision computation.

Example

Here is a simple example computing a 1000-point forward 1D FFT.


// Create some random signal data.
RandomNumberGenerator rand = new RandGenMTwist(427);
DoubleVector data = new DoubleVector(1000, rand);

// Compute the FFT
fft1000.FFTInPlace(data);


The FFT of Real (non-Complex) data results in a FFT signal of complex-conjugate symmetric data. For memory efficiency this is returned to the user in a packed format (making in-place computation possible). To facilitate the unpacking of this data, signal reader classes are supplied that support random-access indexers into the packed data. Continuing with the example above.

// Ask the FFT instance for the correct reader,
// passing in the FFT data.
DoubleSymmetricSignalReader reader
= fft1000.GetSignalReader(data);

// Also the entire result can be unpacked
DoubleComplex[] unpackedfft =
reader.UnpackFullToArray();


The readers are not necessary for the Complex versions of the FFT classes because FFT’s of Complex data is Complex and so no data packing is possible (for memory savings).

Packing Format Notes

As mentioned above, the Fourier transform of a real signal, results in a complex-conjugate symmetric signal. This symmetry is used by CenterSpace to pack the Fourier transform into an array which is the same size as the signal array.

The following table describes the layout of the packed complex-conjugate symmetric signal, of length N, in one dimension.

If we were to unroll the array, where each element in the array contains alternating real and complex values, for the case of N even, we would have an array of length 2*N.

The complexities of the packing in two dimensions increase substantially, and will not be recorded here. All NMath FFT users are encourage to use the readers to unwind packed results. Not only does this reduce coding complexity, if the underlying packing format changes, the readers will still provide the expected functionality.

Finally, when inverting complex-conjugate symmetric signals, using the DoubleSymmetricBackward1DFFT class, the input signals are expect be packed.

-Paul