fourier/internal/fftpack: factor out FFTPACK port to internal package
MODIFY
fourier/fourier.go
fourier/fourier.go
package fourier
+import "gonum.org/v1/gonum/fourier/internal/fftpack"
+
// FFT implements Fast Fourier Transform and its inverse for real sequences.
type FFT struct {
work []float64
@@ -35,7 +37,7 @@ func (t *FFT) Reset(n int) {
t.work = make([]float64, 2*n)
t.real = make([]float64, n)
}
- rffti(n, t.work, t.ifac[:])
+ fftpack.Rffti(n, t.work, t.ifac[:])
}
// FFT computes the Fourier coefficients of the input sequence, seq,
@@ -56,7 +58,7 @@ func (t *FFT) FFT(dst []complex128, seq []float64) []complex128 {
panic("fourier: destination length mismatch")
}
copy(t.real, seq)
- rfftf(len(t.real), t.real, t.work, t.ifac[:])
+ fftpack.Rfftf(len(t.real), t.real, t.work, t.ifac[:])
dst[0] = complex(t.real[0], 0)
if len(seq) < 2 {
return dst
@@ -105,7 +107,7 @@ func (t *FFT) IFFT(dst []float64, coeff []complex128) []float64 {
dst[2*i+1] = real(cv)
dst[2*i+2] = imag(cv)
}
- rfftb(len(dst), dst, t.work, t.ifac[:])
+ fftpack.Rfftb(len(dst), dst, t.work, t.ifac[:])
return dst
}
@@ -150,7 +152,7 @@ func (t *CmplxFFT) Reset(n int) {
t.work = make([]float64, 4*n)
t.real = make([]float64, 2*n)
}
- cffti(n, t.work, t.ifac[:])
+ fftpack.Cffti(n, t.work, t.ifac[:])
}
// FFT computes the Fourier coefficients of a complex input sequence,
@@ -175,7 +177,7 @@ func (t *CmplxFFT) FFT(dst, seq []complex128) []complex128 {
t.real[2*i] = real(cv)
t.real[2*i+1] = imag(cv)
}
- cfftf(len(dst), t.real, t.work, t.ifac[:])
+ fftpack.Cfftf(len(dst), t.real, t.work, t.ifac[:])
for i := range dst {
dst[i] = complex(t.real[2*i], t.real[2*i+1])
}
@@ -204,7 +206,7 @@ func (t *CmplxFFT) IFFT(dst, coeff []complex128) []complex128 {
t.real[2*i] = real(cv)
t.real[2*i+1] = imag(cv)
}
- cfftb(len(dst), t.real, t.work, t.ifac[:])
+ fftpack.Cfftb(len(dst), t.real, t.work, t.ifac[:])
for i := range dst {
dst[i] = complex(t.real[2*i], t.real[2*i+1])
}

RENAME
fourier/array_bounds_checks.go
fourier/internal/fftpack/array_bounds_checks.go
// +build bounds
-package fourier
+package fftpack
import "fmt"

RENAME
fourier/array_no_bounds_checks.go
fourier/internal/fftpack/array_no_bounds_checks.go
// +build !bounds
-package fourier
+package fftpack
// The types in array.go implement Fortran-like arrays for bootstrapping
// the implementation of the FFT functions translated from FFTPACK; they

