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) + ...]

(1)
**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)

where

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.

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