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.

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

Started by ●April 25, 2009

Reply by ●April 25, 20092009-04-25

On Apr 25, 2:05�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.

Reply by ●April 25, 20092009-04-25

On Apr 25, 2:49�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."

Reply by ●April 25, 20092009-04-25

On Apr 25, 3:42�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.

Reply by ●April 25, 20092009-04-25

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

Reply by ●April 25, 20092009-04-25

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

Reply by ●April 25, 20092009-04-25

Hi. Thanks for giving it another try. On Apr 25, 4:32�am, Hendrik van der Heijden <h...@gmx.de> wrote:> 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. >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.