fast Fourier transform library in go

[Matt Jibson]

I have written a fast Fourier transform library in Go. I am aware that
there are FFTW bindings, but I wanted a solution implemented purely in
Go. The library works for all input lengths, and computes lengths not
a power of 2 correctly (i.e., without just zero padding the input).
This initial version's purpose is functionality, not speed. There are
certainly many optimizations that can be made in the future. This is
my first experience writing Go and a package. Corrections and
suggestions are appreciated.

Code and installation instructions at:
https://github.com/mjibson/go-dsp

[Jan Mercl]

On Friday, December 9, 2011 8:35:01 AM UTC+1, Matt Jibson wrote:

Corrections and
suggestions are appreciated.

The code looks, from a glance, good to me, some notes anyway:

• No package comment, but more important - no godocs for the API/exported functions?
  
• The ensureFactors writes to the global/static *factors maps, and it is called from both of the exported Radix2FFT and BluesteinFFT which disqualifies them to be properly invoked concurrently. I suggest using locks or local data.

[Matt Jibson]

The purpose of the *factors maps is to compute, once, any needed
twiddle factors. I think that it is not prone to any race conditions
or other concurrent issues. If it is called concurrently, the worst
that could happen is some things are computed twice, but all
assignments can be run in any order, and the result should be the
same. Still though, locks are probably a good idea. Will do.

I'll look into package comments and godocs. Thanks for the notes.

[Jan Mercl]

On Friday, December 9, 2011 10:12:16 AM UTC+1, Matt Jibson wrote:

If it is called concurrently, the worst
that could happen is some things are computed twice, but all
assignments can be run in any order, and the result should be the
same.

I'm affraid that's not the case. Concurrent writes to a map can have any/strange/unpredictable behavior, corrupting the map's content being the least severe of them, but even that is probably never acceptable.

 

Still though, locks are probably a good idea. Will do.

Good decision.

[roger peppe]

it's nice to see this package.

a few other thoughts in addition to Jan's (please take with a pinch
of salt - i've never used FFT in anger):

- it looks like quite a few of the functions don't deserve
to be exported. the ones i'm thinking of are:
ToComplex, IsPowerOf2, NextPowerOf2, ZeroPad and ZeroPad2,
and perhaps even RadixFFT and BluesteinFFT given
that they're invoked appropriately as necessary.

- FFT_real might be better spelt "FFTFloat". In general variable
names in Go use camelCase rather than underscore separators.
e.g. input_len would conventionally be spelt inputLen.

