Linear convolution

[Oliver Lylloff]

Hello all,

Is anyone aware of a linear convolution implementation? 

The Base.conv and Base.conv2 functions are implemented with fft which makes them circular convolution functions (as far as I know).

I'm looking for something alike Matlabs conv2 or SciPys signal.convolve2d.

Should be straightforward to implement though.

Best,

Oliver

[Tim Holy]

Images.jl's imfilter might be what you want.
--Tim

[Oliver Lylloff]

Thanks Tim.

Can't believe I missed that - been working with Images.jl all day. Nice job by the way, very useful.

Best,

Oliver

[Andreas Noack Jensen]

Both conv and conv2 are linear convolutions but the implementations use the fft. Maybe the documentation could be more clear on that.

--
Med venlig hilsen

Andreas Noack Jensen

[Oliver Lylloff]

Well ok,

Maybe I misunderstood the whole thing then but since fft assumes periodic input then I don't see how it can be a linear convolution. I guess Base.conv2 probably uses zero-padding to reduce wrap-around but in theory it would still be a circular convolution. I'll read up on it :)

Best, 

Oliver

[Toivo Henningsson]

Yes, with sufficient padding, you can compute a linear convolution (of finite length vectors) exactly using a circular convolution. The FFT might introduce a little noise in the result, but that is all.

[Tim Holy]

The main reason to choose one or the other is merely kernel size; for small
kernels, a direct FIR convolution will be many times faster.

--Tim

> >> Andreas Noack Jensen

[Tobias Knopp]

Tim,

a little bit offtopic question but might it make sense to break of the algorithmic parts of Images.jl and put it into some "signal processing" package?
I know that the imagemagick dependency is a soft one but all the filtering stuff is IMHO so basic that it belongs to base, or rather into one "signal" package that could be one of the "default packages" that we hopefully get (see https://github.com/JuliaLang/julia/issues/1906)

Tobi

[Tim Holy]

I'm fine with that. Do you want to start it?

--Tim

[Tobias Knopp]

I don't want to give a definate yes to it but will think a little bit how such a package could look like.
My Cartesian macro foo is currently completely absent so that I would have to learn this first.

Do you know of others efforts/packages that go into that direction?

[Miguel Bazdresch]

There is DSP.jl: http://dspjl.readthedocs.org/en/latest/index.html

[João Felipe Santos]

We have some standard DSP stuff in https://github.com/JuliaDSP/DSP.jl. Our idea is to put everything DSP-related that is not application-specific in there, and then make other application-specific packages depend on it.

--
João Felipe Santos

On Tue, Mar 4, 2014 at 10:22 AM, Tobias Knopp <tobias...@googlemail.com> wrote:

[Miguel Bazdresch]

I was worried about getting "a little noise in the result", so I ran a quick test in Matlab and Julia, and got almost exactly the same error. This is the Matlab code:

Ts=0.01;
t=-10:Ts:10;
s=sinc(t);
sc=Ts*conv(s,s);
sc=sc(1000:3000);
sum((sc-s).*(sc-s))

ans =

    0.3695

So, at least for accuracy, julia's conv implementation seems to be no worse than Matlab's.

[Tobias Knopp]

Thanks for the hint. This seems to be a good start although the available functions already seem to be quite specialized.

The ndimage module from scipy (http://docs.scipy.org/doc/scipy/reference/ndimage.html#module-scipy.ndimage.filters) goes into a direction what in my mind could cover a Signal.jl package.

[Tobias Knopp]

This depends a lot on the input sizes. For full length convolutions, the fft approach should be more accurate because of less additions. For very short kernels this does not hold anymore. But in practice these kinds of errors are mostly negligable.

[João Felipe Santos]

That's funny, because in my opinion the functions in scipy.ndimage.filters are also specialized, as most of the time they seem to be used in image processing (but also in time-series processing) :)

In DSP.jl we mainly have functions for one-dimensional signal processing, though, and it would be nice to add multidimensional support to it. Filters such as max, min, median, and gaussian would also be interesting additions.

--
João Felipe Santos

[Miguel Bazdresch]

It's interesting that calculating the FFT and then the IFFT is more accurate; I wasn't aware of that. Thanks!

-- mb

[Don Gateley]

On Tuesday, March 4, 2014 3:38:15 AM UTC-8, Oliver Lylloff wrote:

Hello all,

Is anyone aware of a linear convolution implementation?

Not clear what you want but if it is to convolve a kernel with a longer stream you are looking for an implementation of overlap-add or overlap-save fast convolution which uses fft internally to accomplish fast.  I am not savvy enough with the Julia ecosystem to know if that exists or what it would be called but is not at all hard to implement.  Wikipedia covers it pretty well and has pseudo-code for it.  http://en.wikipedia.org/wiki/Overlap_add

If all you want is to see the linear extension of convolving one kernel by another just zero extend both to the sum of the two lengths and then fft each, element by element complex multiply them, and ifft.

[Vincent Picaud]

Hi All,
I have a blog post
https://pixorblog.wordpress.com/2016/07/17/direct-convolution/
that presents a detailed implementation of direct convolution in Julia.
Maybe it can help,
Vincent