fourier: new package for fourier analysis and related functions #412

kortschak · 2018-02-04T22:40:32Z

Please take a look.

So far I have only added the backing code. I want to add CFFT[IFB] and then I will dismantle some of the array.go support (remove 1-based indexing and the linear support types). It makes sense to add the API in this PR, so suggestions welcome here.

kortschak · 2018-02-04T23:05:50Z

From an API perspective, I think that given there are work arrays kept, a type holding those that can be reused would be the most sensible approach:

type FFT struct

// NewFFT returns an FFT initialized for work on sequences of length n.
func NewFFT(n int) *FFT

// Len returns the length of the acceptable input.
func (t *FFT) Len() int

// Reset reinitializes the FFT for work on sequences of length n.
func (t *FFT) Reset(n int)

// FFT computes the Fourier coefficients of the input sequence, s,
// placing the result in dst and returning it.
// If the length of s is not t.Len(), FFT will panic.
// If dst is nil, a new slice is allocated and returned. If dst is not nil and
// the length of dst does not equal the length of s, FFT will panic.
func (t *FFT) FFT(dst, s []float64) []float64

// IFFT computes the real perodic sequence from the Fourier coefficients,
// c, placing the result in dst and returning it.
// If the length of s is not t.Len(), IFFT will panic.
// If dst is nil, a new slice is allocated and returned. If dst is not nil and
// the length of dst does not equal the length of c, IFFT will panic.
func (t *FFT) IFFT(dst, c []float64) []float64

I'm not wildly happy about the method names FFT and IFFT but they seem to be what people generally use (I'd prefer Coefficients and Sequence).

kortschak · 2018-02-04T23:53:27Z

fourier/rfft.go

