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

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&#4294967295;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&#4294967295;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&#4294967295;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&#4294967295;am, Hendrik van der Heijden <h...@gmx.de> wrote:
> Hendrik van der Heijden schrieb:
> &#4294967295;> You need to transform back with size n2, restore the cyclic
> &#4294967295;> boundary conditions (wrap-around at n instead of n2) and
> &#4294967295;> then crop w2 to (1:n).
> &#4294967295;>
> &#4294967295;> For a, b of length n and transform size t,
> &#4294967295;> 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 &#4294967295;t >= 2n-1. You already mentioned
> this in your other post.
>
> e = min(n-1, t-n) &#4294967295; &#4294967295; &#4294967295; % 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.