why my fft does not match my convolution ?

15 visualizaciones (últimos 30 días)
Rabih Sokhen
Rabih Sokhen el 3 de Mzo. de 2022
Editada: Paul el 9 de Mzo. de 2022
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])

Respuesta aceptada

Paul
Paul el 4 de Mzo. de 2022
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 comentarios
Paul
Paul el 9 de Mzo. de 2022
Editada: Paul el 9 de Mzo. de 2022
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
Rabih Sokhen el 9 de Mzo. de 2022
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

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Fourier Analysis and Filtering en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by