Please download the new version of ft_demo.m and ft_demo.fig. This has some new options, including the Hilbert transform, convolutions, and derivatives.
>>Y_fft=fftshift(fft(ifftshift(Y))); % Using these shift commands puts the frequencies in the natural order.
The
command fftshift shifts
the entries of a vector so that the first entry becomes the “central”
entry. The inverse operation is ifftshift. For
examply, try >> fftshift([1
2 3 4 5]);
This is useful because the finite fourier
transform (fft)
of a vector treats the first entry as corresponding to the 0th
frequency, but when considering continuous functions we tend to think
of the 0th
frequency as being in the middle.
The Hilbert transform can be implemented in Matlab in four steps: (1) sample the function
(2) apply the
Fourier transform, (3) multiply the negative frequencies by -1 and the
zero frequency
by 0, (4) apply the inverse Fourier transform.
For example, to compute the Hilbert transform of exp(-x^2):
>> X=-10:0.01:10; Y=exp(-X.^2);
% sample the function : Y
consists of samples of exp(-x^2) at the sample points
%in X.
>> Y_fft=fftshift(fft(ifftshift(Y))); % apply the
Fourier transform using the shift operations as explained above.
>> Y2_fft=Y_fft; Y2_fft(find(X<0))=Y2_fft(find(X<0))*(-1);
Y2_fft(find(X==0))=0; % multiply the negative freqs
by -1 and the zero freq by 0
>> Y2=fftshift(ifft(ifftshift(Y2_fft)));
% apply the inverse Fourier
transform
>> figure; plot(X, real(Y2),X,imag(Y2)); % This will plot the Hilbert transform of the original function
The
command find can be used to find a subset of
indices of a vector. Above, we used the
command to multiply all of the negative frequencies by -1 and the zero
frequency by 0. We are taking advantage of the fact that
the indices of the negative entries of X are the same as the indices of
the negative frequencies Y_fft.
You can use the ft_demo.m
program to experiment with the Fourier transform of various functions.
The finite Fourier
transform is most
efficiently implemented for vectors of size 2^n… see the exercise below.
The fft
can be used to approximate the derivative of a sampled function. Use the following steps:
(1)
Sample the function on a uniform grid X
(2)
Apply the Fourier transform
(3)
Multiply the Fourier transform by i*x
(here x is the
grid in the Fourier domain)
(4)
Apply the inverse Fourier transform.
For
example, to compute the
numerical derivative of 1/(1+x^8):
>>
X=-10:0.01:10; Y=1./(1+X.^8); figure, plot(X,Y);
% sample the
function
>> Y_fft=fftshift(fft(ifftshift(Y))); k=X/0.01/20; %
Apply the
Fourier transform. k is the vector of frequencies.
>> Ydiff_fft=Y_fft.*(i*k); % multiply the Fourier coefficients by i*k
>> Ydiff=fftshift(ifft(ifftshift(Ydiff_fft))); %
apply the inverse Fourier
transform
>> Ydf =
-(8*X)./((1+X.^8).^2);
% We compute the derivative
algebraically.
>> figure;
plot(X,real(Ydiff),'b',X,imag(Ydiff),'r',
X,Ydf,'g');
% plot the result. We
show the actual derivative (X,Ydf) as the green curve.
We see that
the scaling is not correct in the Fourier computation of the
derivative. We
correct the scale by comparing the maxima of Ydf and Ydiff
>> scale=max(Ydf)/max(abs(Ydiff)); %This give the correct scale factor.
Often in applications it is important to use an algorithm to compute a
simple example, in order to get the scaling correct.
>> figure;
plot(X,real(scale*Ydiff),'b',X,imag(scale*Ydiff),'r',
X,Ydf,'g');
Try
experimenting with
derivatives of functions. You can use ft_demo.
Up to a
constant factor, the
following procedure can be used to convolve two functions, f and g.
(1)
Sample f and g on a uniform grid
(2)
Apply the fourier
transform
to obtain f_fft and g_fft
(3)
Pointwise multiply h_fft = f_fft
.* g_fft
(4)
Apply the inverse Fourier transform to the convolution h.
Note that
this is procedure will
convolve f and g as periodic functions (or functions on the circle). This convolution will only coincide with a
convolution on the real line, over a certain interval, if f and g
are supported on a sufficiently
small
interval within the sampling grid. Thus
to obtain an accurate convolution of functions defined on the real
line, zero
padding (or sampling on a larger interval) may be necessary (see the example below).
Try
convolving various
functions. You can use ft_demo.
Exercise 1. What are the symmetries of the Hilbert transform and the Fourier transform? Use the ft_demo program to experiment with the Hilbert/Fourier transforms. What can you say about the Hilbert/Fourier transform of a real, pure imaginary, even, odd function? Verify with examples that the Hilbert transform is shift invariant. Is the Fourier transform shift invariant?
Exercise 2. The fast fourier transform (fft) is most efficient for vectors of length 2^n. Use the tic and toc commands to time the fft for vectors of various lengths. For example, try
>>
tic; fft(ones(1,2^17));
toc
How much time elapsed? Now try a vector of length 2^17-1 (note that 2^17-1 is a prime number). Is there a significant difference in the elapsed time? What about 2^17+3? What is the factorization of 2^17+3 (you can use the factor procedure).
Exercise 3. When numerically computing the derivative of a function (using the above fft procedure), is it more important to use a small sample step or to sample the function on a large interval? Try some examples with ft_demo. For example, compute the derivative of 1/(1+x^8) using a large sample step of 0.3, or a small interval of [-1.5, 1.5]. What is your conclusion. Explain.
Exercise 4. Use ft_demo to convolve 1/(1+x^8) with f(x)=(-0.1<=x).*(x<0.1)/(2*0.1). What is this f(x)? What is the effect of convolving with f(x)?