- this code:
if _, present := n2_factors[input_len]; !present {
could more concisely be written as
if n2_factors[input_len] == nil {
because you will never add a nil map to n2_factors.

- you might want to use a sync.RWMutex to
guard the factors map.

[Matt Jibson]

All: thanks for the notes and feedback. I'll work on your suggestions.
Roger:
I agree that currently many of the functions don't deserve to be
exported based on how things are currently organized. However, the
plan for this library is to be a more general DSP library. Some of
these functions will be useful in some sort of dsp utils package that
provides general helper functions to the core methods. But until that
actually happens, I agree that it seems correct to make them not
exported.
Jonathan:
Now that the main FFT functions are written, I'm planning on running
benchmarks to see how useful goroutines would be. Many steps of the
radix-2 algorithm can be run in parallel. I'm hoping to see if
goroutines can produce any useful speedup there, or in other parts
(like when two FFTs can be calculated at the same time in the
Bluestein function).
Look more closely, the radix-2 function is Cooley-Tukey.
What is the best way to easily support float32? I would like to do it
without doing lots of copy-paste to just make two functions that have
different signatures but do the same thing. This is one thing Go
doesn't do well (unless I just totally misunderstand the whole
generics thing, which is likely). Can I use an interface that requires
support for math.SinCos and multiplication? That doesn't seen right...
Yes, I plan on adding more DSP packages. My next interest is audio resampling.
On Fri, Dec 9, 2011 at 3:53 AM, Jonathan Wills <runni...@gmail.com> wrote:
> Thanks for putting this together!
>
> I've also been wanting a native go fft package for a while and have been
> working on my own.  All I have right now is Cooley-Tukey, which it looks
> like you don't have, so if you want to grab any code from my repo feel free
> (it's a bit of a mess though):
>
> https://github.com/runningwild/go-fft
>
> Comments:
> You have a todo that says:  'All FFT functions: use goroutines'
> Does this mean you are going to try to parallelize ffts?  While in some
> cases that is a great idea it should at least be an option.  My typical use
> case for FFTs is running several FFTs on very large inputs, and having one
> running per core at a time, so having my FFTs get split up like that would
> slow me down.
>
> I think there should be a float32 package - basically the same thing just
> with float64 replaced by float32 and complex128 by complex64 (in as
> automated a way as possible).  This can be installed separately (as fft32
> for example), so that if you don't need the precision of float64 then you
> can save on space and get a pretty solid gain in speed (at least that has
> been my experience with fftw).
>
> The name of the package (go-dsp) suggests that you plan on adding more
> things beyond ffts.  What else did you have in mind?
>
>
> -Jonathan

[Kyle Lemons]

What is the best way to easily support float32?

The easiest way by far is to only support float64 :).  Float32 has surprisingly little precision anyway.  Currently you'd have to copy/paste or provide a helper function to repack everything as float64s.

As for the factors, you could potentially do a copy-on-write and avoid the locks.  In general, this requires the "global" reference being a pointer to the data itself.  I am not an expert, but I think it can be summarized thus:

{

  local := *global

  // do stuff

  if needToModify {

    local = copyData(local)

    // change local copy

    global = &copy

  }

  // do stuff

}

[Rob 'Commander' Pike]

For an FFT library, float32 is a perfectly reasonable type to support.

I suggest automating the production of the code with a little scripting magic to avoid the copy-paste problem, but however you do it you will, at least for now, have two copies of the implementation.

-rob

[Jonathan Wills]

On Fri, Dec 9, 2011 at 10:17 AM, Matt Jibson <matt....@gmail.com> wrote: 

Look more closely, the radix-2 function is Cooley-Tukey.

But if it is not a power of two it goes straight to Bluestein.  I guess this is more about optimization, but Bluestein is probably best used for large primes, and cooley-tukey for the rest.  As an example, I believe fftw uses Cooley-Tukey for all factors <= 19 and Bluestein or Rader for everything larger.

[Robert Bloomquist]

I'm not very familiar with FFT algorithms, but here's some more general suggestions.

The FFT() function doesn't do anything but add extra overhead to the computeFFT() function.  I suggest renaming computeFFT() to FFT() and deleting the original FFT() function.

In the computeFFT() function, it looks like when the length of the slice x is 1 or 0, you simply copy the slice into r and return it.  Originally I was wondering why you didn't simply return x, but I suppose since slices have reference semantics, it's desirable to return an actual copy.  In that case, you can just use the built-in copy() function.  Also, r is really only needed in the branch where lx <= 1; you can just return the values from the function calls directly rather than assigning them to r.  With this, you also don't need the last else statement.

It looks like the factor maps for both the radix-2 and the Bluestein algorithm are initialized no matter which algorithm is actually used.  This seems like a waste.  You could just as well initialize only the factors you need in each function, rather than doing extra work for stuff you don't need.  I thought about whether to make the factor maps local, but decided that this way you don't duplicate work calculating the same factors over and over.

In the ZeroPad() function, you create a slice of complex128s, copy the current values, and initialize the rest to zero.  In Go, all values are zero-initialized, so the extra initialization should not be necessary.

I've forked your copy on Github and implemented all these changes, as well as cleaned up various other things.  I also added in some package comments.  Hope you find them useful.

[Andrew Gerrand]

Great library! Please do document this.

Here's how it appears in godoc at present:

  http://gopkgdoc.appspot.com/pkg/github.com/mjibson/go-dsp/fft

This blog post describes some best practices for documenting Go code:

  http://blog.golang.org/2011/03/godoc-documenting-go-code.html

Cheers,

Andrew