Reply by Felipe G. Nievinski●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.
Reply by Hendrik van der Heijden●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 Hendrik van der Heijden●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 Felipe G. Nievinski●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 Felipe G. Nievinski●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 Felipe G. Nievinski●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 Felipe G. Nievinski●April 25, 20092009-04-25
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.