why my fft does not match my convolution ?

Question

Rabih Sokhen 2022-3-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1662500-why-my-fft-does-not-match-my-convolution

编辑： Paul 2022-3-9

采纳的回答： Paul

在 MATLAB Online 中打开

hy guys, i am trying to verify that g(w)=fft( h_t . x_t) = conv( h_f,x_f)

h_t is my impulse response in the time domain

x_t is my filter in the time domain

h_f and x_f are the Fourier transform of h_t and x_t

I have tried to write the following code, however, I don't have the same value.

Any suggestion how to fix it?

thank you in advance

code:

clear all
clc
dt=0.01; 
t=-15:dt:15; 
df = 1/dt; 
f = linspace(-df/2, df/2, length(t));
w=2*pi*f;
w0=10;
Q=5;
a=w0/sqrt(1-1/(4*Q^2));
b=w0*sqrt(1-1/(4*Q^2));
h_t=a*sin(b*t).*exp(-w0*t/(2*Q)).*(t>=0);
h_f1=fftshift(fft(ifftshift(h_t)))*dt;
figure
subplot(231)
plot(t,real(h_t),'r')
grid on
title('h_t')
xlabel('temps');
legend('h_t=w_0/(1- (1/4Q^2)^0^.^5) sin (w_0(1- (1/4Q^2)^0^.^5 t) e^-^w^0^t^/^2^Q')
subplot(234)
plot(f,real(h_f1))
legend('h_f= TF(h_t)');
grid on
title('h_f=TF (h_t)')
ylabel('amplitude')
xlabel('frequence');
xlim([-10 10])
subplot(232)
dnn=1/2.355 *4; %largeur du faisceau gaussien
porte_t=exp(-0.5*(t./dnn).^20);
plot(t,real(porte_t),'r')
grid on
title('porte_t=hyper gaussien')
xlabel('temps');
subplot(235)
porte_f=(fftshift(fft(ifftshift(porte_t))))*dt;
plot(f,real(porte_f),'-b')
grid on
xlim([-5 5]);
title('porte_f= TF( porte_t)')
xlabel('frequence');
subplot(233)
g_t=h_t.*porte_t;
plot(t,real(g_t),'r')
title('g_t=h_t . porte_t')
grid on
xlabel('temps');
subplot(236)
g_f1=fftshift(fft(ifftshift(g_t),(length(h_f1)+length(porte_f)-1)))*dt;
g_f2=conv(h_f1,porte_f);
ff = linspace(-df/2, df/2, length(g_f1));
plot(ff,real(g_f1),'k',ff,real(g_f2),'b')
grid on
xlabel('frequence');
legend('g_f_1=TF(g_t)','g_f_2=conv(h_f, porte_f)')
xlim([-15 15])

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2022-3-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1662500-why-my-fft-does-not-match-my-convolution#answer_909369

在 MATLAB Online 中打开

The DFT of the product of two finite duration sequences is the normalized circular convolution of their DFTs. For example:

x = 1:5;
y = 5:-1:1;
X = fft(x);
Y = fft(y);
XY = fft(x.*y)
XY = 
  35.0000 + 0.0000i  -4.7361 - 3.4410i  -0.2639 - 0.8123i  -0.2639 + 0.8123i  -4.7361 + 3.4410i
cconv(X,Y,5)/5
ans = 
  35.0000 + 0.0000i  -4.7361 - 3.4410i  -0.2639 - 0.8123i  -0.2639 + 0.8123i  -4.7361 + 3.4410i

Is this what you're trying to show?

8 个评论
显示 6更早的评论隐藏 6更早的评论

Rabih Sokhen 2022-3-8

hy Paul

i multiply by (dt) to normalis

and i use ( ifftshift )bcz the time in symetric to the referenc zero

my code is :

code:

clear all

clc

dt=0.01;

t=-15:dt:15;

df = 1/dt;

f = linspace(-df/2, df/2, length(t));

w=2*pi*f;

w0=10;

Q=5;

a=w0/sqrt(1-1/(4*Q^2));

b=w0*sqrt(1-1/(4*Q^2));

h_t=a*sin(b*t).*exp(-w0*t/(2*Q)).*(t>=0);

h_f1=fftshift(fft(ifftshift(h_t)))*dt;

figure

subplot(231)

plot(t,real(h_t),'r')

grid on

title('h_t')

xlabel('temps');

legend('h_t=w_0/(1- (1/4Q^2)^0^.^5) sin (w_0(1- (1/4Q^2)^0^.^5 t) e^-^w^0^t^/^2^Q')

subplot(234)

plot(f,real(h_f1))

legend('h_f= TF(h_t)');

grid on

title('h_f=TF (h_t)')

ylabel('amplitude')

xlabel('frequence');

xlim([-10 10])

subplot(232)

dnn=1/2.355 *4; %largeur du faisceau gaussien

porte_t=exp(-0.5*(t./dnn).^20);

plot(t,real(porte_t),'r')

grid on

title('porte_t=hyper gaussien')

xlabel('temps');

subplot(235)

porte_f=(fftshift(fft(ifftshift(porte_t))))*dt;

plot(f,real(porte_f),'-b')

grid on

xlim([-5 5]);

title('porte_f= TF( porte_t)')

xlabel('frequence');

subplot(233)

g_t=h_t.*porte_t;

plot(t,real(g_t),'r')

title('g_t=h_t . porte_t')

grid on

xlabel('temps');

subplot(236)

g_f1=fftshift(fft(ifftshift(g_t)))*dt;

g_f2=fftshift(cconv(h_f1,porte_f,length(g_f1))/length(g_f1));

ff = linspace(-df/2, df/2, length(g_f1));

plot(ff,real(g_f1),'k',ff,real(g_f2),'b')

grid on

