Thursday, September 15, 2011

Fourier Series

Understanding the Fourier Transform is not trivial. Here is a simpler approach (only for periodic signals) from the same author termed the Fourier Series.
Joseph Fourier [Wikipedia] demonstrated that any periodic function could be decomposed in a series of sine and cosine functions.

1- The square wave/signal [Wikipedia]

Imagine an image of vertical strips alternating black and white colors as shown in Fig.1 . The script to draw this image is here [see post]. Then, draw a horizontal line (Tool #5) and compute a Analyze > Plot Profile (shortcut: Ctrl +K), you get what is called a square wave (also called step function).
Fig. 1: Vertical white and black strips

According to Fourier, this kind of curve is the sum of sine and cosine functions. The formula - in this case - is the following:

f(t)=4*[ 1/1*sin(1*t) + 1/3*sin(3*t) + 1/5*sin(5*t) + ...]
Note: This is not magic, the coefficients of the Fourier Series can be determined but this is beyond the scope of this post. Here, there is no cosine because this is an odd function.

2- How can we check this formula in ImageJ?

The simplest way is to create an image where the X-axis corresponds to t (or time) and the Y-axis to the frequencies.

Thus, the algorithm can be divided in two parts:
  • first, we draw each sine function in a separate image row
  • second, we sum all these rows to get the final result. This has the advantage to sum the whole sine functions or just a subset (by defining a region of interest).
Refresher: the formula of a sine function is
f(t) = A*sin(B*t + C)
A is the amplitude
2*π/|B| = T is the period and f (the frequency) is equal to 1/T
C is the phase shift.

2-1- Drawing a sine function in ImageJ
To draw such a curve, we have to convert the x-coordinate of the pixel in a value comprised between 0 and 2π, this is simply done by dividing x by the image width and multiplying the result by 2π.
Now, open the recorder Plugins > Macro > Record...
- Create a 32-bit image entitled sineFunction of 256x256 pixels with a black background.
- Go to Process > Math > Macro... and type the following formula:

v=sin(x / w * 2 * PI)
click on the OK button and the following image is drawn.
- Finally, for computing a plot profile, draw a horizontal line (Tool #5 and Shift + click and drag) and go to Analyze > Plot Profile (or Ctrl + K).

Fig. 2: Drawing of a sine function and the corresponding plot profile
2-2- Drawing Fourier Series in ImageJ
Now, go back to the Fourier Series, in the Recorder, click on the 'Create' button, clean up your code and keep two lines: (i) the image creation and (ii) the formula.
- Each sine term is function of the frequencies (y-coordinate in our case). We have to slightly modify the formula in Process > Math < Macro...  leading to:

v = sin ( y * x / w * 4 * 2 * PI) / y

Moreover, we need to restrict the drawing to odd frequencies:

if (y%2==1) v = sin ( y * x / w * 4 * 2 * PI) / y

Et voilà. The script becomes

+++ IJ snippet +++
+++ end of IJ snippet +++

Fig.3: Fourier Series of a square wave

Each row of the resulting image contains a sine term. The frequency 0 is at the top.

Where is the square wave ?
According to the formula (1), we need to sum all these rows to get the square wave. To do that, choose the Rectangle Tool (Tool #1), select the whole image (Ctrl + A), and move the bottom edge of the rectangle in such a way you select the first 10 rows. Then, run Analyze > Plot Profile to display the sum of all the rows contained in the rectangular selection area.
Fig. 4: Plot calculated from the selected area (bottom left window). In ImageJ, a plot calculated from a rectangular area sums all the horizontal lines. The resulting plot approximates our square wave.

Try again with different selection area, to see the evolution of the curve.

Note: For high frequencies, there are some strange oscillating lines. This is an artifact known as Gibbs phenomenon.

2-3 The power spectrum
Now to summarize the result, a diagram (termed the power spectrum) is calculated displaying the amplitude in function of the frequencies.

Fig. 5: Power spectrum of the square wave.

Fig. 6: Power spectrum of the square wave calculated with FFT of ImageJ

3-Animating the curves

Here is a small script to follow the evolution of the sines sum...

+++ IJ snippet +++
+++ end of IJ snippet +++

No comments:

Post a Comment