PDA

View Full Version : convolution output length



hcrisp
02-24-2011, 08:36 AM
I am using CONVOL1D to calculate the autocorrelation of a signal. I see there is a keyword DIRECT which overrides the default method (FFT) to compute the convolution directly. However the output vector lengths are not the same.



x = randomu(s, 512)
y1 = convol1d(x, x)
y2 = convol1d(x, x, /direct)
info, y1, y2
; Y1 (DWPLOTINTERPDATA) DOUBLE = Array(512)
; Y2 (DWPLOTINTERPDATA) DOUBLE = Array(511)

x = randomu(s, 5437)
y1 = convol1d(x, x)
y2 = convol1d(x, x, /direct)
info, y1, y2
; Y1 (DWPLOTINTERPDATA) DOUBLE = Array(10935)
; Y2 (DWPLOTINTERPDATA) DOUBLE = Array(10873)

The DIRECT method seems to give the correct length, while it appears that the FFT method is adding extra zeros to the end (probably due to POW235 padding). How can I reliably remove these extra zeros? Something like this?



npts = 5437
x = randomu(s, npts)
y1 = convol1d(x, x)
cnpts = (npts * 2) - 1
y1 = y1(0:cnpts-1)
info, y1
; Y1 (DWPLOTINTERPDATA) DOUBLE = Array(10873)

hcrisp
02-24-2011, 10:58 AM
CORR1D has related behavior, but the zeros are in the middle not at the end. How do I deal with this?



x1 = randomn(s, 4096)
y1 = CORR1D(x1,x1)
x2 = [x1, 0.]
y2 = CORR1D(x2,x2)
info, y1, y2
; Y1 (DWPLOTINTERPDATA) FLOAT = Array(8192)
; Y2 (DWPLOTINTERPDATA) FLOAT = Array(8640)
plot, y1
oplot, y2, color=230

omega
03-03-2011, 07:33 AM
Looking at the files convol1d.pro and corr1d.pro in the lib directory of the math option may shed some light on the implementation details. They appear to contain detailed comments on what is computed and how.

hth
-w

hcrisp
03-03-2011, 11:29 AM
That I have done already, and it looks like the developer (Mike Pulverenti, 5/6/94) left out the portion of code which would have removed the padded zeros when the FFTCOMP method is used. Thus, unless VNI modifies the file, you must do it manually.

For CORR1D I have found this to work:



x = randomn(s, 4097)
acorr = corr1d(x, x)
npts = N_ELEMENTS(x)
cnpts = (npts * 2) - 1
nonzero_ids = [LINDGEN(npts), REVERSE(N_ELEMENTS(acorr) - LINDGEN(npts-1) - 1)]
acorr = acorr(nonzero_ids)