fourier

README

Fourier

• Fast Fourier Transform implementation via Cooley-Tukey (Radix-2 DIT).
• Convolution engine which performs partitioned convolution in the frequency domain using the overlap-add method.
• Windowing functions for creating impulse responses. (e.g. Hann, Lanczos, etc)
• Functions for creating common types of FIR filters. (e.g. low-pass, high-pass, etc)

This library was written for use in a real-time audio context. Convolver allocates all of its buffers up-front and Forward/Inverse (FFT/IFFT) operate in-place. This is to avoid allocations in the hot-path. I've used this library to implement convolution reverb and perform various types of filtering.

Examples

Fast Fourier Transform

buf := make([]complex128, 8)
for i := range buf {
    buf[i] = complex(float64(i + 1), 0)
}
// buf [(1+0i) (2+0i) (3+0i) (4+0i) (5+0i) (6+0i) (7+0i) (8+0i)]
fourier.Forward(buf)
// buf [(36+0i) (-4+9.65685424949238i) (-4+4i) (-4+1.6568542494923797i)...
fourier.Inverse(buf)
// buf [(1+0i) (2+0i) (3+0i) (4+0i) (5+0i) (6+0i) (7+0i) (8+0i)] (+/- small error)

Convolution

var (
    ir       = []float64{1, 1}
    conv, _  = fourier.NewConvolver(8, ir)
    in       = []float64{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
    out      = make([]float64, len(in)+len(ir)-1)
)

_ = conv.Convolve(out, in, len(out))

// out [1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 16] (+/- small error)

Filtering

var (
    blockSize  = 256
    sampleRate = 44100.0
    cutoff     = 200.0
    kernel     = make([]float64, 32)
)

// Build a filter kernel that filters frequencies higher than 200Hz at 44.1kHz sampling rate.
filter.MakeLowPass(kernel, window.Lanczos, cutoff/sampleRate)

conv, _ := fourier.NewConvolver(blockSize, kernel)

conv.Convolve(out, in, len(out))

Documentation

Overview

Package fourier provides a Fast Fourier Transform implementation and Convolver that performs partioned convolution in the frequency domain.

Index

• func Forward(v []complex128) error
• func Inverse(v []complex128) error
• type Convolver
  • func NewConvolver(desiredBlockSize int, ir []float64, opts ...ConvolverOption) (*Convolver, error)
  • func (c *Convolver) Convolve(out, in []float64, numSamples int) error
  • func (c *Convolver) SetImpulseResponse(ir []float64) error
• type ConvolverOption
  • func ForChannel(channel, numChannels int) ConvolverOption

Examples

• Convolver (Chunks)
• Convolver (Simple)

Functions

func Forward

func Forward(v []complex128) error

Forward performs a forward FFT via Cooley-Tukey Radix-2 DIT. The buffer length is required to be a power of two.

func Inverse

func Inverse(v []complex128) error

Inverse performs an inverse FFT via Cooley-Tukey Radix-2 DIT. The buffer length is required to be a power of two.

Types

type Convolver

type Convolver struct {
	// contains filtered or unexported fields
}

Convolver performs partioned convolution using the overlap-add method. It is designed to convolve very long input streams with a FIR filter.

Example (Chunks)

var (
	blockSize = 8
	ir        = []float64{1, 1}
	conv, _   = NewConvolver(blockSize, ir)
	in        = []float64{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
	out       = make([]float64, len(in)+len(ir)-1)
)

for i := 0; i < len(out); i += blockSize {
	var (
		inBegin = min(i, len(in))
		inEnd   = min(i+blockSize, len(in))
		outEnd  = min(i+blockSize, len(out))

		inChunk  = in[inBegin:inEnd]
		outChunk = out[i:outEnd]
	)

	// We are deriving the input and output chunks here, but they would be
	// presented to you via the callback mechanisms in a streaming audio
	// scenario.
	_ = conv.Convolve(outChunk, inChunk, blockSize)
}

// Round to nearest integer (removing error) for pretty printing
const roundTo = 1e10
for i := range out {
	out[i] = math.Round(out[i]*roundTo) / roundTo
}

fmt.Println(out)

Output:

[1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 16]

Example (Simple)

var (
	blockSize = 8
	ir        = []float64{1, 1}
	conv, _   = NewConvolver(blockSize, ir)
	in        = []float64{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
	out       = make([]float64, len(in)+len(ir)-1)
)

// Convolve the entire input buffer. We're using the output buffer length
// here in order to capture the full convolved output of INPUT_LENGTH +
// IR_LENGTH - 1. Convolving more samples will result in the values
// eventually dropping to zero.
_ = conv.Convolve(out, in, len(out))

// Round to nearest integer (removing error) for pretty printing
const roundTo = 1e10
for i := range out {
	out[i] = math.Round(out[i]*roundTo) / roundTo
}

fmt.Println(out)

Output:

[1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 16]

func NewConvolver

func NewConvolver(desiredBlockSize int, ir []float64, opts ...ConvolverOption) (*Convolver, error)

NewConvolver returns a new Convolver.

desiredBlockSize is maximum size consumed from the input (loaded into an input segment) for each convolution operation (FFT -> multiply -> IFFT). Ideally, it's maximum number of samples you plan on processing for each call to Convolve.

The length of the impulse response is limited to 1920000 samples (20s * 96kHz). Exceeding this maximum length will result in truncation of the IR.

func (*Convolver) Convolve

func (c *Convolver) Convolve(out, in []float64, numSamples int) error

Convolve convolves an a chunk of input against the loaded impulse response.

func (*Convolver) SetImpulseResponse

func (c *Convolver) SetImpulseResponse(ir []float64) error

SetImpulseResponse sets the impulse response used in convolution.

type ConvolverOption

type ConvolverOption func(*Convolver) error

ConvolverOption is a configuration option for Convolver.

func ForChannel

func ForChannel(channel, numChannels int) ConvolverOption

ForChannel configures a Convolver to target a specific channel when the input buffer contains multiple interleaved channels.