why my fft does not match my convolution ?
Show older comments
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])
Accepted Answer
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)
cconv(X,Y,5)/5
Is this what you're trying to show?
8 Comments
Rabih Sokhen
on 4 Mar 2022
Rabih Sokhen
on 8 Mar 2022
Paul
on 8 Mar 2022
Can you verify again that circular convolution is really what you want to show? I'm only asking because the original Question also mentions a filter, and the output of a filter is the linear convolution of its impulse response and input. Can you post the problem as it was given mathematically?
Rabih Sokhen
on 8 Mar 2022
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.
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)
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)
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)
X(w) = fourier(x(t),t,w)
G(w) = H(w)*X(w);
ifourier(G(w),w,t)
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
on 9 Mar 2022
More Answers (0)
Categories
Find more on Transforms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