RENAME
fourier/cfft.go
fourier/internal/fftpack/cfft.go
// Paul N Swarztrauber, placed in the public domain at
// http://www.netlib.org/fftpack/.
-package fourier
+package fftpack
import "math"
-// cffti initializes the array work which is used in both cfftf
-// and cfftb. the prime factorization of n together with a
+// Cffti initializes the array work which is used in both Cfftf
+// and Cfftb. the prime factorization of n together with a
// tabulation of the trigonometric functions are computed and
// stored in work.
//
@@ -22,14 +22,14 @@ import "math"
// Output parameters:
//
// work A work array which must be dimensioned at least 4*n.
-// the same work array can be used for both cfftf and cfftb
+// the same work array can be used for both Cfftf and Cfftb
// as long as n remains unchanged. Different work arrays
// are required for different values of n. The contents of
-// work must not be changed between calls of cfftf or cfftb.
+// work must not be changed between calls of Cfftf or Cfftb.
//
// ifac A work array containing the factors of n. ifac must have
// length 15.
-func cffti(n int, work []float64, ifac []int) {
+func Cffti(n int, work []float64, ifac []int) {
if len(work) < 4*n {
panic("fourier: short work")
}
@@ -113,8 +113,8 @@ outer:
}
}
-// cfftf computes the forward complex Discrete Fourier transform
-// (the Fourier analysis). Equivalently, cfftf computes the
+// Cfftf computes the forward complex Discrete Fourier transform
+// (the Fourier analysis). Equivalently, Cfftf computes the
// Fourier coefficients of a complex periodic sequence. The
// transform is defined below at output parameter c.
//
@@ -128,13 +128,13 @@ outer:
// to be transformed.
//
// work A real work array which must be dimensioned at least 4*n.
-// in the program that calls cfftf. The work array must be
-// initialized by calling subroutine cffti(n,work,ifac) and a
+// in the program that calls Cfftf. The work array must be
+// initialized by calling subroutine Cffti(n,work,ifac) and a
// different work array must be used for each different
// value of n. This initialization does not have to be
// repeated so long as n remains unchanged thus subsequent
// transforms can be obtained faster than the first.
-// the same work array can be used by cfftf and cfftb.
+// the same work array can be used by Cfftf and Cfftb.
//
// ifac A work array containing the factors of n. ifac must have
// length of at least 15.
@@ -147,18 +147,18 @@ outer:
//
// where i=sqrt(-1)
//
-// This transform is unnormalized since a call of cfftf
-// followed by a call of cfftb will multiply the input
+// This transform is unnormalized since a call of Cfftf
+// followed by a call of Cfftb will multiply the input
// sequence by n.
//
// The n elements of c are represented in n pairs of real
// values in r where c[j] = r[j*2]+r[j*2+1]i.
//
// work Contains results which must not be destroyed between
-// calls of cfftf or cfftb.
+// calls of Cfftf or Cfftb.
// ifac Contains results which must not be destroyed between
-// calls of cfftf or cfftb.
-func cfftf(n int, r, work []float64, ifac []int) {
+// calls of Cfftf or Cfftb.
+func Cfftf(n int, r, work []float64, ifac []int) {
if len(r) < 2*n {
panic("fourier: short sequence")
}
@@ -174,8 +174,8 @@ func cfftf(n int, r, work []float64, ifac []int) {
cfft1(n, r[:2*n], work[:2*n], work[2*n:4*n], ifac[:15], -1)
}
-// cfftb computes the backward complex Discrete Fourier Transform
-// (the Fourier synthesis). Equivalently, cfftf computes the computes
+// Cfftb computes the backward complex Discrete Fourier Transform
+// (the Fourier synthesis). Equivalently, Cfftf computes the computes
// a complex periodic sequence from its Fourier coefficients. The
// transform is defined below at output parameter c.
//
@@ -189,13 +189,13 @@ func cfftf(n int, r, work []float64, ifac []int) {
// to be transformed.
//
// work A real work array which must be dimensioned at least 4*n.
-// in the program that calls cfftb. The work array must be
-// initialized by calling subroutine cffti(n,work,ifac) and a
+// in the program that calls Cfftb. The work array must be
+// initialized by calling subroutine Cffti(n,work,ifac) and a
// different work array must be used for each different
// value of n. This initialization does not have to be
// repeated so long as n remains unchanged thus subsequent
// transforms can be obtained faster than the first.
-// The same work array can be used by cfftb and cfftb.
+// The same work array can be used by Cfftf and Cfftb.
//
// ifac A work array containing the factors of n. ifac must have
// length of at least 15.
@@ -208,18 +208,18 @@ func cfftf(n int, r, work []float64, ifac []int) {
//
// where i=sqrt(-1)
//
-// This transform is unnormalized since a call of cfftf
-// followed by a call of cfftb will multiply the input
+// This transform is unnormalized since a call of Cfftf
+// followed by a call of Cfftb will multiply the input
// sequence by n.
//
// The n elements of c are represented in n pairs of real
// values in r where c[j] = r[j*2]+r[j*2+1]i.
//
// work Contains results which must not be destroyed between
-// calls of cfftf or cfftb.
+// calls of Cfftf or Cfftb.
// ifac Contains results which must not be destroyed between
-// calls of cfftf or cfftb.
-func cfftb(n int, r, work []float64, ifac []int) {
+// calls of Cfftf or Cfftb.
+func Cfftb(n int, r, work []float64, ifac []int) {
if len(r) < 2*n {
panic("fourier: short sequence")
}

RENAME
fourier/fftpack_test.go
fourier/internal/fftpack/fftpack_test.go
// functions by Paul N Swarztrauber, placed in the public
// domain at http://www.netlib.org/fftpack/.
-package fourier
+package fftpack
import (
"math"
@@ -22,7 +22,7 @@ func TestRfft(t *testing.T) {
for _, test := range rfftTests {
work := make([]float64, 2*test.n)
ifac := make([]int, 15)
- rffti(test.n, work, ifac)
+ Rffti(test.n, work, ifac)
var failed bool
if !floats.EqualApprox(work, test.wantiwork, 1e-6) {
failed = true
@@ -69,7 +69,7 @@ func TestRfft(t *testing.T) {
y[test.n-1] = sum1 - sum2
}
- rfftf(test.n, x, work, ifac)
+ Rfftf(test.n, x, work, ifac)
var rftf float64
for i := 0; i < test.n; i++ {
rftf = math.Max(rftf, math.Abs(x[i]-y[i]))
@@ -95,7 +95,7 @@ func TestRfft(t *testing.T) {
}
y[i] = 2 * sum
}
- rfftb(test.n, x, work, ifac)
+ Rfftb(test.n, x, work, ifac)
var rftb float64
for i := 0; i < test.n; i++ {
rftb = math.Max(rftb, math.Abs(x[i]-y[i]))
@@ -106,8 +106,8 @@ func TestRfft(t *testing.T) {
t.Errorf("unexpected rftb value for n=%d: got:%f want:0", test.n, rftb)
}
- rfftb(test.n, y, work, ifac)
- rfftf(test.n, y, work, ifac)
+ Rfftb(test.n, y, work, ifac)
+ Rfftf(test.n, y, work, ifac)
cf := 1.0 / fn
var rftfb float64
for i := 0; i < test.n; i++ {
@@ -284,7 +284,7 @@ func TestCfft(t *testing.T) {
for _, test := range cfftTests {
work := make([]float64, 4*test.n)
ifac := make([]int, 15)
- cffti(test.n, work, ifac)
+ Cffti(test.n, work, ifac)
var failed bool
if !floats.EqualApprox(work, test.wantiwork, 1e-6) {
failed = true
@@ -304,7 +304,7 @@ func TestCfft(t *testing.T) {
x, y1 := cmplxSeries(test.n)
cx := cmplxAsFloat(x)
- cfftf(test.n, cx, work, ifac)
+ Cfftf(test.n, cx, work, ifac)
x = floatAsCmplx(cx)
var cftf float64
@@ -321,7 +321,7 @@ func TestCfft(t *testing.T) {
y2 := updatedCmplxSeries(x)
cx = cmplxAsFloat(x)
- cfftb(test.n, cx, work, ifac)
+ Cfftb(test.n, cx, work, ifac)
x = floatAsCmplx(cx)
var cftb float64
@@ -335,8 +335,8 @@ func TestCfft(t *testing.T) {
}
cx = cmplxAsFloat(x)
- cfftf(test.n, cx, work, ifac)
- cfftb(test.n, cx, work, ifac)
+ Cfftf(test.n, cx, work, ifac)
+ Cfftb(test.n, cx, work, ifac)
x = floatAsCmplx(cx)
var cftfb float64

RENAME
fourier/rfft.go
fourier/internal/fftpack/rfft.go
// Paul N Swarztrauber, placed in the public domain at
// http://www.netlib.org/fftpack/.
-package fourier
+package fftpack
import "math"
-// rffti initializes the array work which is used in both rfftf
-// and rfftb. The prime factorization of n together with a
+// Rffti initializes the array work which is used in both Rfftf
+// and Rfftb. The prime factorization of n together with a
// tabulation of the trigonometric functions are computed and
// stored in work.
//
@@ -22,14 +22,14 @@ import "math"
// Output parameters:
//
// work A work array which must be dimensioned at least 2*n.
-// The same work array can be used for both rfftf and rfftb
+// The same work array can be used for both Rfftf and Rfftb
// as long as n remains unchanged. different work arrays
// are required for different values of n. The contents of
-// work must not be changed between calls of rfftf or rfftb.
+// work must not be changed between calls of Rfftf or Rfftb.
//
// ifac A work array containing the factors of n. ifac must have
// length of at least 15.
-func rffti(n int, work []float64, ifac []int) {
+func Rffti(n int, work []float64, ifac []int) {
if len(work) < 2*n {
panic("fourier: short work")
}
@@ -110,7 +110,7 @@ outer:
}
}
-// rfftf computes the Fourier coefficients of a real perodic sequence
+// Rfftf computes the Fourier coefficients of a real perodic sequence
// (Fourier analysis). The transform is defined below at output
// parameter r.
//
@@ -124,13 +124,13 @@ outer:
// to be transformed.
//
// work a work array which must be dimensioned at least 2*n.
-// in the program that calls rfftf. the work array must be
+// in the program that calls Rfftf. the work array must be
// initialized by calling subroutine rffti(n,work,ifac) and a
// different work array must be used for each different
// value of n. This initialization does not have to be
// repeated so long as n remains unchanged. Thus subsequent
// transforms can be obtained faster than the first.
-// The same work array can be used by rfftf and rfftb.
+// The same work array can be used by Rfftf and Rfftb.
//
// ifac A work array containing the factors of n. ifac must have
// length of at least 15.
@@ -150,15 +150,15 @@ outer:
// r[n-1] = the sum from i = 0 to i = n-1 of
// (-1)^i*r[i]
//
-// This transform is unnormalized since a call of rfftf
-// followed by a call of rfftb will multiply the input
+// This transform is unnormalized since a call of Rfftf
+// followed by a call of Rfftb will multiply the input
// sequence by n.
//
// work contains results which must not be destroyed between
-// calls of rfftf or rfftb.
+// calls of Rfftf or Rfftb.
// ifac contains results which must not be destroyed between
-// calls of rfftf or rfftb.
-func rfftf(n int, r, work []float64, ifac []int) {
+// calls of Rfftf or Rfftb.
+func Rfftf(n int, r, work []float64, ifac []int) {
if len(r) < n {
panic("fourier: short sequence")
}
@@ -624,7 +624,7 @@ func radfg(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) {
}
}
-// rfftb computes the real perodic sequence from its Fourier
+// Rfftb computes the real perodic sequence from its Fourier
// coefficients (Fourier synthesis). The transform is defined
// below at output parameter r.
//
@@ -638,13 +638,13 @@ func radfg(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) {
// to be transformed.
//
// work A work array which must be dimensioned at least 2*n.
-// in the program that calls rfftb. The work array must be
+// in the program that calls Rfftb. The work array must be
// initialized by calling subroutine rffti(n,work,ifac) and a
// different work array must be used for each different
// value of n. This initialization does not have to be
// repeated so long as n remains unchanged thus subsequent
// transforms can be obtained faster than the first.
-// The same work array can be used by rfftf and rfftb.
+// The same work array can be used by Rfftf and Rfftb.
//
// ifac A work array containing the factors of n. ifac must have
// length of at least 15.
@@ -662,15 +662,15 @@ func radfg(ido, ip, l1, idl1 int, cc, c1, c2, ch, ch2, wa []float64) {
// 2*r(2*k-1)*cos(k*i*2*pi/n)
// -2*r(2*k)*sin(k*i*2*pi/n)
//
-// This transform is unnormalized since a call of rfftf
-// followed by a call of rfftb will multiply the input
+// This transform is unnormalized since a call of Rfftf
+// followed by a call of Rfftb will multiply the input
// sequence by n.
//
// work Contains results which must not be destroyed between
-// calls of rfftf or rfftb.
+// calls of Rfftf or Rfftb.
// ifac Contains results which must not be destroyed between
-// calls of rfftf or rfftb.
-func rfftb(n int, r, work []float64, ifac []int) {
+// calls of Rfftf or Rfftb.
+func Rfftb(n int, r, work []float64, ifac []int) {
if len(r) < n {
panic("fourier: short sequence")
}

fourier/internal/fftpack: factor out FFTPACK port to internal package