Math 584, Mathematics of Medical Imaging

Matlab Lesson 3: Hilbert transform, shift invariant filters, and sampling issues

 

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.

The Finite Fourier transform

The Finite Fourier transform is implemented in matlab using the "fft" command:
>> X=-10:0.25:10; Y=1./(1+X.^2); % sample the function
>> Y_fft=fft(Y);
% This computes the Finite Fourier transform of the samples Y. If we plot Y_fft then we find that the positive frequencies come first and then the negative frequencies.
>> xaxis=1:length(Y_fft); figure, plot(xaxis,real(Y_fft),xaxis, imag(Y_fft))% With this ordering the Fourier tranform has both  real and imaginary parts


>>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.

  >>figure, plot(xaxis,real(Y_fft),xaxis, imag(Y_fft)) % Now we see that zero frequency is in the center. Also there is no imaginary part. Do you know why?

 

 

The Finite Fourier transform and Hilbert transform

 

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.

 

Derivatives with fft

 

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.

 

Convolution with fft

 

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.

 

 

Exercises

 

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