Documentation
¶
Overview ¶
Package fft provides support for 1 dimensional fast fourier transform and spectra.
Package fft is part of http://zikichombo.org
Fast fourier transforms provide frequency domain representation of time domain signals. Package fft supports a fairly efficient fourier transform implementation with a lot of simplicity and convenience.
Features ¶
Package fft is designed primarily for repeated usage on a data stream. As such, the incremental interface supports the creation of buffers which automatically have size and capacity which guarantee no copying or allocations during executation without the user needing to worry about the details of the sizes and capacities.
Package fft provides an efficient real-only interface for even transform sizes with spectra represented in a half complex format like fftw. An odd transform size real-only interface is also supported, but is not as efficient as the even transform size real-only interface.
The interface guarantees the Parseval equation, which states the sum of squares of the amplitudes in the time domain equals the sum of squares of the frequency coeficients in the frequency domain. It also guarantees that the inverse transform is an inverse without the user needing to worry about scaling.
Non Features ¶
Package fft is designed exclusively for 1d data. Support for matrices is a non-goal of this package.
Algorithms ¶
The implementation uses in place decimation in time binary radix 2 Cooley Tuckey for sizes of powers of 2, and Bluestein algorithm otherwise. Twiddle factors are pre-computed, and attention is paid to minimize allocations and copying, both internally and in the interface. Package fft provides O(N Log N) transforms for all inputs and allows for in place transforms. The Real only interface uses the complex interface for half-sized inputs together with some O(N) pre/post processing.
Index ¶
- func BinRange(sFreq freq.T, n, b int) (l, u freq.T)
- func Dilate(d []complex128, p, q int)
- func Do(d []complex128)
- func FreqAt(sFreq freq.T, n int, i float64) freq.T
- func FreqBin(sFreq, aFreq freq.T, n int) int
- func Inv(d []complex128)
- func InvTo(dst, src []complex128) []complex128
- func Ny(n int) int
- func R2Size(n int) int
- func To(dst, src []complex128) []complex128
- func WinSize(fS, fA, fE freq.T) (n, c int)
- func ZeroPad(d []complex128, n int) []complex128
- func ZeroPadTo(dst, src []complex128, n int) []complex128
- type HalfComplex
- func (h HalfComplex) Cmplx(i int) complex128
- func (h HalfComplex) FromCmplx(d []complex128)
- func (h HalfComplex) FromPolar(mag, ph []float64)
- func (h HalfComplex) Imag(i int) float64
- func (h HalfComplex) Len() int
- func (a HalfComplex) MulElems(b HalfComplex) HalfComplex
- func (h HalfComplex) Real(i int) float64
- func (h HalfComplex) SetCmplx(i int, c complex128)
- func (h HalfComplex) SetImag(i int, v float64)
- func (h HalfComplex) SetReal(i int, v float64)
- func (h HalfComplex) ToCmplx(d []complex128)
- func (h HalfComplex) ToPolar(mag, ph []float64)
- type Real
- type S
- func (s *S) At(i int) complex128
- func (s *S) CopyFrom(t *S)
- func (s *S) FoldReal()
- func (s *S) FromHalfComplex(hc HalfComplex) error
- func (s *S) FromRect(d []complex128) error
- func (s *S) ItpPeaks(dst []float64) []float64
- func (s *S) ItpQMag(f float64) float64
- func (s *S) Mag(i int) float64
- func (s *S) MagDb(i int) float64
- func (s *S) N() int
- func (s *S) Ny() int
- func (s *S) PeakItpQ(i int) (idx float64, mag float64, phase float64)
- func (s *S) Peaks() []int
- func (s *S) PeaksTo(dst []int) []int
- func (s *S) Phase(i int) float64
- func (s *S) PlotMag(b image.Rectangle) *image.RGBA
- func (s *S) PlotMagTo(b image.Rectangle, p string) error
- func (s *S) Power() float64
- func (s *S) Rect(dst []complex128) []complex128
- func (s *S) SetMag(i int, m float64)
- func (s *S) SetPhase(i int, p float64)
- func (s *S) ToHalfComplex(dst HalfComplex) HalfComplex
- type T
- func (t *T) AutoCorr(d []complex128) error
- func (t *T) Cap() int
- func (t *T) Do(d []complex128) error
- func (t *T) Inv(d []complex128) error
- func (t *T) InvTo(dst, src []complex128) ([]complex128, error)
- func (t *T) N() int
- func (t *T) Scale(v bool)
- func (t *T) To(dst, src []complex128) ([]complex128, error)
- func (t *T) Win(c []complex128) []complex128
Examples ¶
Constants ¶
This section is empty.
Variables ¶
This section is empty.
Functions ¶
func BinRange ¶
BinRange gives the range of frequencies in a DFT frequency bin b where the DFT is based on an n-sample window at rate sFreq.
func Dilate ¶
func Dilate(d []complex128, p, q int)
Dilate changes the frequency basis by a factor of n/m. This is a pitch shift relative to the quantized frequency domain, hence there is informaion loss, and some values may just be clobbered...
func Do ¶
func Do(d []complex128)
Do performs an in-place FFT on d, if possible.
As fft implementations require different slice lengths and capacities depending on the length of d, the operation Do will only be in-place if d has len and cap equal to
New(len(d)).Win(d)
func FreqAt ¶
FreqAt gives the frequency at a floating point index i for A FT for data of length n with sample frequency sFreq.
func FreqBin ¶
FreqBin gives the DFT frequency bin of frequency a in a window of n-sample at rate sFreq
func InvTo ¶
func InvTo(dst, src []complex128) []complex128
InvTo is like To but for inverse transform
func Ny ¶
Ny gives the index of the first frequency bin at or above the Nyquist limit of a FT of size n.
func R2Size ¶
R2Size returns the least power of 2 not less than n, the size for which fourier transforms use the fast radix-2 algorithm.
func To ¶
func To(dst, src []complex128) []complex128
To performs a FFT on src placing the results in dst if dst has approprioate capacity returning either dst or the results in a new slice.
func WinSize ¶
WinSize gives the smallest positive window size (n, c) (n samples covering c cycles of a) of a DFT applied to input signal with sample rate fS for frequency fA such that a cycle appears after n samples. The cycle is approximate for period functions of frequency fA but exact for some frequency within an error range around that (at worst within [fA .. fA + fE]).
Special case:
if it overflows, it returns -1, -1
For example, to find a good window size for detecting a=440 Hz in an s=17.6Khz sample, we'd call WinSize(440Hz, 17.6KHz, 0.001Hz), giving n, and c. The window size n in any signal of sampling rate s has c full cycles of any signal at frequency a, with very little error because (sc % a) is small. As a result, using a DFT over n samples will often have less error and noise in it, since the DFT assumes a cyclic signal (actually stationary, cyclic applied infinitely satisfies).
func ZeroPad ¶
func ZeroPad(d []complex128, n int) []complex128
ZeroPad returns a zero-padding of the Fourier Transform d for time interpolation. If we divide a transform into a 0(dc) component, components not exceeding the Nyquist limit, and higher frequencies as "negative" frequencies, then the zero padding goes in-between the non-negative frequencies and the negative frequencies. The inverse FT will then generate time-interpolation of the data for which d is a FT.
func ZeroPadTo ¶
func ZeroPadTo(dst, src []complex128, n int) []complex128
ZeroPadTo is like ZeroPad but allows specifying a destination vector. If dst doesn't have sufficient capacity, then a new one is created and returned.
Types ¶
type HalfComplex ¶
type HalfComplex []float64
HalfComplex is a format for storing complex spectra of real data of length N in a []float64 of length N. Such spectra have Hermitian symmetry.
Since the spectrum is so symmetric, all the information will fit, provided we reflect around the points correctly.
The format is (like in fftw):
r0, r1, ..., r_{n/2}, i_{(n+1)/2 - 1}, ..., i2, i1
for complex values rn + in, n an integer in 0...N/2
Some things to note
Due to the symmetry, i0 is always 0 and i_{n/2} is always 0 if N is even.
The number of reals is 2 greater than the number of imaginaries when N is even
The number of reals is 1 greater than the number of imaginaries when N is odd because i_{n/2} doesn't exist.
func (HalfComplex) Cmplx ¶
func (h HalfComplex) Cmplx(i int) complex128
Cmplx returns the complex128 representation of element i.
func (HalfComplex) FromCmplx ¶
func (h HalfComplex) FromCmplx(d []complex128)
FromCmplx places a complex spectrum of a real sequence in d.
FromCmplx panics if len(h) != len(d).
FromCmplx does not check that d is in the symmetric form of a DFT of real data.
func (HalfComplex) FromPolar ¶
func (h HalfComplex) FromPolar(mag, ph []float64)
FromPolar places a complex spectrum of a real sequence in h.
FromCmplx panics if len(h) != len(d).
FromCmplx does not check that d is in the symmetric form of a DFT of real data.
func (HalfComplex) Imag ¶
func (h HalfComplex) Imag(i int) float64
Imag returns the imaginary part of the complex number at i.
func (HalfComplex) Len ¶
func (h HalfComplex) Len() int
Len returns the number of complex numbers in h, which does not include elements with symmetric pairs.
func (HalfComplex) MulElems ¶
func (a HalfComplex) MulElems(b HalfComplex) HalfComplex
MulElems computes the elementwise multiplication of a and b, placing the result in a and returning it. MulElems panics if a.Len() != b.Len().
func (HalfComplex) Real ¶
func (h HalfComplex) Real(i int) float64
Real returns the real part of the complex number at i.
func (HalfComplex) SetCmplx ¶
func (h HalfComplex) SetCmplx(i int, c complex128)
SetCmplx sets the complex number i to c in h.
func (HalfComplex) SetImag ¶
func (h HalfComplex) SetImag(i int, v float64)
SetImag sets the imaginary part of the complex number at i to v. Since all imaginary parts at complex number 0 and h.Len()/2 are 0, if i == 0 or h.Len()/2, then SetImag is a no-op.
func (HalfComplex) SetReal ¶
func (h HalfComplex) SetReal(i int, v float64)
SetReal sets the real part of the complex number at i.
func (HalfComplex) ToCmplx ¶
func (h HalfComplex) ToCmplx(d []complex128)
ToCmplx fills d with complex numbers stored in h.
if len(h) != len(d), then ToCmplx panics.
func (HalfComplex) ToPolar ¶
func (h HalfComplex) ToPolar(mag, ph []float64)
ToPolar fills mag, phase with the magnitude and phase of complex numbers stored in h
type Real ¶
type Real struct {
// contains filtered or unexported fields
}
Real computes an FFT for a Real only data.
For even length transforms, the implementation uses a complex FFT of size N/2 and some pre/post processing. For odd length transforms, the implementtion uses a complex FFT of size N.
func (*Real) Do ¶
func (r *Real) Do(d []float64) HalfComplex
Do performs a DFT on real data in d.
d must be the size specified in the call to NewReal() which created r, or Do will panic.
Do operates in place, overwriting d. Do returns d overwritten as (i.e. cast to) a HalfComplex object.
func (*Real) Inv ¶
func (r *Real) Inv(hc HalfComplex) []float64
Inv computes the inverse transform of a real data from a HalfComplex object.
Inv operates in place but returns the same data as hc, cast to a []float64.
type S ¶
type S struct {
// contains filtered or unexported fields
}
S provides convenience wrappers around a ft spectrum.
func NewS ¶
func NewS(d []complex128) *S
NewS creates a spectrum from the data in d. Once NewS returns, d is free to be used or gc'd.
func NewSHalfComplex ¶
func NewSHalfComplex(hc HalfComplex) *S
NewSHalfComplex creates a new spectrum object from a HalfComplex object.
func (*S) At ¶
func (s *S) At(i int) complex128
At returns the complex representing the spectrum value at symmetric index i.
func (*S) FoldReal ¶
func (s *S) FoldReal()
FoldReal takes the spectrum and makes all negative frequency bins complex conjugates of the corresponding positive frequency bin to guarantee real output of inverse fft.
func (*S) FromHalfComplex ¶
func (s *S) FromHalfComplex(hc HalfComplex) error
FromHalfComplex makes s contain spectrum from hc.
FromHalfComplex returns a non-nil error if s doesn't contain the same number of elements as hc.
func (*S) FromRect ¶
func (s *S) FromRect(d []complex128) error
FromRect resets s to use the complex spectrum d.
func (*S) ItpPeaks ¶
ItpPeaks interpolates all the peaks in the spectrum, returning their interpolated indices, magnitudes and phases. Quadratic interpolation on log scale magnitudes is used, as in PeakItpQ, but the returned magnitudes are not log scale.
The returned slice contains interpolated indices at n*3 positions and interpolated magnitudes at n*3+1 positions, and interpolated phases at n*3+2 position.
For example, if there are two peaks at 1.23 and 30.91 with magnitudes 10 and 100, and phases pi/49, 8pi/9 then the returned slice would be
{1.23, 10, pi/49, 30.91, 100, 8pi/9}
func (*S) ItpQMag ¶
ItpQMag uses quadratic interpolation to find the magnitude at index i. Linear interpolation is used if s.N() == 2. ItpQMag panics if s.N() < 2.
func (*S) Ny ¶
Ny returns the index of the first bin at or above the Nyquist from the index of the input from NewS().
func (*S) PeakItpQ ¶
PeakItpQ performs interpolation of spectrum peaks, giving a floating point index, magnitude, and phase. To retrieve the frequency at the index, FreqAt is available. Peaks often can correspond to sinusoidal waves which are off-center of the frequency bin. The peak shape is modelled as a parabola from neighboring points.
If i is <= 1 or at the end of the Nyquist limit, then no interpolation takes place and the bin information is returned.
func (*S) Peaks ¶
Peaks returns the indices of the non-negative frequency bins which are higher than one of their two neighbors and not less than either neighbor. If there is only one element, that element is returned. Endpoints are treated as though they are strictly higher than beyond the endpoint.
func (*S) Power ¶
Power returns the total power of the spectrum, the sum of squares of magnitudes. Power assumes s represents real data.
func (*S) Rect ¶
func (s *S) Rect(dst []complex128) []complex128
Rect puts the spectrum in rectangular complex (real + imag) form in dst. If dst doesn't have capacity for the data, then a new slice is allocated and returned. Otherwise, the results are placed in dst and returned.
func (*S) ToHalfComplex ¶
func (s *S) ToHalfComplex(dst HalfComplex) HalfComplex
ToHalfComplex places the spectrum s in dst.
If dst doesn't have capacity for the data, then a new slice is allocated and returned. Otherwise, the results are placed in dst and returned.
type T ¶
type T struct {
// contains filtered or unexported fields
}
T maintains state for efficient repeated computation on windows of data.
Example ¶
var d = []complex128{(1 + 0i), (1 + 0i), (1 + 0i)}
tr := New(len(d))
w := tr.Win(nil)
copy(w, d) // for correct capacity.
tr.Do(w)
tr.Inv(w)
Output:
func New ¶
New creates a new T which can compute repeated windows more efficiently than repeatedly calling Do or To
- n is the size of the data windows to transform
- inv is whether or not to use the "inverse" transform
func (*T) AutoCorr ¶
func (t *T) AutoCorr(d []complex128) error
AutoCorr computes the (circular) autocorrelation of the input d, returning an error if d isn't proper size/capacity as in T.Do(). To obtain normal (actually tapered) autocorrelation, zero padding must be added.
func (*T) Inv ¶
func (t *T) Inv(d []complex128) error
Inv performs an in-place inverse transform on d.
func (*T) InvTo ¶
func (t *T) InvTo(dst, src []complex128) ([]complex128, error)
InvTo perform an inverse FFT on src, placing the results in dst and leaving src untouched. To returns dst, e
A non-nil error can occur if src and dst are not aligned with t. A new dst is allocated and returned if dst is nil.
func (*T) Scale ¶
Scale can be set to turn on or off scaling. By default, scaling is on.
Scaling in this fft implementation divides the results of both forward and backwards transforms by sqrt(t.N()), which enforces Parsevals power equivalence between powers input and output to transforms and enforces that an inverse is an inverse, rather than inverse on a different scale.
func (*T) To ¶
func (t *T) To(dst, src []complex128) ([]complex128, error)
To perform a FFT on src, placing the results in dst and leaving src untouched. To returns dst, e
A non-nil error can occur if src and dst are not aligned with t. A new dst is allocated and returned if dst is nil.
func (*T) Win ¶
func (t *T) Win(c []complex128) []complex128
Win returns a slice of data with dimensions set so no error occurs if used as argument to Do() or Inv(), nor if used as destination in To() or InvTo().
These errors ensure not only that the input size makes sense but also the capacity, which might vary from length due to internal zero padding.
Returned windows are of length t.N() and capacity t.Cap() and contain all the data in c
Source Files