Filter löschen
Filter löschen

Error in convolution of two gaussians using FFT and IFFT

3 Ansichten (letzte 30 Tage)
Adam
Adam am 10 Okt. 2014
Bearbeitet: Matt J am 17 Okt. 2014
I am trying two convolve two gaussian functions with the form of f(x)=a*exp(-(x-b)^2/2*c^2). a is the height, b is the position of the curve's center, and c controls the width of the "bell" shape.
To perform a convolution, you can take the inverse fourier transform of two functions that undergo a fourier transform and are then multiplied: F^-1(F(f1(x))*F(f2(x)))
My code is as follows:
xstep=-10:0.05:25;
c = 5/(2*sqrt(2*log(2)));
c2 = .2/(2*sqrt(2*log(2)));
y1 = exp(-(xstep-5).^2/(2*c^2)); %actual
y2 = .5*exp(-(xstep).^2/(2*c2^2)); %instrument response
Fy1 = fft(y1);
Fy2 = fft(y2);
con = ifft(Fy1.*Fy2); % This is my convolution
figure(2)
hold on
plot(xstep,y1)
plot(xstep,y2,'r')
plot(xstep,con,'--k')
axis([-10,25,0,2])
hold off
This produces a function that is shifted and much larger than it should be. The correct results should appear as such:
test = (1/(sqrt(2*pi*(c^2+c2^2))))*exp((-(xstep-(5+3)).^2)/(2*(c^2+c2^2)));
plot(xstep,test)

Antworten (2)

SK
SK am 11 Okt. 2014
xstep starts from -10, so the peak of y1 is actually at 15 and the peak of y2 is at 10 (fft() doesn't know about your x-axis). The convolution has a peak at 10+15 = 25. But the x-axis on your plot starts at -10, ie: 0 is -10, so 25 will be 15. Hence your graph.

Matt J
Matt J am 11 Okt. 2014
Bearbeitet: Matt J am 11 Okt. 2014
You forgot to weight your (discrete) convolution by the sampling interval
con = ifft(Fy1.*Fy2)*0.05;
which partly account for the scaling. You also haven't weighted y1 and y2 by 1/sqrt(2*pi*c^2) and1/sqrt(2*pi*c2^2). So there is no reason I can see that the multiplier in con will be (1/(sqrt(2*pi*(c^2+c2^2)))).
The shifts are fine, since shifts in convolution are addiitive. Since y1 is shifted by 10 from the beginning of the array and y2 is shifted by 15 from the beginning of the array, then con will be shifted by 10+15=25 from the beginning of the array, which is what your plots show.
  2 Kommentare
Adam
Adam am 17 Okt. 2014
I have found that you are correct. My assumption of the "test" value was that y1 and y2 were scaled by 1/(sqrt(2pi)*sigma). However, they were scaled by 1 and 0.5. You mentioned that the final result should be scaled by 0.05. Is there anything else that should contribute to the final scaling?
Matt J
Matt J am 17 Okt. 2014
Bearbeitet: Matt J am 17 Okt. 2014
Well, you won't be able to get the exact result of a continuous convolution. You will lose accuracy (scaling and otherwise) due to the discretization of the signals. That should be negligible, however, if your sampling is fine enough.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Fourier Analysis and Filtering finden Sie in Help Center und File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by