fourier: add quarterwave functions #465

kortschak · 2018-04-14T22:53:59Z

Please take a look.

At this point the change of .*FFT -> .*Coefficients and .*IFFT -> .*Sequence methods names should make sense.

kortschak · 2018-04-15T01:11:57Z

OK, this is ready. PTAL

vladimir-ch · 2018-04-15T20:19:26Z

"github.com/gonum/floats" sneaked into fftpack_test.go (already in the initial PR).

vladimir-ch

Just a quick look at what's in the PR.

vladimir-ch · 2018-04-15T11:43:37Z

fourier/internal/fftpack/cosq.go

+	}
+}
+
+// Subroutine Cosqb computes the fast fourier transform of quarter


Delete "Subroutine" here and also at other places in this PR, shown by golint.

Done - and other doc fixes.

vladimir-ch · 2018-04-15T20:21:45Z

fourier/fourier_sincos.go

+	ifac [15]int
+}
+
+// NewQuarterWave returns a QuarterWave initialized for work on sequences of length n.


NewQuarterWaveFFT

vladimir-ch · 2018-04-15T20:57:42Z

fourier/fourier.go

-// the length of dst does not equal t.Len()/2+1, FFT will panic.
-func (t *FFT) FFT(dst []complex128, seq []float64) []complex128 {
+// the length of dst does not equal t.Len()/2+1, Coefficients will panic.
+func (t *FFT) Coefficients(dst []complex128, seq []float64) []complex128 {


I didn't make it in time for the initial PR, so I'll comment here. I think the type would be better named DFT because it computes the Discrete Fourier Transform (the problem being solved) while the Fast Fourier Transform is the (family of) algorithm that computes it (implementation detail that solves the problem). It would be also more consistent with the transforms introduced in this PR. Although I do know that DFT and FFT are often used interchangeably.

Unfortunately discoverability will likely trump correctness, though this is strictly still correct I suppose.

I think I'm with @kortschak , while DFT may be correct, everyone will be searching for FFT and clearly having that name would be useful. If there is a potential name clash I could be convinced otherwise, but it seems like bowing to commonality here .

We could also call this DFT, but somehow include FFT in the package comment to help searching .

kortschak · 2018-04-26T04:37:03Z

ping

btracey · 2018-04-23T15:24:59Z

fourier/fourier.go

@@ -40,15 +40,15 @@ func (t *FFT) Reset(n int) {
 	fftpack.Rffti(n, t.work, t.ifac[:])
 }

-// FFT computes the Fourier coefficients of the input sequence, seq,
+// Coefficients computes the Fourier coefficients of the input sequence, seq,


Can we add in a comment that's something like "(transforming the time series into the frequency spectrum)". I know it basically says that, but having something that overly clearly states "this is FFT and not IFFT" will help avoid confusion with our somewhat abnormal naming.

btracey · 2018-04-23T19:28:10Z

fourier/fourier.go

-// the length of dst does not equal t.Len()/2+1, FFT will panic.
-func (t *FFT) FFT(dst []complex128, seq []float64) []complex128 {
+// the length of dst does not equal t.Len()/2+1, Coefficients will panic.
+func (t *FFT) Coefficients(dst []complex128, seq []float64) []complex128 {


I think I'm with @kortschak , while DFT may be correct, everyone will be searching for FFT and clearly having that name would be useful. If there is a potential name clash I could be convinced otherwise, but it seems like bowing to commonality here .

We could also call this DFT, but somehow include FFT in the package comment to help searching .

btracey · 2018-04-26T14:44:27Z

fourier/fourier.go

@@ -74,15 +74,15 @@ func (t *FFT) FFT(dst []complex128, seq []float64) []complex128 {
 	return dst
 }

-// IFFT computes the real perodic sequence from the Fourier coefficients,
+// Sequence computes the real perodic sequence from the Fourier coefficients,


again, it would be nice to have something that makes it overly clear this is what people normally think of as ifft

btracey · 2018-04-26T15:05:53Z

fourier/fourier_sincos.go

@@ -0,0 +1,238 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.


why not just call this file sincos.go

It's called this for historical reasons. I wanted it to be in the fourier.go file initially, but that was in the air, so I wrote it into a new file - same with the tests. I actually think the name it has now is more accurate and have grown to think the separation is good.

If it should be separated, I would have sincos.go and quarter.go. Does this sound OK?

That sounds fine.

btracey · 2018-04-26T15:28:03Z

fourier/fourier_sincos.go

+	return dst
+}
+
+// QuarterWaveFFT implements Fast Fourier Transform for quarter wave data.


I think it would be nice to provide a reference here

Yes. The best I can find is this https://pdfs.semanticscholar.org/presentation/6e0d/2cbafb339cdf7ba90fdc36997d6b91dfffad.pdf.

btracey · 2018-04-28T15:17:21Z

fourier/fourier_sincos.go

+	fftpack.Cosqi(n, t.work, t.ifac[:])
+}
+
+// CosCoefficients computes the Fast Fourier Transform of quarter wave data for


Should this say Discrete? DST for instance uses discrete while this uses fast.

They are sort of the same in this domain. The docs reflect the documentation in fftpack.

btracey · 2018-04-28T15:22:21Z

fourier/fourier_test.go

+
+			want := make([]float64, n)
+			for i := range want {
+				want[i] = rand.Float64()


These should have their own RNG in the test. Also, should these maybe use NormFloat so there are both positive and negative numbers, as well as numbers > 1? Sorry, missed this in the last PR.

Same for the other tests

Done for the source, but changing the norm floats causes failure that I don't have time to look into right now.

btracey · 2018-04-28T15:25:03Z

fourier/fourier_test.go

+	t.Run("Reset DCT", func(t *testing.T) {
+		dct := NewDCT(1000)
+		for n := 2; n <= 2000; n++ {
+			dct.Reset(n)


Maybe add one test at the end where you reset to a smaller size?

Functionally this is what happens. In each case the dct is initially created with 1000, and then increased from the smallest value it can take.

Good point.

btracey · 2018-04-28T15:54:57Z

fourier/internal/fftpack/cosq.go

@@ -0,0 +1,220 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.


This file looks good.

btracey

Next chunk done. Haven't found anything that would cause it to break with negative numbers sadly

btracey · 2018-04-29T15:40:57Z

fourier/fourier_test.go

+	t.Run("Reset DCT", func(t *testing.T) {
+		dct := NewDCT(1000)
+		for n := 2; n <= 2000; n++ {
+			dct.Reset(n)


Good point.

btracey · 2018-04-29T15:58:48Z

fourier/internal/fftpack/cost.go

@@ -0,0 +1,143 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.


This file looks good.

btracey · 2018-04-29T16:15:55Z

fourier/internal/fftpack/sinq.go

@@ -0,0 +1,179 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.


This file looks good.

btracey · 2018-04-29T16:20:28Z

fourier/internal/fftpack/sint.go

+	default:
+		x[0] = 0
+		for k := 0; k < n/2; k++ {
+			kc := n - k


maybe kc := n - k - 1 like you did in the other file so the index just works below.

Also in cost.go.

btracey · 2018-04-29T16:20:59Z

fourier/internal/fftpack/sint.go

+			t1 := xh[k] - xh[kc-1]
+			t2 := was[k] * (xh[k] + xh[kc-1])
+			x[k+1] = t1 + t2
+			x[kc+1-1] = t2 - t1


Either way with the above we don't need kc+1-1 :)

btracey · 2018-04-29T16:22:44Z

fourier/internal/fftpack/sint.go

@@ -0,0 +1,146 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.


Looks good minus the comments below.

kortschak · 2018-04-29T23:41:52Z

Next chunk done. Haven't found anything that would cause it to break with negative numbers sadly

The reporting of failures is pretty sparse, so it may just be a small delta on the expectation. I just have not had time to look into it.

kortschak · 2018-04-29T23:56:19Z

PTAL

btracey · 2018-05-01T14:53:35Z

fourier/fourier_sincos.go

@@ -0,0 +1,238 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.


That sounds fine.

btracey · 2018-05-03T16:10:39Z

fourier/internal/fftpack/fftpack_test.go

@@ -616,3 +616,827 @@ var cfftTests = []struct {
 		wantiifac: []int{2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 	},
 }
+
+func TestSint(t *testing.T) {


Could you please add comments to the sections of the tests to make it more clear what's going on? (This goes for the other added tests here)

Is it OK if I just say that these tests reiterate the tests that are part of FFTPACK?

It turns out I do have that, in the header for the public domain attribution.

What I'd write is essentially this, but it's more than we write for lapack and blas, which are arguably more complex in their tests because of the structure of those suites.

// Broadly, the tests each compare the prepared work and // prime factor slices to golden values for a variety of // lengths known to excercise all parts of the forward and // backward transforms. Then for each of these lengths a // synthetic sequence with known frequency spectrum is // analysed with the forward transform and compared to the // spectrum and the spectrum is synthesised and compared to // the sequence. Finally, the the synthetic spectrum is // synthesised to the sequence and back again, being compared // to the original spectrum. In all cases the transforms // are normalised by the appropriate constant specified in // the documentation.

Thanks. Added some comments on what I would add.

I'll add those, but it does seem odd to me to do it here when the tests that exist in testblas and to a lesser extent testlapack largely don't have any commentary.

Do you have particular functions in mind? I did a quick survey and most of the float64 blas tests are only one section, and of the blas/lapack functions that have multiple sections typically those sections have some amount of explanation text.

This is always going to be a subjective assessment, but here is a list of tests in blas and lapack that look like they are commented to a similar level. There are others where there is commentary, but at least when I read it, I need to do a lot of work to understand what is being done - this is obviously an interaction between my level of linalg experience and the commentary. I've left out tests where there are tables or the tests are trivial.

blas: Zgbmv, Zhbmv, Ztbmv, Ztpsv, Ztrsv
lapack: Debgd2, Dgelqf, Dgeqp3, Dlags2, Dlange, Dlaln2, Dlanst, Dlansy, Dlanv2, Dlapll, Dlaqr23, Dlarfx, Dlartg, Dlasr, Dlatrd, Dorg2r, Dorgbr, Dorgtr, Dpotf2, Dtrtri2, Dtrtri.

It's worth noting that many of the tests explain what they do in the Errorf format string, though not necessarily how they get there - I think the fftpack tests do both (intentionally).

Thanks for the list. I'll go through and try to improve.

btracey

I think it would be nice to spread the testing comments out something along the lines of what I suggested. Otherwise this looks good.

btracey · 2018-05-04T15:06:48Z

fourier/internal/fftpack/fftpack_test.go

@@ -616,3 +616,827 @@ var cfftTests = []struct {
 		wantiifac: []int{2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
 	},
 }
+
+func TestSint(t *testing.T) {


Thanks. Added some comments on what I would add.

btracey · 2018-05-04T15:09:25Z

fourier/internal/fftpack/fftpack_test.go

+		work := make([]float64, 5*(test.n)/2)
+		ifac := make([]int, 15)
+
+		Sinti(test.n-1, work, ifac)


// Compute the transform of golden sequence and compare with known transform.

// Compute the work and factor slices and compare to known values.

btracey · 2018-05-04T15:18:15Z

fourier/internal/fftpack/fftpack_test.go

+			continue
+		}
+
+		fn := float64(test.n)


// Construct a sequence with known frequency spectrum and compare the computed spectrum

btracey · 2018-05-04T15:20:50Z

fourier/internal/fftpack/fftpack_test.go

+			t.Errorf("unexpected sintt value for n=%d: got:%f want:0", test.n, sintt)
+		}
+
+		Sint(test.n-1, x, work, ifac)


// Check that the transform is its own inverse up to a constant

kortschak · 2018-05-05T23:13:24Z

PTAL

kortschak · 2018-05-07T22:12:09Z

Thanks for the review

btracey · 2018-05-07T22:13:15Z

Thanks for all the translation work! Another big chunk of functionality no longer missing!

fourier/internal/fftpack: add quarterwave fourier functions

baed6d7

kortschak changed the title ~~Fourier/quarterwave~~ fourier: add quarterwave functions Apr 14, 2018

kortschak force-pushed the fourier/quarterwave branch from 72503a7 to da4de00 Compare April 15, 2018 01:09

kortschak added 2 commits April 15, 2018 10:40

fourier: add quarterwave API

0333e8f

fourier: fix up documentation and use English method names

d5211d1

kortschak force-pushed the fourier/quarterwave branch from da4de00 to d5211d1 Compare April 15, 2018 01:11

kortschak requested review from sbinet, vladimir-ch and btracey April 15, 2018 01:11

vladimir-ch reviewed Apr 15, 2018

View reviewed changes

kortschak added 2 commits April 16, 2018 07:19

fftpack: fix documentation

36c85fc

fourier: fix doc comment

72f48fd

btracey reviewed Apr 28, 2018

View reviewed changes

kortschak added 2 commits April 29, 2018 11:03

fourier: address pr comments

4bd1ea3

fourier: update error formats for new names

e78c829

btracey reviewed Apr 29, 2018

View reviewed changes

fftpack: address review comments

952740a

btracey reviewed May 3, 2018

View reviewed changes

fourier: split fourier_sincos and rename

103e716

btracey reviewed May 4, 2018

View reviewed changes

fftpack: add commentary

84faf85

vladimir-ch mentioned this pull request May 7, 2018

blas/testblas: document implementation of some L2 complex tests #479

Merged

btracey approved these changes May 7, 2018

View reviewed changes

kortschak merged commit ea4bf2f into master May 7, 2018

kortschak deleted the fourier/quarterwave branch May 7, 2018 22:13

vladimir-ch mentioned this pull request May 16, 2018

testlapack: add implementation comments to some tests #494

Closed

21 tasks

		@@ -0,0 +1,238 @@
		// Copyright ©2018 The Gonum Authors. All rights reserved.

		@@ -0,0 +1,220 @@
		// Copyright ©2018 The Gonum Authors. All rights reserved.

		@@ -0,0 +1,143 @@
		// Copyright ©2018 The Gonum Authors. All rights reserved.

		@@ -0,0 +1,179 @@
		// Copyright ©2018 The Gonum Authors. All rights reserved.

		@@ -0,0 +1,146 @@
		// Copyright ©2018 The Gonum Authors. All rights reserved.

fourier: add quarterwave functions #465

fourier: add quarterwave functions #465

Conversation

kortschak commented Apr 14, 2018 • edited

kortschak commented Apr 15, 2018

vladimir-ch commented Apr 15, 2018

vladimir-ch left a 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak Apr 15, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak commented Apr 26, 2018

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

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

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak Apr 29, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

btracey left a 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kortschak Apr 29, 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 Apr 29, 2018

kortschak commented Apr 29, 2018

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

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

btracey left a 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

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 May 5, 2018

kortschak commented May 7, 2018

btracey commented May 7, 2018

kortschak commented Apr 14, 2018 •

edited

kortschak Apr 15, 2018 •

edited

kortschak Apr 29, 2018 •

edited

kortschak Apr 29, 2018 •

edited