why my fft does not match my convolution ?

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])

Answer 1

[Paul]

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?

[Rabih Sokhen]

hello Paul

yes this is exactly what i want to do

i tried to apply your example in my code but i did not get the same amplitude,

this is my modification:

code:

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')

can you help me solve this plz

thank you in advance

[Paul]

Hard to say without seeing the complete code. Having said that ...

Why does g_f1 apply ifftshift() to g_t? Why does g_f1 include a muliply by dt?

[Rabih Sokhen]

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]

@Rabih Sokhen

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]

I'm not very familiar with circular or linear convolution, i am trying to verify that g(w)=fft( h_t . x_t) = conv( h_f,x_f) by considering a impulse response h_t and by filtering h_t with x_t.

[Paul]

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]

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]

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