+			ai2 = dc2*ai2 + ds2*ar2
+			ar2 = ar2h
+			for ik := 1; ik <= idl1; ik++ {
+				ch2m.set(ik, l, ch2m.at(ik, l)+ar2*c2m.at(ik, j))


These two lines suggest that a twoArray.add(i, j int, v float64) method would be a worthwhile optimisation.

These are done.

kortschak · 2018-02-04T23:53:39Z

fourier/rfft.go

+	}
+	for j := 2; j <= ipph; j++ {
+		for ik := 1; ik <= idl1; ik++ {
+			ch2m.set(ik, 1, ch2m.at(ik, 1)+c2m.at(ik, j))


kortschak · 2018-02-04T23:55:29Z

fourier/rfft.go

+			ai2 = dc2*ai2 + ds2*ar2
+			ar2 = ar2h
+			for ik := 1; ik <= idl1; ik++ {
+				c2m.set(ik, l, c2m.at(ik, l)+ar2*ch2m.at(ik, j))


kortschak · 2018-02-04T23:55:43Z

fourier/rfft.go

+
+	for j := 2; j <= ipph; j++ {
+		for ik := 1; ik <= idl1; ik++ {
+			ch2m.set(ik, 1, ch2m.at(ik, 1)+ch2m.at(ik, j))


kortschak · 2018-02-05T02:18:11Z

~~Damn. Not quite correct yet.~~

btracey · 2018-02-07T03:53:02Z

I'll get to this when I can, I'm on travel at the moment.

kortschak · 2018-02-07T05:34:53Z

No worries. Despite my best efforts, I managed to sneak some off-by-one errors into this - when you read the subtle one-based indexing removal, also read the fix off-by-one bugs commit. The sign bug goes with the use of the sign parameter commit. Rebased away.

kortschak · 2018-02-14T22:12:42Z

Apologies for the unhelpful ordering - I revised the author name, and github insists on presenting PR changes in chronological order in a way that made that revision show differences that were not otherwise there.

btracey · 2018-02-27T06:40:16Z

I understand the desire for having the types to reuse the work arrays. From what I understand, this is how FFTW works as well. Do you think it's desirable to also have some high level functions for ease of use? I haven't looked at all of the intended functions, but I wonder if we could get away with Float64 and Complex128 as types that hold the temporary memory, and the could have FFT be a function that allocates the temporary memory.

Just a thought. I'll take a more detailed look soon hopefully.

kortschak · 2018-02-27T06:52:49Z

They are an easy addition, so I think we should leave that for when someone asks for it.

kortschak · 2018-03-04T21:03:36Z

I have put together the implementation and API for DCT, DST and QW-FFT (tests will come later), but I won't include them in this PR. However, after seeing how the API looks I'm going to push a little harder towards English for the method names. I have two proposals (in preference order):

FFT -> Coefficients and IFFT -> Sequence.
FFT -> Analysis and IFFT -> Synthesis (or the verbs).

This gives us nicer API for the new routines mentioned above, e.g. for QW-FFT, SinCoefficients and CosCoefficients (the methods are on the same type, QuarterWaveFFT).

btracey · 2018-03-09T05:21:09Z

fourier/rfft.go

+	rffti1(n, work[n:2*n], ifac[:15])
+}
+
+func rffti1(n int, wa []float64, ifac []int) {


Checked this one and looks good.

btracey · 2018-03-09T05:31:43Z

fourier/rfft.go

+	rfftf1(n, r[:n], work[:n], work[n:2*n], ifac[:15])
+}
+
+func rfftf1(n int, c, ch, wa []float64, ifac []int) {


This looks good too

vladimir-ch · 2018-04-12T12:30:51Z

fourier/fourier_example_test.go

+	// freq=0.375 cycles/period, magnitude=3, phase=3.142
+}
+
+func Example_fFT2() {


Lower-case f in fFT2

That's to give it package scope. See the example names section https://blog.golang.org/examples

Oh, there's always something to learn. Thanks.

kortschak · 2018-04-12T21:34:58Z

Thanks. There are a few of things to do; ~~I forgot the comment that the cfft backend code is refactored to allow direction being set by sign~~(done), the change s/since/;/g (@btracey are you happy with this?) and ~~a lot of delicate history revision including creating internal/fftpack (I'm not confident about whether this is possible in this PR)~~(done).

On the subject of the internal/fftpack package: I know we've done this kind of thing in mathext, but I'm really not convinced it's a good idea. The package boundary is to define separation of concerns and modularity. I don't see that division here except as a consequence of history. Though maybe it does help in debugging with other implementations. Comments?

The implementation of FFT and IFFT here are a direct translation of the RFFTI/RFFTF/RFFTB and CFFTI/CFFTF/CFFTB functions from FFTPACK. The Fortranisms will be reduced, but not eliminated (column major will remain). API will follow.

kortschak · 2018-04-13T04:04:46Z

I've done as much history condensation as I think is achievable while retaining the capacity to identify bugs later.

btracey · 2018-04-13T05:07:28Z

Semicolon sounds great.

I think it's a separation of concerns. The code base isn't really "native go", just translated fortran, while the main package is novel code. The FFTPack code has additional "license" documentation that goes with it that the other code does not have. I think it's also good to separate for other implementations, in particular it would be nice to clean up a lot of the FFTPack code (i.e. using complex128 where appropriate) at which point those alternate functions could live in the main package, preferably with comment styles that match gonum rather than matching netlib. In short, the FFTPACK code feels out of place, and putting it in internal is a clear signal that it's somehow different.

kortschak · 2018-04-13T05:23:58Z

OK. I'll try to rebase it into a new package and at the same time add the semicolons. I'll ping when it's done since that push won't alert you.

kortschak · 2018-04-13T06:50:02Z

Actually, the history revision is non-trivial to the extent that I think I will just add another commit moving the fftpack code out to the internal package. There is too much breakage opportunity doing it the way I initially wanted.

kortschak · 2018-04-13T12:49:17Z

PTAL

kortschak · 2018-04-14T22:16:34Z

Thanks for the review.

kortschak requested review from sbinet, vladimir-ch and btracey February 4, 2018 22:40

kortschak commented Feb 4, 2018

View reviewed changes

kortschak force-pushed the fourier branch from 27ff5cd to 756394f Compare February 5, 2018 02:17

kortschak force-pushed the fourier branch 8 times, most recently from dd0e55e to 2c4ea68 Compare February 6, 2018 05:25

kortschak force-pushed the fourier branch 7 times, most recently from 2a1a53f to 31ce9e8 Compare February 14, 2018 22:10

btracey reviewed Mar 9, 2018

View reviewed changes

vladimir-ch reviewed Apr 12, 2018

View reviewed changes

kortschak dismissed btracey’s stale review via 391e516 April 13, 2018 03:18

kortschak added 13 commits April 13, 2018 13:17

fourier: remove oneIntArray

48bc73d

fourier: remove 1-based slicing

07e7c0b

fourier: change direction int to bool

40e355b

fourier: remove trivial one-based indexing

61d4279

fourier: remove subtle one-based indexing

268053a

fourier: remove refactoring support types and make arrays zero-based

7a594dd

fourier: remove fortran scars

d2347a8

fourier: make bound checks optional

9e20d47

fourier: add direct update optimisation

2943e17

fourier: use sign parameter instead of code duplication

2eac324

fourier: tighten bounds

c80b8a4

fourier: improve documentation wording

d627f71

kortschak force-pushed the fourier branch from 391e516 to cfaec5f Compare April 13, 2018 03:59

kortschak force-pushed the fourier branch from cfaec5f to 24397c5 Compare April 13, 2018 05:28

fourier: add API and code examples

4dec74c

kortschak force-pushed the fourier branch from 24397c5 to 4dec74c Compare April 13, 2018 06:31

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

18eefc9

btracey approved these changes Apr 14, 2018

View reviewed changes

kortschak merged commit 7b0fe87 into master Apr 14, 2018

kortschak deleted the fourier branch April 14, 2018 22:15

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fourier: new package for fourier analysis and related functions #412

fourier: new package for fourier analysis and related functions #412

kortschak commented Feb 4, 2018

kortschak commented Feb 4, 2018 •

edited

kortschak Feb 4, 2018 •

edited

kortschak Feb 6, 2018 •

edited

kortschak Feb 4, 2018

kortschak Feb 4, 2018

kortschak Feb 4, 2018

kortschak commented Feb 5, 2018 •

edited

btracey commented Feb 7, 2018

kortschak commented Feb 7, 2018 •

edited

kortschak commented Feb 14, 2018 •

edited

btracey commented Feb 27, 2018

kortschak commented Feb 27, 2018

kortschak commented Mar 4, 2018 •

edited

btracey Mar 9, 2018

btracey Mar 9, 2018

vladimir-ch Apr 12, 2018

btracey Apr 12, 2018

vladimir-ch Apr 12, 2018

kortschak commented Apr 12, 2018 •

edited

kortschak commented Apr 13, 2018

btracey commented Apr 13, 2018

kortschak commented Apr 13, 2018

kortschak commented Apr 13, 2018

kortschak commented Apr 13, 2018

kortschak commented Apr 14, 2018

fourier: new package for fourier analysis and related functions #412

fourier: new package for fourier analysis and related functions #412

Conversation

kortschak commented Feb 4, 2018

kortschak commented Feb 4, 2018 • edited

kortschak Feb 4, 2018 • edited

Choose a reason for hiding this comment

kortschak Feb 6, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak commented Feb 5, 2018 • edited

btracey commented Feb 7, 2018

kortschak commented Feb 7, 2018 • edited

kortschak commented Feb 14, 2018 • edited

btracey commented Feb 27, 2018

kortschak commented Feb 27, 2018

kortschak commented Mar 4, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak commented Apr 12, 2018 • edited

kortschak commented Apr 13, 2018

btracey commented Apr 13, 2018

kortschak commented Apr 13, 2018

kortschak commented Apr 13, 2018

kortschak commented Apr 13, 2018

kortschak commented Apr 14, 2018

kortschak commented Feb 4, 2018 •

edited

kortschak Feb 4, 2018 •

edited

kortschak Feb 6, 2018 •

edited

kortschak commented Feb 5, 2018 •

edited

kortschak commented Feb 7, 2018 •

edited

kortschak commented Feb 14, 2018 •

edited

kortschak commented Mar 4, 2018 •

edited

kortschak commented Apr 12, 2018 •

edited