COMPGROUPS.NET | Search | Post Question | Groups | Stream | About | Register

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

• Email
• Follow

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

```
 0
Reply fgnievinski1 (40) 4/25/2009 7:49:57 AM

See related articles to this posting

```"Felipe G. Nievinski" <fgnievinski@gmail.com> wrote in message <488f6733-cd24-4427-8658-ae1497660ccd@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?).

Bruno
```
 0
Reply b.luong5955 (6403) 4/25/2009 8:05:06 AM

```On Apr 25, 2:05=A0am, "Bruno Luong" <b.lu...@fogale.findmycountry>
wrote:
> "Felipe G. Nievinski" <fgnievin...@gmail.com> wrote in message <488f6733-=
>
>> 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 =3D 2^nextpow2(n)
ya2 =3D fft(xa,n2)
yb2 =3D fft(xb,n2)
it's well known that
xa2 =3D ifft(ya2);  xa2(n:n2) =3D []
xb2 =3D ifft(yb2);  xb2(n:n2) =3D []
yields xa2=3D=3Dxa, xb2=3D=3Dxb. In other words, zero-padding is harmless f=
or
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.
```
 0
Reply fgnievinski1 (40) 4/25/2009 8:49:06 AM

```Felipe G. Nievinski schrieb:
> 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)

This is cyclic convolution, not linear convolution.
Linear convolution of xa (lenght n) and xb (length n)
would yield a lenght 2n-1 result.

> 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);

If you have a nice optimized FFT (like FFTW), power of two is not
necessarily the fastest transform size.

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

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.

Hendrik vdH
```
 0
Reply hvdh (28) 4/25/2009 8:49:45 AM

```"Felipe G. Nievinski" <fgnievinski@gmail.com> wrote in message <973fbe3c-6067-4d9f-8530-0f888e5d3a17@s38g2000prg.googlegroups.com>...
xb2=3D=3Dxb. In other words, zero-padding is harmless f=
> or
> 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?

Quite fair. What I meant is there is a bunch of details you need to take care in a discrete case to make it works, and it is not a "piece of cake" (not to me anyway).

IMO the padding need more than the next power two to avoid wrap around because both length of vector adds. Anyway I have not work on the detail, so I will stop my comment here.

Certainly someone must done this already (I believe one of Matlab functions CONV, CONV2 or CONVN does just that, can't recall which one)

Good luck

Bruno
```
 0
Reply b.luong5955 (6403) 4/25/2009 9:30:16 AM

```% 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)

% Bruno
```
 0
Reply b.luong5955 (6403) 4/25/2009 9:42:02 AM

```On Apr 25, 2:49=A0am, 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) =3D 0, b(n+1:t) =3D 0:
>
> lw =3D ifft(fft(xa,n2).*fft(xb,n2), n2)
>
> cw =3D lw(1:n)
> cw(1:n2-n) +=3D lw(n+1:n2)
>
> cw should now be same as w, except for rounding errors.
>
Putting everything together:

n =3D 6
rand('seed',0)
xa =3D rand(n,1)
xb =3D rand(n,1)

ya =3D fft(xa)
yb =3D fft(xb)
z =3D ya.*yb
w =3D ifft(z)

n2 =3D 2^nextpow2(n)
lw =3D ifft(fft(xa,n2).*fft(xb,n2),n2)
w2 =3D lw(1:n)
w2(1:n2-n) =3D w2(1:n2-n) + lw(n+1:n2)

unfortunately I still get w2 different than w:

w =3D
0.89221
1.0851
1.3708
1.7607
1.4493
1.2075

w2 =3D
1.8637
1.9155
0.65666
0.93028
1.1919
1.2075

lw =3D
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."

```
 0
Reply fgnievinski1 (40) 4/25/2009 10:08:27 AM

```On Apr 25, 3:42=A0am, "Bruno Luong" <b.lu...@fogale.findmycountry>
wrote:
> % Try this, it is not newpower2 of n but 2*n to avoid wrap around
> n =3D 6
> xa =3D rand(n,1)
> xb =3D rand(n,1)
>
> n2 =3D 2^nextpow2(2*n);
>
> w=3Dconv(xa,xb)
>
> w2 =3D ifft(fft(xa,n2).* fft(xb,n2),n2);
> w2 =3D w2(1:2*n-1)
>
Putting everything together:

n =3D 6
rand('seed',0)
xa =3D rand(n,1)
xb =3D rand(n,1)

ya =3D fft(xa)
yb =3D fft(xb)
z =3D ya.*yb
w =3D ifft(z)

%n2 =3D 2^nextpow2(n)
% it is not newpower2 of n but 2*n to avoid wrap around:
n2 =3D 2^nextpow2(2*n)
w2 =3D ifft(fft(xa,n2).* fft(xb,n2),n2)
w2 =3D w2(1:2*n-1)

w3 =3D conv(xa,xb)

w3 =3D=3D 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.
```
 0
Reply fgnievinski1 (40) 4/25/2009 10:10:36 AM

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

```
 0
Reply hvdh (28) 4/25/2009 10:22:29 AM

```This "wraparound" occurs because the convolution Fourier formula in the discrete case applied for periodic function. If the zero padding is not large enough, the head of B (after flipping) get back (wraparound) and multiplied with the head of A when flipped-B vector will be shifted during the convolution.

It seems you can get away be doing TWO convolutions with two flipped vector (on A, or on B exclusively). But then why not pad enough zeros and let FFT do its magic in a single shot?

But again, I'm not expert on this topic, may be someone find a clever method that you are after.

Bruno
```
 0
Reply b.luong5955 (6403) 4/25/2009 10:31:02 AM

```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
```
 0
Reply hvdh (28) 4/25/2009 10:32:20 AM

```Hi. Thanks for giving it another try.

On Apr 25, 4:32=A0am, Hendrik van der Heijden <h...@gmx.de> wrote:
> Hendrik van der Heijden schrieb:
> =A0> You need to transform back with size n2, restore the cyclic
> =A0> boundary conditions (wrap-around at n instead of n2) and
> =A0> then crop w2 to (1:n).
> =A0>
> =A0> For a, b of length n and transform size t,
> =A0> with zero padding a(n+1:t) =3D 0, b(n+1:t) =3D 0:
>
> Some corrections (also replaced n2 with t):
>
> I forgot one important thing: lw must not wrap-around,
> therefore you need =A0t >=3D 2n-1. You already mentioned
> this in your other post.
>
> e =3D min(n-1, t-n) =A0 =A0 =A0 % extra elems to wrap around
> cw =3D lw(1:n)
> cw(1:e) +=3D 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 =3D 6
rand('seed',0)
xa =3D rand(n,1)
xb =3D rand(n,1)

ya =3D fft(xa)
yb =3D fft(xb)
z =3D ya.*yb
w =3D ifft(z)

%n2 =3D 2^nextpow2(n)
% it is not newpower2 of n but 2*n to avoid wrap around:
n2 =3D 2^nextpow2(n + n - 1)
% extra elems to wrap around:
ne =3D min(n-1,n2-n)

lw =3D ifft(fft(xa,n2).*fft(xb,n2),n2)
w2 =3D lw(1:n)
%w2(1:ne) =3D w2(1:ne) + lw(n+1:n+1+ne) % (slightly wrong).
w2(1:ne) =3D w2(1:ne) + lw(n+1:n+ne)

Now w2 =3D=3D 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.
```
 0
Reply fgnievinski1 (40) 4/25/2009 5:26:29 PM

11 Replies
32 Views

Similar Articles

12/7/2013 6:36:11 PM
page loaded in 177835 ms. (0)