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

Started by 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...

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.

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

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.
```