Forums

convolution with FFT: same resut with or without zero-padding?

Started by Felipe G. Nievinski April 25, 2009
Hi. Given
    n = 6
    xa = rand(n,1)
    xb = rand(n,1)
If I do
    ya = fft(xa)
    yb = fft(xb)
    z = ya.*yb
    w = ifft(z)
I get the result that I want (w); I'm interested in obtaining exactly
the same result zero-padding input vectors to the next power of two,
for faster FFT (my n is actually ~ 1e6 and I have several such input
vectors); I tried doing
    n2 = 2^nextpow2(n)
    ya2 = fft(xa,n2)
    yb2 = fft(xb,n2)
    z2 = ya2.*yb2
    w2 = ifft(z2, n)
    w2 = ifft(z2(1:n))
    w2 = ifft(z2);  w2 = w2(1:n)
but none of these yields w2 == w. This is supposed to be a piece of
cake...

Thank you for your time.
Felipe.

On Apr 25, 2:05&#2013266080;am, "Bruno Luong" <b.lu...@fogale.findmycountry>
wrote:
> "Felipe G. Nievinski" <fgnievin...@gmail.com> wrote in message <488f6733-cd24-4427-8658-ae1497660...@j9g2000prh.googlegroups.com>... > >> I'm interested in obtaining exactly >> the same result zero-padding input vectors to the next power of two, >> for faster FFT (my n is actually ~ 1e6 and I have several such input >> vectors); > > This is an illusion (of a cake?). >
Let me try to elaborate on that. Given the intermediary results I posted originally, n2 = 2^nextpow2(n) ya2 = fft(xa,n2) yb2 = fft(xb,n2) it's well known that xa2 = ifft(ya2); xa2(n:n2) = [] xb2 = ifft(yb2); xb2(n:n2) = [] yields xa2==xa, xb2==xb. In other words, zero-padding is harmless for the FFT/IFFT of xa & xb individually, and potentially beneficial in terms of their processing speed. I'm looking for a similiar expression for the convolution of xa & xb, taking advantage of faster FFT for input with power of two length (possibly ya2 & yb2?). Is this fair to expect? Thanks, Felipe.
On Apr 25, 2:49&#2013266080;am, Hendrik van der Heijden <h...@gmx.de> wrote:
> Felipe G. Nievinski schrieb: > > ... > This is cyclic convolution, not linear convolution.
Thanks for clarifying the distinction between the two types of convolution.
> If you have a nice optimized FFT (like FFTW), power of two is not > necessarily the fastest transform size. >
I'll certainly bear that in mind; at this point, the time spent trying to get this issue right might not be worth the possible speed gains.
> You need to transform back with size n2, restore the cyclic > boundary conditions (wrap-around at n instead of n2) and > then crop w2 to (1:n). > > For a, b of length n and transform size t, > with zero padding a(n+1:t) = 0, b(n+1:t) = 0: > > lw = ifft(fft(xa,n2).*fft(xb,n2), n2) > > cw = lw(1:n) > cw(1:n2-n) += lw(n+1:n2) > > cw should now be same as w, except for rounding errors. >
Putting everything together: n = 6 rand('seed',0) xa = rand(n,1) xb = rand(n,1) ya = fft(xa) yb = fft(xb) z = ya.*yb w = ifft(z) n2 = 2^nextpow2(n) lw = ifft(fft(xa,n2).*fft(xb,n2),n2) w2 = lw(1:n) w2(1:n2-n) = w2(1:n2-n) + lw(n+1:n2) unfortunately I still get w2 different than w: w = 0.89221 1.0851 1.3708 1.7607 1.4493 1.2075 w2 = 1.8637 1.9155 0.65666 0.93028 1.1919 1.2075 lw = 1.0853 1.0368 0.65666 0.93028 1.1919 1.2075 0.77848 0.87867 I'm sure we're close to getting this right; I just can't seem to get my head around the "cyclic boundary conditions". Would you care to give it another try? Your help is appreciated. Felipe. PS: there seems to be related content at <http://en.wikipedia.org/wiki/ Bluestein's_FFT_algorithm>, starting at "The use of zero-padding for the convolution ... deserves some additional comment."
On Apr 25, 3:42&#2013266080;am, "Bruno Luong" <b.lu...@fogale.findmycountry>
wrote:
> % Try this, it is not newpower2 of n but 2*n to avoid wrap around > n = 6 > xa = rand(n,1) > xb = rand(n,1) > > n2 = 2^nextpow2(2*n); > > w=conv(xa,xb) > > w2 = ifft(fft(xa,n2).* fft(xb,n2),n2); > w2 = w2(1:2*n-1) >
Putting everything together: n = 6 rand('seed',0) xa = rand(n,1) xb = rand(n,1) ya = fft(xa) yb = fft(xb) z = ya.*yb w = ifft(z) %n2 = 2^nextpow2(n) % it is not newpower2 of n but 2*n to avoid wrap around: n2 = 2^nextpow2(2*n) w2 = ifft(fft(xa,n2).* fft(xb,n2),n2) w2 = w2(1:2*n-1) w3 = conv(xa,xb) w3 == w2 works beatifuly, but both are different than my original w. :- ( Furthermore, this would require doubling my already rather large input vectors... :-| My convolution problem seems a little convoluted, indeed! :-D Thanks a lot for your help! Still on the lookout for a definitive solution. Felipe.
Hendrik van der Heijden schrieb:
> You need to transform back with size n2, restore the cyclic > boundary conditions (wrap-around at n instead of n2) and > then crop w2 to (1:n). > > For a, b of length n and transform size t, > with zero padding a(n+1:t) = 0, b(n+1:t) = 0: > > lw = ifft(fft(xa,t).*fft(xb,t), t) > > cw = lw(1:n) > cw(1:t-n) += lw(n+1:t) > > cw should now be same as w, except for rounding errors.
(I replaced n2 with t in the quote above.) I forgot one important thing: lw must not wrap-around, therefore you need n2 >= 2n-1. You already mentioned this in your other post. Are you sure you really want circular (not linear) convolution? Hendrik vdH
Hendrik van der Heijden schrieb:
 > You need to transform back with size n2, restore the cyclic
 > boundary conditions (wrap-around at n instead of n2) and
 > then crop w2 to (1:n).
 >
 > For a, b of length n and transform size t,
 > with zero padding a(n+1:t) = 0, b(n+1:t) = 0:

Some corrections (also replaced n2 with t):

I forgot one important thing: lw must not wrap-around,
therefore you need  t >= 2n-1. You already mentioned
this in your other post.

e = min(n-1, t-n)	% extra elems to wrap around
cw = lw(1:n)
cw(1:e) += lw(n+1: n+1+e)


Are you sure you really want circular convolution?
Usually, one does it the other way round: zero-pad a
circular convolution to make it work like a linear one.


Hendrik vdH
Hi. Thanks for giving it another try.

On Apr 25, 4:32&#2013266080;am, Hendrik van der Heijden <h...@gmx.de> wrote:
> Hendrik van der Heijden schrieb: > &#2013266080;> You need to transform back with size n2, restore the cyclic > &#2013266080;> boundary conditions (wrap-around at n instead of n2) and > &#2013266080;> then crop w2 to (1:n). > &#2013266080;> > &#2013266080;> For a, b of length n and transform size t, > &#2013266080;> with zero padding a(n+1:t) = 0, b(n+1:t) = 0: > > Some corrections (also replaced n2 with t): > > I forgot one important thing: lw must not wrap-around, > therefore you need &#2013266080;t >= 2n-1. You already mentioned > this in your other post. > > e = min(n-1, t-n) &#2013266080; &#2013266080; &#2013266080; % extra elems to wrap around > cw = lw(1:n) > cw(1:e) += lw(n+1: n+1+e) >
> Are you sure you really want circular convolution? > Usually, one does it the other way round: zero-pad a > circular convolution to make it work like a linear one. >
I just had to make one minor modification: n = 6 rand('seed',0) xa = rand(n,1) xb = rand(n,1) ya = fft(xa) yb = fft(xb) z = ya.*yb w = ifft(z) %n2 = 2^nextpow2(n) % it is not newpower2 of n but 2*n to avoid wrap around: n2 = 2^nextpow2(n + n - 1) % extra elems to wrap around: ne = min(n-1,n2-n) lw = ifft(fft(xa,n2).*fft(xb,n2),n2) w2 = lw(1:n) %w2(1:ne) = w2(1:ne) + lw(n+1:n+1+ne) % (slightly wrong). w2(1:ne) = w2(1:ne) + lw(n+1:n+ne) Now w2 == w, thanks!!! Given the requirement to do an FFT on a vector twice as large as my original data, and also the fact that FFTw (used in MATLAB) is not necessarily slower for non-power of two input, I think I'll stick with my original data and give up on zero-padding it. I thought that zero- padding for convolution would be trivial, but based on this thread it seems fair to say that it's tricky. Thanks once again, Felipe.