Convolution property of DFT
in [DSP]
|
From: third_person on 29 Jun 2010 12:26 Hi, I'm trying fast convolution property but there seems to be some mistake (with the answer). Here is the Matlab code for it, clear all; clc; % Test Vector Convolution a = [1 2 3 4 5]; b = [10 20 30 40 50]; c= conv(a,b) A = [1 2 3 4 5 zeros(1,5)]; B = [10 20 30 40 50 zeros(1,5)]; d= ifft(fft(A) .* fft(B)) c - d(1:9) The results are: c = 10 40 100 200 350 440 460 400 250 d = 10.0000 40.0000 100.0000 200.0000 350.0000 440.0000 460.0000 400.0000 250.0000 0.0000 ans = (for the difference b/w the two) 1.0e-012 * -0.0906 -0.0497 -0.0284 -0.0284 0 0.1137 0 0 0 Shouldn't the result be smaller than eps (2.2204e-016)?
From: Raymond Toy on 29 Jun 2010 14:17 On 6/29/10 12:26 PM, third_person wrote: > Hi, I'm trying fast convolution property but there seems to be some mistake > (with the answer). > > Here is the Matlab code for it, > > clear all; clc; > > % Test Vector Convolution > > a = [1 2 3 4 5]; > b = [10 20 30 40 50]; > c= conv(a,b) > A = [1 2 3 4 5 zeros(1,5)]; > B = [10 20 30 40 50 zeros(1,5)]; > d= ifft(fft(A) .* fft(B)) > c - d(1:9) > > The results are: > > c = > > 10 40 100 200 350 440 460 400 250 > > d = > > 10.0000 40.0000 100.0000 200.0000 350.0000 440.0000 460.0000 > 400.0000 250.0000 0.0000 > > ans = (for the difference b/w the two) > > 1.0e-012 * > > -0.0906 -0.0497 -0.0284 -0.0284 0 0.1137 0 0 0 > > > Shouldn't the result be smaller than eps (2.2204e-016)? > Why do you think it should smaller than eps? Do you think fft and ifft have no roundoff? I don't know what the actual roundoff should be but a difference of 1e-13 seems fairly reasonable. Ray
From: Greg Heath on 30 Jun 2010 09:53 On Jun 29, 12:26 pm, "third_person" <third_person(a)n_o_s_p_a_m.ymail.com> wrote: > Hi, I'm trying fast convolution property but > there seems to be some mistake > (with the answer). > > Here is the Matlab code for it, > > clear all; clc; > > % Test Vector Convolution > > a = [1 2 3 4 5]; > b = [10 20 30 40 50]; > c= conv(a,b) > A = [1 2 3 4 5 zeros(1,5)]; > B = [10 20 30 40 50 zeros(1,5)]; > d= ifft(fft(A) .* fft(B)) > c - d(1:9) > > The results are: > > c = > > 10 40 100 200 350 440 460 400 250 > > d = > > 10.0000 40.0000 100.0000 200.0000 350.0000 440.0000 460.0000 > 400.0000 250.0000 0.0000 > > ans = (for the difference b/w the two) > > 1.0e-012 * > > -0.0906 -0.0497 -0.0284 -0.0284 0 0.1137 0 0 0 > > Shouldn't the result be smaller than eps (2.2204e-016)? Typically, 1. A and B are zeropadded to length(A)+length(B)-1 2. If A an B are real, d = real(ifft(fft(A) .* fft(B)) is used because ifft is notorious for creating spurious imaginary roundoff error However, the results I obtained below surprised me (Note the change in notation) clear all, clc a = [1 2 3 4 5]'; b = [10 20 30 40 50]'; c= conv(a,b) ; A = [1 2 3 4 5 zeros(1,4)]'; B = [10 20 30 40 50 zeros(1,4)]'; C = ifft(fft(A) .* fft(B)); D = [c C(1:9)] % D = % % 10 10 % 40 40 +1.2632e-014i % 100 100 -6.3159e-015i % 200 200 % 350 350 -6.3159e-015i % 440 440 -7.7816e-015i % 460 460 % 400 400 -6.3159e-015i % 250 250 +1.4097e-014i A = [1 2 3 4 5 zeros(1,5)]'; B = [10 20 30 40 50 zeros(1,5)]'; C = ifft(fft(A) .* fft(B)); D = [c C(1:9)] % D = % % 10 10 % 40 40 % 100 100 % 200 200 % 350 350 % 440 440 % 460 460 % 400 400 % 250 250 I can't explain it. Can someone else? Greg
From: robert bristow-johnson on 30 Jun 2010 12:03 On Jun 30, 9:53 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Jun 29, 12:26 pm, "third_person" > > > > <third_person(a)n_o_s_p_a_m.ymail.com> wrote: > > Hi, I'm trying fast convolution property but > > there seems to be some mistake > > (with the answer). > > > Here is the Matlab code for it, > > > clear all; clc; > > > % Test Vector Convolution > > > a = [1 2 3 4 5]; > > b = [10 20 30 40 50]; > > c= conv(a,b) > > A = [1 2 3 4 5 zeros(1,5)]; > > B = [10 20 30 40 50 zeros(1,5)]; > > d= ifft(fft(A) .* fft(B)) > > c - d(1:9) > > > The results are: > > > c = > > > 10 40 100 200 350 440 460 400 250 > > > d = > > > 10.0000 40.0000 100.0000 200.0000 350.0000 440.0000 460.0000 > > 400.0000 250.0000 0.0000 > > > ans = (for the difference b/w the two) > > > 1.0e-012 * > > > -0.0906 -0.0497 -0.0284 -0.0284 0 0.1137 0 0 0 > > > Shouldn't the result be smaller than eps (2.2204e-016)? > > Typically, > 1. A and B are zeropadded to length(A)+length(B)-1 > 2. If A an B are real, d = real(ifft(fft(A) .* fft(B)) > is used because ifft is notorious for creating > spurious imaginary roundoff error > > However, the results I obtained below surprised me > (Note the change in notation) > > clear all, clc > > a = [1 2 3 4 5]'; > b = [10 20 30 40 50]'; > c= conv(a,b) ; > A = [1 2 3 4 5 zeros(1,4)]'; > B = [10 20 30 40 50 zeros(1,4)]'; > C = ifft(fft(A) .* fft(B)); > D = [c C(1:9)] > > % D = > % > % 10 10 > % 40 40 +1.2632e-014i > % 100 100 -6.3159e-015i > % 200 200 > % 350 350 -6.3159e-015i > % 440 440 -7.7816e-015i > % 460 460 > % 400 400 -6.3159e-015i > % 250 250 +1.4097e-014i > > A = [1 2 3 4 5 zeros(1,5)]'; > B = [10 20 30 40 50 zeros(1,5)]'; > C = ifft(fft(A) .* fft(B)); > D = [c C(1:9)] > > % D = > % > % 10 10 > % 40 40 > % 100 100 > % 200 200 > % 350 350 > % 440 440 > % 460 460 > % 400 400 > % 250 250 > > I can't explain it. Can someone else? i don't understand what's troubling you, Greg. is it the extremely tiny imaginary values that result (presumably from roundoff) when N=9 that don't when N=10? r b-j
From: Greg Heath on 30 Jun 2010 16:10 On Jun 30, 12:03 pm, robert bristow-johnson <r...(a)audioimagination.com> wrote: > On Jun 30, 9:53 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > > > > > > > On Jun 29, 12:26 pm, "third_person" > > > <third_person(a)n_o_s_p_a_m.ymail.com> wrote: > > > Hi, I'm trying fast convolution property but > > > there seems to be some mistake > > > (with the answer). > > > > Here is the Matlab code for it, > > > > clear all; clc; > > > > % Test Vector Convolution > > > > a = [1 2 3 4 5]; > > > b = [10 20 30 40 50]; > > > c= conv(a,b) > > > A = [1 2 3 4 5 zeros(1,5)]; > > > B = [10 20 30 40 50 zeros(1,5)]; > > > d= ifft(fft(A) .* fft(B)) > > > c - d(1:9) > > > > The results are: > > > > c = > > > > 10 40 100 200 350 440 460 400 250 > > > > d = > > > > 10.0000 40.0000 100.0000 200.0000 350.0000 440..0000 460.0000 > > > 400.0000 250.0000 0.0000 > > > > ans = (for the difference b/w the two) > > > > 1.0e-012 * > > > > -0.0906 -0.0497 -0.0284 -0.0284 0 0.1137 0 0 0 > > > > Shouldn't the result be smaller than eps (2.2204e-016)? > > > Typically, > > 1. A and B are zeropadded to length(A)+length(B)-1 > > 2. If A an B are real, d = real(ifft(fft(A) .* fft(B)) > > is used because ifft is notorious for creating > > spurious imaginary roundoff error > > > However, the results I obtained below surprised me > > (Note the change in notation) > > > clear all, clc > > > a = [1 2 3 4 5]'; > > b = [10 20 30 40 50]'; > > c= conv(a,b) ; > > A = [1 2 3 4 5 zeros(1,4)]'; > > B = [10 20 30 40 50 zeros(1,4)]'; > > C = ifft(fft(A) .* fft(B)); > > D = [c C(1:9)] > > > % D = > > % > > % 10 10 > > % 40 40 +1.2632e-014i > > % 100 100 -6.3159e-015i > > % 200 200 > > % 350 350 -6.3159e-015i > > % 440 440 -7.7816e-015i > > % 460 460 > > % 400 400 -6.3159e-015i > > % 250 250 +1.4097e-014i > > > A = [1 2 3 4 5 zeros(1,5)]'; > > B = [10 20 30 40 50 zeros(1,5)]'; > > C = ifft(fft(A) .* fft(B)); > > D = [c C(1:9)] > > > % D = > > % > > % 10 10 > > % 40 40 > > % 100 100 > > % 200 200 > > % 350 350 > > % 440 440 > > % 460 460 > > % 400 400 > > % 250 250 > > > I can't explain it. Can someone else? > > i don't understand what's troubling you, Greg. is it the extremely > tiny imaginary values that result (presumably from roundoff) when > N=9 that don't when N=10? I usually get imaginary valued roundoff when I use ifft and the result should be real. I was intrigued that, in contrast, the OP got purely real roundoff. Then I noticed that he used one more zero than necessary in the zeropadding. So, I removed the extra zero and got purely imaginary roundoff. Satisfied that my understanding was validated, I put the extra zero back in to see if that was the cause of the real valued roundoff... Much to my surprise, my calculation resulted in no roundoff error. I find this puzzling, even intriguing, but certainly not troublesome. Greg
|
Pages: 1 Prev: Sigma Delta - Decimation Filter Design Next: Noise Normalisation |