ECE 401 MP5: Discrete Fourier Transform

In this file, we're going to use the inverse DFT to create images from MRI spectral measurements. Then we'll use the DFT to filter out noise in an image.

0. Data

The data are provided for you in the file data.hdf5. Let's load it, and see what it contains.

1. MRI

MRI (magnetic resonance images) are acquired in the frequency domain. It works roughly like this:

  1. Every atom in your body spins, all the time, emitting weak radio signals. A strong magnetic field will align a few percent of them all in the same direction, so that their radio signals can be measured.
  2. If the external magnetic field slightly different in different locations, then the synchronized spin will be at slightly different frequencies in different locations. By holding this difference for a carefully controlled amount of time, we can set it up so that your atoms are spinning with slightly different phases in different locations. The phase at location $(n_1,n_2)$ might be tweaked to be $\phi=-\omega_1n_1-\omega_2n_2$.
  3. If the density of hydrogen atoms at location $(n_1,n_2)$ is given by $x[n_1,n_2]$, then the total signal collected is
$$X(\omega_1,\omega_2) = \sum_{n_1}\sum_{n_2} x[n_1,n_2] e^{-j\omega_1 n_1}e^{-j\omega_2 n_2}$$
  1. After measuring $X(\omega_1,\omega_2)$ at $N_1$ different values of $\omega_1$, and $N_2$ different values of $\omega_2$, then we compute the inverse DFT in order to calculate the image.

The image you are provided, mri_dft, is a set of complex-valued frequency samples from a real MRI image (a public domain image from Wikipedia, donated by Johns Hopkins University). Let's try taking the inverse DFT, to see if we can see the image.

The aliasing we see, here, is caused by capturing only $N_1\times N_2$ frequency samples, but computing the DFT using a spatial range $M_1\times M_2$ such that $M_1>N_1$ and/or $M_2>N_2$. By doing that, we accidentally discover that the DFT is periodic in space, not just in frequency:

$$x[N+n] = \frac{1}{N}\sum_k X[k]e^{j\frac{2\pi k(n+N)}{N}} = \frac{1}{N}\sum_k X[k]e^{j\frac{2\pi kn}{N}} =x[n]$$

In order to solve this problem, we need to figure out how many of the frequency samples in mri_dft are actually nonzero-valued. To get an idea, let's try plotting the first row:

It looks like maybe the odd-numbered spectral samples are zero, and only the even-numbered spectral samples are valid. Let's try taking the inverse fft of only the even-numbered samples.

That's the right size, but now the whole image has been shifted by about $N_1/2$ rows, and by about $N_2/2$ columns! One way to solve this is by multiplying the DFT by a shift operator, before we inverse DFT:

$$X_{\mbox{corrected}}[k_1,k_2] = \begin{cases} X[k_1,k_2] e^{-j\frac{2\pi k_1d_1}{N_1}}e^{-j\frac{2\pi k_2d_2}{N_2}} & 0\le k_1<\frac{N_1}{2},~0\le k_2<\frac{N_2}{2}\\ X[k_1,k_2] e^{-j\frac{2\pi (k_1-N_1)d_1}{N_1}}e^{-j\frac{2\pi (k_2-N_2)d_2}{N_2}} &\frac{N_1}{2}\le k_1<N_1,~\frac{N_2}{2}\le k_2<N_2 \end{cases} $$

where $d_1$ is the number of rows that we want to shift, and $d_2$ is the number of columns that we want to shift. Notice that the second line is actually equal to the first line, in theory (because $e^{-j\frac{2\pi N_1d_1}{N_1}}=1$), but implementing it in the way shown above will help to avoid artifacts caused by floating-point roundoff error. Notice that the formula above only lists two of the four cases; I think you can figure out what the other two cases should be.

2. Noisy Image

Filtering out noise is easier in the time domain, if the noise is wideband noise (meaning that there is noise in a wide range of frequencies). Sometimes, though, the noise is in a narrow range of frequencies. For example, consider this image (this is the image "Clifton Beach 5, South Arm, Tasmania, Australia" by JJ Harrison, distributed on Wikimedia under Gnu Free Documentation License):