proposal: signal processing package

[Kevin Zimmerman]

I've been advised to begin a discussion here on gonum-dev regarding my desire to start designing/implementing a signal processing package in gonum.

I've personally had a need for the convolution operation in Go and somewhat expected to find it in gonum but didn't. Does anyone have any comments about the need for such a package and what the design of it should be.

I've never contributed before so this is all new to me.

Thanks in advance,

kczimm

[Kevin Zimmerman]

For instance, the scipy.signal package has a convolve function and so does MATLAB. The function calls are very similar between both programming languages. I imagined I'd find a function in gonum like,

func Convolve(a, b []float64) []float64

but didn't. 

[Sebastien Binet]

Kevin,

FWIW, I think having a "signal" package (with, most probably, a "fft" sub package) makes sense.

(indeed, my next-next-next blog post should be on reproducing the results/paper from this year's Physics Nobel prize, based on this repo: https://github.com/cranmer/ligo-binder.

the notebook can be run there: https://beta.mybinder.org/repo/cranmer/ligo-binder)

and I think many people (everybody?) on this list have already agreed that (ultimately) a fft package should be part of the gonum tool box.

I must say that I have been just a consumer (at best) of FFT packages (in python) and thus I don't feel confident enough to shepherd this endeavour.

Perhaps we could use the gonum/exp repo to host our API experiments?

In any case, I believe we should:

- start small (convolution sounds like a more manageable piece than FFT)

- clearly say what we want to achieve in our first attempt

- list the needed functionalities for our first attempt

and try to see how/whether it meshes with the core interfaces of Gonum (I surmise: gonum/mat + gonum/stat)

what scipy.signal provides is a rather huge laundry list of functionalities:

https://docs.scipy.org/doc/scipy/reference/signal.html

without getting too much ahead of ourselves, I think the first order of business is: what type(s) should we use.

it's not uncommon in signal processing to involve complex data.

do we arange our signal processing package(s) to work on f64, c128, f64+c128, ... ?

cheers,

-s

[Kevin Zimmerman]

Unfortunately, I'm in a similar situation as you with only be a consumer of FFT packages. However, I wouldn't mind making an attempt towards an fft package (or sub-package as you suggested) in the gonum/exp repo with assistance from others.

I agree in general about starting small with a limited scope at first. Scipy's signal package is very large, as you mentioned, and structured sub-optimally IMHO. The most basic set of functions that I think a gonum/signal (or gonum/exp/signal) package should have are the following:

- convolve

- deconvolve (maybe)

- correlate

- filter

- fft

- ifft

And in the spirit of keeping the scope small, start with only the one-dimensional case.

With regards to your described first order of business, I think we are obligated, due to the nature of signal processing in general, to arrange the signal processing package(s) to work on both real and complex data.

-kz

[Scott Cotton]

Probably too little too late, but in case not and for future viewers, Here is one convolution package in go.

Scott  

[Brendan Tracey]

Also note that gonum now does have an FFT package.

https://godoc.org/gonum.org/v1/gonum/fourier

[Dan Kortschak]

Bringing this back to gonum-dev.

https://groups.google.com/d/topic/gonum-dev/iGeEn43k9H4/discussion

Comments in-line.

On Thu, 2018-08-16 at 16:37 -0700, William Cotton wrote:
> I do not have any feel for where the discussion is most appropriate,
> just
> following guidance of y'all.

Yeah, no worries.

> Re performance claims, micro bench in zikichombo.org/dsp/fft is a
> start
> with go test -bench.  If you look at the constants, for real_test.go
> and
> t_test.go you will see they are for the same problem size.  A real
> interface which just puts float64s in the complex128s would just add
> to the
> cost of the t_test.go.
> 
> Although not exactly measuring what happens or would happen in gonum
> fft,
> it is at least a point of reference.

Because of the differences in implementation the performance
characteristics of the approaches need to be assessed. The cost of a
copy from float64 to the real part of complex128 is likely to be very
minor when compared to the cost of the transform operation, but this
assertion is why there needs to be a comparison.

Can you explain to me what the difference between the half complex and
the real FFT is in your package?

> Unfortunately, I am not familiar with FFTPACK implementations, so we
> have
> some disjointness in our points of reference to consider in the
> discussion.

Yes, that's why we're having this discussion.

>   That said, it sounds like you took my comment about relations
> between HC
> MDCT and FFT roughly as intended (intensions being based on perusal
> of FFTW
> and relevant MDCT/DCTIV relations)

If you are prepared to do some initial work on what that would entail,
I'm happy to pick up the work when I have time. Otherwise feel free to
send a PR adding MDCT to either or both of fourier and internal/fftpack
(where I think we would be better off starting for efficiency of
operation's sake).

> Apologies, I did not understand anything regarding your point about
> interfaces and docs and examples. Perhaps we could back up and try

I suspect you are using the term interface in the more general sense
(the I in API), rather than the Go type sense (type I interface {...),
since there does not appear to be a commonality of methods between the
difference concrete types. Maybe I'm misunderstanding here.

> to
> re-sync the topic.  I found the spectra functionality in
> zc.org/dsp/fft
> useful but of course it is also too audio specific for gonum to
> prioritise.  Just thought maybe there were some aspects you might
> like,
> such as not using fftshift, "FoldReal".  I guess gonum might have
> fancier
> ways of interpolating the peaks, but what zc has works well for
> audio.

I guess what I'd like to see is some examples of use. The method names
in zikichombo's dsp packages are pretty terse and documentation is
sparse, so getting a feel of how to use it is non-trivial.

> On a larger note, I find the "standard" dsp libraries like Matlab etc
> often
> have implicitly broken APIs, such as fftshift.  But they have become
> so
> standard that people start to think of them as "mathematical" ideas
> just
> because they are an artefact of matlab.  Maybe I'm just too much an
> oddball
> in this respect, and that's fine, but I thought it worth voicing this
> anyway.  Feel free to ignore.

Opinions are not ignored, but may be disagreed with. The ShiftIdx
operation could be renamed as CenteredIdx (though that would likely
reduce the readability of the API in the case of the gonum package),
but it's as mathematically inappropriate as PCA. The justification of
centring is that it's multiply useful (the centring operation has uses
in image processing as well as as a basis for folding, which is merely
an additional operation on top of that). This is a case of less is
more; we have very few methods, but they can be composed in such a way
that multiple uses are possible. There is certainly a case for thinking
about adding a dsputils (or better named) package that does some of the
operation that zikichombo has based on the fourier package's base
types. Your guidance in that direction would be very helpful.

Dan

> 
> On 17 August 2018 at 00:43, Dan Kortschak <notifi...@github.com>
> wrote:
> 
> > 
> > I think this discussion is still at the level of what should be
> > happening
> > at gonum-dev. Until there's a reasonable set of concrete proposals
> > can we
> > keep it there. There is a lot here to take in.
> > 
> > Some of the things here need to be addressed though to avoid being
> > lost.
> > 
> > What I will say though is that there is a general approach in Gonum
> > to not
> > overly specialise implementations to any particular field, doing so
> > reduces
> > flexibility and usability in others. We follow the *Do less. Enable
> > more*
> > <https://blog.golang.org/open-source> axiom that is ingrained in Go
> > practice. A frequency type appears to do to much to me.
> > 
> > I would like to keep the reduced input length restrictions we have
> > (FFTPACK does not require even length input and this is the basis
> > for our
> > code).
> > 
> > I'm not sure of the interface argument, there are commonly named
> > methods
> > in the dsp packages, but they cannot be an interface as they take
> > different
> > types. So I'm not entirely sure how that would work. Maybe some
> > examples or
> > extended documentation in dsp would help here.
> > 
> > To that end, and the performance claims, it would be good to see
> > some
> > benchmark comparisons (both to see difference in use and perf). I
> > doubt the
> > inplaceness of the operations is a significant contributor to the
> > cost of
> > them, but that is something that can be addressed by doing some
> > rework on
> > the internal package.
> > 
> > What I take from the comment above is a desire for MDCT (which we
> > can get
> > from the FFT code that we have) and the HalfComplex (which we can
> > probably
> > derive reasonably easily from the complex FFT code that we have).
> > 
> > —
> > You are receiving this because you authored the thread.
> > Reply to this email directly, view it on GitHub
> > <https://github.com/gonum/gonum/issues/569#issuecomment-413706233>,
> > or mute
> > the thread
> > <https://github.com/notifications/unsubscribe-auth/ASR014Pe_PgeGMm3
> > UwlNzoDJwwoD7Dvhks5uRfWCgaJpZM4V5KWA>
> > .
> > 
> 
> 
> -- 
> Scott Cotton
> President, IRI France SAS
> http://www.iri-labs.com
> 
> 

[Scott Cotton]

Thanks, responses inline.

Although a comparison of this cost would be informative, I do not think it 

would be authoritative as it relates to memory cache, which in turn can be mostly

a function of the context in which it is called.   You could end up with benchmarks

that say the cost is relatively high or low and calling contexts which reverse that.

For this reason, I believe when there is a choice of a copy or move of memory

or not, some effort or consideration should be directed toward eliminating the operation 

regardless.  

Can you explain to me what the difference between the half complex and
the real FFT is in your package?

the relationship is as follows.  HalfComplex is a data layout which places 

N complex numbers in N floats by exploiting the symmetry properties of 

the spectrum of a real-only fft input. It is thus a memory layout which is 

appropriate for representing a spectrum.

Rather than hiding it in the zc API, I decided to make it a Go type which maintains 

the memory layout but provides easy access to the N complex via methods.

type HalfComplex []float64

in zc/fft.Real the method signature implementing the transform for float64 input

is:

   Do(d []float64) HalfComplex

and the inverse transform is

   Inv(HalfComplex) []float64

As HalfComplex uses the same memory as the slice, the result is that memory-wise,

the transform is in-place.  So, to summarise, zc.org/dst/fft.Real is an "object" which

does FFT transforms for a fixed size and inverse transforms on real-only data.  .Do is the forward

transform, .Inv is the inverse transform.  fft.HalfComplex is a type for representing the 

frequency domain/spectrum of these operations.

As for the other names, especially in convolution, yes they are quite terse, but 

the package names also inform the reader or coder.  Once you factor those in,

in cases where the package name is actually expressed in the code, it is not so terse.

In convol/ packages.  the naming convention is

convol.T  does linear convolution of two arguments A,B each of fixed size

convol.K  does the same, but uses a fixed convolution "kernel" for

A and it can be called repeatedly for different values of B, so long as

they are of a fixed size.

convol.Ola does blockwise overlap-add convolution for a given kernel.

It repeatedly calls convol.K for a fixed kernel on a sequence of arguments

B0 B1 B2 ..., which, when concatenated into B*, are equal to what convol.K(B*) 

would be when A is the "kernel".  It is used for when B* is large or a stream.  

 

 

> Unfortunately, I am not familiar with FFTPACK implementations, so we
> have
> some disjointness in our points of reference to consider in the
> discussion.

Yes, that's why we're having this discussion.

>   That said, it sounds like you took my comment about relations
> between HC
> MDCT and FFT roughly as intended (intensions being based on perusal
> of FFTW
> and relevant MDCT/DCTIV relations)

If you are prepared to do some initial work on what that would entail,
I'm happy to pick up the work when I have time. Otherwise feel free to
send a PR adding MDCT to either or both of fourier and internal/fftpack
(where I think we would be better off starting for efficiency of
operation's sake).

We'll see when I can get to it.

 

> Apologies, I did not understand anything regarding your point about
> interfaces and docs and examples. Perhaps we could back up and try

I suspect you are using the term interface in the more general sense
(the I in API), rather than the Go type sense (type I interface {...),
since there does not appear to be a commonality of methods between the
difference concrete types. Maybe I'm misunderstanding here.

> to
> re-sync the topic.  I found the spectra functionality in
> zc.org/dsp/fft
> useful but of course it is also too audio specific for gonum to
> prioritise.  Just thought maybe there were some aspects you might
> like,
> such as not using fftshift, "FoldReal".  I guess gonum might have
> fancier
> ways of interpolating the peaks, but what zc has works well for
> audio.

I guess what I'd like to see is some examples of use. The method names
in zikichombo's dsp packages are pretty terse and documentation is
sparse, so getting a feel of how to use it is non-trivial.

See above.  It will have to do for the moment.  I focus on documentation/examples

only after the underlying implementation is settled in terms of organisation.  When

you sum up development time for a given goal, for me this approach is (much) more efficient.

One open issue for zc/dsp is that HalfComplex forced a change in organisation of 

implementation.  As it is "alpha", I feel it is ok to publish what is there even though

this change in organisation makes the package less approachable to readily exploit.

> On a larger note, I find the "standard" dsp libraries like Matlab etc
> often
> have implicitly broken APIs, such as fftshift.  But they have become
> so
> standard that people start to think of them as "mathematical" ideas
> just
> because they are an artefact of matlab.  Maybe I'm just too much an
> oddball
> in this respect, and that's fine, but I thought it worth voicing this
> anyway.  Feel free to ignore.

Opinions are not ignored, but may be disagreed with. The ShiftIdx
operation could be renamed as CenteredIdx (though that would likely
reduce the readability of the API in the case of the gonum package),
but it's as mathematically inappropriate as PCA. The justification of
centring is that it's multiply useful (the centring operation has uses
in image processing as well as as a basis for folding, which is merely
an additional operation on top of that). This is a case of less is
more; we have very few methods, but they can be composed in such a way
that multiple uses are possible.

This can become difficult when the space of use cases have conflicting needs

in terms of composability of operations they rely on.  That's not so say I can propose

anything better in terms of general approach.  But it is not a silver bullet philosophy.

 

 

There is certainly a case for thinking
about adding a dsputils (or better named) package that does some of the
operation that zikichombo has based on the fourier package's base
types. Your guidance in that direction would be very helpful.

I'm not sure I'm the best qualified person for that.  DSP is a deep and broad field

the FFT/convolution part of it is a major cornerstone, and even in that space

it takes effort as evidenced here to synchronise on points of reference and have

enough of a feel for things to organise implementations.  My knowledge of filtering

is not up to par with a real filtering expert, and my knowledge of other important areas

of DSP is debutante.  I would make many mistakes in how to organise a dsputil package 

as a result.  

I think it would be best to find someone with broad DSP background and knowledge of Go

to help guide an effort to provide something worthy of being a general resource for dsp in Go.

 ZC/dsp currently only has resources for the parts of DSP most relevant to audio, and

in some of those parts we would definitely benefit from more experience and know-how. 

Scott

[Scott Cotton]

As a follow up, I'm happy to help coordinate in whatever way y'all see most fit,

As far as placing a banner in gonum for generic dsp related package my feeling is 

still that such an effort would benefit from the engagement of someone or some people

with broad DSP background, and that lacking that or lacking a sustained effort to find such

folks, a less conspicuous package name may be more appropriate.

Do you coordinate package naming and such here or in another venue?

Best,

Scott

   

--
You received this message because you are subscribed to the Google Groups "gonum-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to gonum-dev+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Dan Kortschak]

Here would be the appropriate place.

> > email to gonum-dev+...@googlegroups.com.

[Scott Cotton]

On my end I've found some (so far nominal) interest of a dsp professor via ZikiChombo.

Perhaps coordinating zc/dsp with a go-num dsp package would be best to put in an experimental package like the apache-arrow discussion?  

Also on my end, I am new to go-num numeric code, and I am a bit uncertain about how it would be easiest to coordinate zc/dsp with

go-num.  To the extent that zc/dsp code is or will have overlapping functionality with a gonum dsp package, the more inter-op between the two efforts the less overall effort in the long run.  But it is also harder for me to juggle the two concerns than one.  

 

If a place is made in the repo for the go-num dsp package, it would be great to work with some folk from go-num to help discern what parts of zc/dsp are or will be appropriate and to help me best make use of the existing go-num types and packages.

Scott