xlabel('frequence');

legend('g_f_1=TF(g_t)','g_f_2=conv(h_f, porte_f)')

xlim([-15 15])

Paul 2022-3-8

编辑：Paul 2022-3-9

在 MATLAB Online 中打开

Ok. One more clarification before proceeding. In the code above, the signal x_t is called porte_t and is defined as

porte_t=exp(-0.5*(t./dnn).^20);

and h_t is defined as

h_t=a*sin(b*t).*exp(-w0*t/(2*Q)).*(t>=0);

with both comptuet for values of t as:

t=-15:dt:15;

Are the values of t defined that way because porte_t = 0 for abs(t) > 15 and h_t = 0 for t > 15? Or do the definitions for porte_t and h_t apply for -inf < t < inf, but -15 < t < 15 is assumed sufficient for the analysis?

Also, if h(t) is the impulse response of the filter and x(t) is the input to the filter, the output of the filter g(t) is given by

g(t) = convolution(h(t),x(t))

Here, convolution() means the convolution integral. It can be approximated with the Matlab function conv() and using properties convolution.

In the frequency domain the Fourier transform of g is given by

G(w) = product(H(w),X(w)), i.e., the multiplication of the Fourier transforms of h(t) and x(t). This operation can also be approximated using the fft() function, again being careful to use fft() properly for this operation and assumng that both x(t) and h(t) are non-zero over a finite duration, or if you're willing to approximate them as such.

Paul 2022-3-9

编辑：Paul 2022-3-9

在 MATLAB Online 中打开

Here is an example that illustrates the problem as I understand for different functions h(t) and x(t) than in the Question. Maybe it can be adpated to the actual qeustion

Let the filter impulse response be

syms t
h(t) = sin(5*t)*exp(-t)*heaviside(t);

Let the filter input be

x(t) = exp(-t^2);

figure;

fplot([h(t) x(t)],[-5 10])

The output of the filter is given by the convolution integral

syms tau

g(t) = int(x(tau)*h(t-tau),tau,-inf,inf)

g(t) =

We see that Matlab cannot find a closed form expression for g(t). However, we see that the upper bound on the integral is t, because h(t) = 0 for t < 0. We can instead compute the integral numerically.

g(t) = vpaintegral(x(tau)*h(t-tau),tau,-inf,t)

g(t) =

Plot g(t) at some points

figure

plot(-4:.1:5,double(g(-4:.1:5)))

Now compute g(t) using conv() function. In order to do this we need to account for the fact that conv() assumes that the first point in each sequence is at the orgin. But x(t) is non-zero for for t<0. So we will have to shift x(t) to the right, use conv(), and then shift the result back to the left. Looking at the graphs of h(t) and x(t), we will assume that x(t) is zero for t<3 and that h(t) is zero for t>6. Define functions for h(t) and x(t)

hfunc = @(t) (sin(5*t).*exp(-t).*(t>=0));
xfunc = @(t) (exp(-t.^2));

Evaluate samples of h(t) and x(t-3) at samples of t

dt = 0.01;
tval = 0:dt:6;
hval = hfunc(tval);
xvalshifted = xfunc(tval-3);

Now take the convolution sum

gvalshifted = conv(hval,xvalshifted);

Because gval is the convolution sum and we want to approximate the convolution integral, we have to multiply by dt

gvalshifted = gvalshifted*dt;

Now, the first point of gvalshifted corresponds to t = 0. But because gvalshifted was computed with the right-shifted version of x(n*dt), we know from the time-shift property that the convolution of x(nt) and h(nt) takes the same values as gvalshifted, but the associated time vector is left-shifted by the same amount.

gvalt = gvalshifted;

Now compare to the convolution integral

figure

plot(-4:.1:5,double(g(-4:.1:5)),-3 + (0:numel(gvalt)-1)*dt,gvalt,'o')

That plot shows that the scaled convolution sum well-approximates the "exact" result based on the convolution integral.

In the frequency domain, it is true the the Fourier transform of g(t) is product of the Fourier transforms of h(t) and x(t). Even though we can get expressions for H(w) and X(w), Matlab can't find the closed form expression for g(t).

syms w

H(w) = fourier(h(t),t,w)

H(w) =

X(w) = fourier(x(t),t,w)

X(w) =

G(w) = H(w)*X(w);

ifourier(G(w),w,t)

ans =

So again we resort to the discrete-time domain. One approach would be to sample h(t) and x(t), compute the Discrete Time Fourier Transform (DTFT) of each, multiply those together, and then take iDTFT. A simpler approach is to do the same operation with the DFT, as computed by fft(). Again, we'll have shift x(t). More importanty, we have to zero-pad the samples of h(nT) and x(nT) so that the product of their fft's represents linear convolution, not circular convolution, of the underlying sequences. Again, the result has be to scaled by dt to approximate the continuous time result.

nfft = numel(xvalshifted) + numel(hval) - 1;
gvalf = ifft(fft(xvalshifted,nfft).*fft(hval,nfft))*dt;

Plot and verify the result

figure

plot(-4:.1:5,double(g(-4:.1:5)),-3 + (0:numel(gvalf)-1)*dt,gvalf,'o')

As shown, the continuous time convolution integral of h(t) and x(t) can be approximated by the conv() of samples of h(t) and x(t), or by the iDFT of the product of the zero-padded DFTs of samples of h(t) and x(t), as long as h(t) and x(t) are both of finite duration or can be approximated as such.

Rabih Sokhen 2022-3-9

the values of t are defined that way because the definitions for porte_t and h_t apply for -inf < t < inf, but -15 < t < 15 is assumed sufficient for the analysis

thank you Paul, i understand more now

请先登录，再进行评论。

why my fft does not match my convolution ?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

why my fft does not match my convolution ?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论