CIS 558 / Linguistics 525 -- Homework 2 -- Spring 1999

CIS 558 / Linguistics 525
Computer Analysis and Modeling of Biological Signals and Systems

Homework 2

The MATLAB fft function implements the Discrete Fourier Transform given by the equation

$\begin{displaymath} X(k+1)~~ =~~ \sum_{n=0}^{N-1}~x(n+1)~e^{-i(2\pi/N)kn} \end{displaymath}$

where N = length(x).

In MATLAB terms, this means that for a real (or complex) input vector x of length N, fft returns a complex vector X of length N, where

 
N = length(x);
X = 1:N; % just to create a vector of the correct size
t = 2*pi*(0:(N-1))/N;
for k = 0:(N-1),
  X(k+1) = x*exp(-i*k*t).';
end

(assuming here that x is a row vector, i.e. a matrix of dimensions 1xN).

The MATLAB implementation is quite efficient, using a variety of algorithms depending on the properties of the input.

The use of k+1 on the left in the assignment above reflects the fact that MATLAB vectors begin with 1 rather than 0. We could equivalently have written

for k = 1:N,
 X(k) = x*exp(-i*(k-1)*t)';
end

This problem of whether things start from 0 or from 1 can be very confusing, and often leads to errors. In any case, for data drawn from a time function sampled at SFREQ hertz, the points of the MATLAB fft for indices m from 1 to N/2 represent frequencies

f = (m-1)*SFREQ/N

Thus if we transform 8 samples of a signal sampled 1000 times a second, the first 4 fft points represent frequencies 0*1000/8 = 0, 1*1000/8 = 125, 2*1000/8 = 250, 3*1000/8 = 375.

The fft value at frequency 0 (MATLAB index 1), being the inner product of x with exp(0*t), which is a vectors of 1's, is just the sum of the values of x. This is often referred to as ``DC,'' for ``direct current,'' reflecting the EE heritage of the terminology. If the input x is real, then this sum will obviously be real as well.

The fft value at frequency
SFREQ/2 (MATLAB index (N/2)+1) is the inner product of x with
exp(-i*(N/2)*t).
Graph the real and imaginary parts of this expression in MATLAB, given that t is as above the vector (0:(N-1))*2*pi/N
and you will see that its real part consists of alternating 1 and -1 values, as befits a cosine wave sampled at integer multiples of pi, and its imaginary part is always 0, as befits a sine wave sampled in the same way. Thus this component (like the one at 0 frequency) will always be a real number (for a real-valued input vector). Another way to think about this case is that for a sampled sinusoid with only two samples per period, we cannot distinguish between amplitude and phase, as a little thought will convince you.

What about the rest of the MATLAB fft components from N/2+2 to N? The equation given above determines their values, but we may observe that (for real input x) they are symmetrical, in the sense that for m from 0 to N/2-1,
X(N-m) = conj(X(m+2))
Thus for a length-8 fft,
X(8) = conj(X(2)) ,
X(7) = conj(X(3))
and
X(6) = conj(X(4)).
The values from N/2+1 to N can be treated as negative frequencies, and a simple way to do this, putting frequency 0 in the center of the array, is the MATLAB function fftshift(), which simply swaps the left and right halves of an even-length vector. Treating the component at N/2+1 as a negative frequency is fine, since only the real part is non-zero, and cosine is an even function, that is, one such that f(-t) = f(t).

The MATLAB inverse fft function, ifft(), implements the transform given by the equation

$\begin{displaymath} x(n+1)~~ =~~ (1/N) \sum_{k=0}^{N-1}~X(k+1)~e^{i(2\pi/N)kn} \end{displaymath}$

where N = length(x).

Note the two differences (other than the change of direction): the scaling factor of 1/N, and the sign change in the exponent.

Sine and cosine version of the DFT.
1. Create a set of orthogonal basis vectors for a 16-dimensional space, using sampled sine and cosine waves. Explain each step.
2. Pick three integer frequencies between 1 and 7 inclusive; pick three amplitudes and phases as you please; then create a 16-element vector f by adding up three unit-amplitude sinusoids with the chosen amplitudes and phases.
3. Transform f into a new 16-element vector F whose elements are the inner products of f with the 16 (sinusoidal) orthogonal basis vectors. Do this by setting up an appropriate 16x16 transform matrix.
4. From the transformed vector F, create an amplitude spectrum (that is, a vector showing the amplitude at each of the available frequencies), and a phase spectrum (that is, a vector showing the phase at each of the available frequencies). (Hint: given sine weight sw and cosine weight cw -- and assuming our sinusoids are cosines expressed as amplitude*cos(t-phase) -- the amplitude is sw*sqrt(1+((cw/sw)^2)) and the phase is atan2(cw,sw) ). What is the relationship of these spectra to the recipe that you used to create f in the first place?
5. How would you translate back from F to a new version of f, fnew? Do it. Did it work?
6. Use Matlab's functions fft(), abs(), and angle() to create the amplitude and phase spectrum for f ( amplitude = fftshift(abs(fft(f))), phase = fftshift(angle(fft(f))) ). How does this compare with your version? If there are differences, where do they come from?
Implementations of sft(): the Slow Fourier Transform.
In each case below, verify that your function works correctly by comparing the transform of a random length-8 input vector (from rand(1,8)) with the results obtained from the built-in MATLAB fft().
1. Write a MATLAB function sft() that directly implements the definition of the discrete Fourier transform given above.
2. Write a second version that first sets up a transform matrix (with rows corresponding to the various values of k) and then multiplies this matrix by the input to perform the transform. Then either:
  - Calculate the determinant of the transform matrix (using MATLAB det()). What does the result mean? How would you have to change the transform definition to make the result equal to 1? OR
  - Determine the effect that the matrix has on the (geometric) length of input vectors (hint: try some random inputs and compare the length of the corresponding outputs). Why do the lengths change? How could you change the matrix so that they don't?
3. Write a MATLAB function isft() that directly implements an inverse discrete Fourier transform. Verify that for a random vector x, isft(sft(x)) == x.
Magnitude and phase.
1. Write a MATLAB function to convert fft output to a magnitude and phase form. Do not use abs() and angle(), but write the relationships out in full, to be sure that you understand them; then you can use abs() and angle() to check your result. Test your function by analyzing a signal that is a length-16 vector created by summing three sinusoids with frequencies 1, 3 and 5 (assuming a sampling frequency of 16), with random amplitudes and phases. In MATLAB, this will be something like:
```
T = 2*pi*(0:15)/16;
a = rand(1,3);
p = 2*pi*rand(1,3) - pi;
x = a(1)*cos(T+p(1)) + a(2)*cos(3*T+p(2)) + a(3)*cos(5*T+p(3));
```
  Now verify that you can recover the frequencies, amplitudes and phases correctly from fft(x). Be aware that fft components with magnitudes that are effectively zero, but not precisely zero due to round-off error in the calculation, will have some phase or another with a value that is not very interesting.
2. Write a MATLAB function to display the magnitude and phase spectrum of a real signal, showing only the non-redundant values. Make sampling frequency a parameter, and number the frequency axis correctly.
Symmetry properties:
1. Construct a real, random, symmetric signal (vector x) of length 64. In MATLAB terms, ``symmetric'' here means that in the vector xx, xx[33+n] = xx[33-n] for n from 1 to 31. The point of this construction is related to the concept that we are looking at one period of a periodic signal, where we think of the first sample as being at time 0; we then ask whether the resulting (infinite-time) signal is even (symmetric, x(t) = x(-t)) or odd (antisymmetric, x(t) = -x(-t)). Look at your signal's DFT (using Matlab's FFT function). Describe what you see.
2. Now do the same for a random antisymmetric real signal, and also for a random real signal that is neither symmetric nor antisymmetric.
What do you think would happen if you did these three experiments on an imaginary signal? Try it.
Shift properties:
1. Construct a real 16-element vector f with Fourier magnitude $\vert F(k)\vert = \vert\frac{N}{2} - k\vert$ , and random phase. Do this by first building the magnitude and phase vectors in the frequency domain (set the 16 magnitudes according to the equation above; set the phases by a call to rand(), except obeying the symmetry required for a real time-domain vector, then convert to a vector of complex coefficients), and then inverse-fourier transform (using matlab's ifft function: f = ifft(F) ).
2. Check that the whole construction worked by Fourier transforming the resulting function and viewing the Fourier magnitude (by applying Matlab's abs function to the fft of f ).
3. Now, circular-shift the signal f by $\frac{N}{2}$ samples, and again compute the Fourier magnitude and phase (in Matlab, you could implement a circular shift of N/2 samples as f1 = [f(9:16); f(1:8)], assuming that f is a column vector). How do they compare to the original? In particular, carefully compare the phase to that of the original (using Matlab's angle() function).
4. What happens to the Fourier magnitude with circular shifts of other sizes?
5. What happens to the Fourier magnitude and phase if you reverse f? (in Matlab, frev = f(16:-1:1) )
6. Create a new signal g[n] = f[n] * (-1)ⁿ. That is, multiply every other sample of f by -1. Again, compute the Fourier magnitude and phase. How do these compare to those of the original? Why?
Two-dimensional DFT
In this question, use the MATLAB function fft2() to calculate the 2-dimensional dft of various matrices. This function computes the 1-D dft of each column of the input, then computes the 1-D dft of each row of the result.
1. Make a 64x64 image containing an impulse (sample of value 1) at location (1,1). What is its DFT magnitude? Use the Matlab function fft2 to compute this.
2. Compute the DFT magnitude of a 64x64 image $f[n,m] = \cos( 2\pi (8n + 16m)/64)$ . Why does it look like it does? What if the image was constructed using the $\sin$ function?
3. Construct an image f[n,m] = e^{-[(n-32)²+(m-32)²]/32}, and compute the magnitude of its DFT.
4. What is the the DFT of the pointwise product of the previous two images? Why?
5. Construct a 256x256 image f(n,m) = [(n-128)² + (m-128)²]^-0.7. Note that for f(128,128), by the equation given, the value is 0^-0.7 = 1/0. Rather than complain, MATLAB will cheerfully set this to Inf ("infinity"), which will have a bad influence on the rest of your efforts. So set f(128,128) to some plausible and finite value, say 1.
  Circular-shift the image f by 128 pixels rightward and upward. Call this image M. Create a random image of the same size containing values between 0 and $2\pi$ . Call this P. Consider these as the magnitude and phase of some image. Compute the image. What does it look like? What is the magnitude of its DFT?
6. Now combine M with the phase of the DFT of the einstein image. Invert the transform. What does this image look like?

Mark Liberman
2/12/2003