Fourier Transform

Home Fourier Series Fixed Point Arithmatic

In the previous section of Fourier Series a data set was described as a series of sin and cos functions and the whole exercise revolved around finding the coefficients $$(a,b$$ and $$ c)$$ of the series assuming the data set had a well defined period.

Here is a Fourier Series for continuos time function

$$ \\ \\ f(t) = c+\sum\limits_{n=1}^{\infty} a_n sin(\omega_n t) + \sum\limits_{n=1}^{\infty} b_n cos(\omega_n t)\right) \hdots \hdots \hdots(2.1)\\ \\ \\ a_n=1/\pi \int\limits_{t=0}^{2\pi} f(t) sin(\omega_n t) dt \\ b_n=1/\pi \int\limits_{t=0}^{2\pi} f(t) cos(\omega_n t) dt \\ c={ 1\over {2\pi}}{\int\limits_{t=0}^{2\pi} f(t) dt} $$

One drawback in the Fourier Series is that it assumes the complex wave has a predefined period. The formula requires data for one complete cycle and expects the behavior to be repeated. But who knows what the wave looked like before and what will it look after ward or is it really a periodic wave per se? After all, nature is too complex to have a repeated pattern. And what are the consequences of taking few samples and basing a judgment on the entire function?

To answer these questions we need a formula that does not depend upon the function being periodic and this is what Fourier Transform is all about, it is based upon the following principle, a sample from a lot tells you something about a lot and more samples from the lot tells you more about the lot, so let’s modify the Fourier Series and remove the periodic dependency and see the consequences of it on the original function.

In Fourier Transform we may obtain a frequency spectrum where coefficient values are not exact but as long as the relative strength is preserved and there is a way to bring back the origional function the objective is achieved.

With this goal in mind let's proceed with the Transform, but first, some words about complex numbers.

Complex Number Arithmatic

Think of complex number as a vector on a coordinate system represented by two numbers, x and y. It is customary to show the coordinate system y-axis as the imaginary axis and x-axis as the real axis. We need to add the imaginary operator $$\sqrt{-1}$$ with the y in order for the number to be used in algebraic expressions and follow rules of vector arithmatic.

$$ (x1+\sqrt{-1}y1) \times (x2+\sqrt{-1}y2) = (x1x2 -y1y2)+\sqrt{-1}(x2y1 + x1y2) \\ $$

$$ (x1+\sqrt{-1}y1) + (x2+\sqrt{-1}y2) = (x1 + x2)+\sqrt{-1}(y1 + y2) \\ $$

$$ (x1+\sqrt{-1}y1) - (x2+\sqrt{-1}y2) = (x1 - x2)+\sqrt{-1}(y1 - y2) \\ $$

$$ {(x1+\sqrt{-1}y1) \over (x2+\sqrt{-1}y2)} = {(x1+\sqrt{-1}y1) \over (x2+\sqrt{-1}y2)} \times {(x2+\sqrt{-1}y2) \over (x2+\sqrt{-1}y2)} = {{(x1x2 -y1y2)+\sqrt{-1}(x2y1 + x1y2)} \over {x2^2 + y2^2}} \\ $$

In the first step, we will consolidate the three coefficients $$a,b$$ and $$c$$ into a single coefficient using Euler's identity.

$$ \\ \\ e^{j\theta}=cos(\theta) + jsin(\theta) \\ \\ e^{-j\theta}=cos(\theta) - jsin(\theta) \\ \\ cos(\theta) = {e^{j\theta}+e^{-j\theta} \over 2 } \\ \\ sin(\theta) = {e^{j\theta}-e^{-j\theta} \over 2j } \\ \\ $$

Replacing the sin and cosine functions of Equation 2.1 with its equivalent exponent notation

(replacing $$1/j$$ with $$-j$$), we get,

$$ \\ \\ f(t) = a_0 +\sum\limits_{n=1}^{\infty} {-ja_n\over 2} \times (e^{jn\omega_0 t} - e^{-jn\omega_0 t}) + \sum\limits_{n=1}^{\infty} {b_n\over 2} \times (e^{jn\omega_0 t} + e^{-jn\omega_0 t})\\ \\ f(t) = a_0 +\sum\limits_{n=1}^{\infty} {-ja_n\over 2} \times (e^{jn\omega_0 t}) + \sum\limits_{n=1}^{\infty} {ja_n\over 2} \times (e^{-jn\omega_0 t}) + \sum\limits_{n=1}^{\infty} {b_n\over 2} \times (e^{jn\omega_0 t}) + \sum\limits_{n=1}^{\infty} {b_n\over 2} \times (e^{-jn\omega_0 t})\\ \\ f(t) = a_0 +\sum\limits_{n=1}^{\infty} {-ja_n + b_n\over 2} \times (e^{jn\omega_0 t}) + \sum\limits_{n=1}^{\infty} {ja_n + b_n\over 2} \times (e^{-jn\omega_0 t})\\ \\ f(t) = a_0 +\sum\limits_{n=1}^{\infty} {(A_n)} \times (e^{jn\omega_0 t}) + \sum\limits_{n=1}^{\infty} {(B_n)} \times (e^{-jn\omega_0 t})\\ \\ A_n = {-ja_n + b_n\over 2}, B_n = {ja_n + b_n\over 2} \\ $$

We can merge the coefficients $$a_n$$ and $$a_0$$ into $$a_n$$ since for n=0 the two exponents are same and simplify the Fourier series

$$ \\ \\ f(t) = \sum\limits_{n=0}^{\infty} {(A_n)} e^{jn\omega_0 t} + \sum\limits_{n=1}^{\infty} {(B_n)} e^{-jn\omega_0 t}\\ \\ $$

We could modify the second summation so that both summation combined would cover the entire spectrum of $$-\infty$$ to $$+\infty$$, but for that we need to redefine the negative values of the exponent n as follows:

$$ \\ \\ f(t) = \sum\limits_{n=0}^{\infty} {(A_n)} e^{jn\omega_0 t} + \sum\limits_{n=-\infty}^{-1} {(B_n)} e^{jn\omega_0 t}\\ \\ $$

Combining the two summation for the entire range of the values for $$n=-\infty$$ to $$+\infty$$ we can describe the Fourier series equivalent of the original function into the complex exponent form as (replacing $$A_n and B_n$$ with $$X_n$$)

$$ \\ \\ f(t) = \sum\limits_{n=-\infty}^{+\infty} {X(n)} e^{jn\omega_0 t} \cdots Eq(2.2)\\ \\ $$

And the coefficient itself can be modified using Euler's identity

$$ \\ \\ X(n) = { 1\over 2\pi} \int\limits_{0}^{2\pi} {f(t)} e^{-jn\omega_0 t}dt \cdots X(n) = { 1\over T} \int\limits_{0}^{T} {f(t)} e^{-jn\omega_0 t}dt \cdots X(n) = { 1\over T} \int\limits_{-T\over 2}^{T \over 2} {f(t)} e^{-jn\omega_0 t}dt \cdots Eq(2.3)\\ \\ $$

Note: Shifting the integral halfway will merely introduce mirroing of the frequencies on the negative side

You may say that now there is a certain amount of elegance in the representation of the Fourier series in Equation 2.2 and 2.3, that we have taken a long series and reduced it to a compact looking function. But we certainly could not use this formula for actual calculations. We need the coefficients $$a_n$$ and $$b_n$$ and for that we still have to go through the same computational process that we did with the long series. Hidden in the coefficient of $$X(n)$$ are $$a_0$$, $$a_n$$ and $$b_n$$ and the $$e^{jn\omega_0 t}$$ hides the $$cos(\omega t)$$ and $$sin (\omega t)$$. Still, the Equation 2.2 expresses the relationship between an input and output.

Removing the periodic dependency

To get complete spectrum of the frequencies in a function we need to observe the function for eternity. Does it mean the job can never be completed? Having a period that approaches infinity has other problems too. It would make the fundamental harmonics near 0 as shown in the Equation 3.12.

$$ \\ \\ f_0 = { 1\over T_{\rightarrow \infty}} \rightarrow { 1 \over \infty} \rightarrow 0 \\ $$

The argument presented in the previous section leads us to one conclusion that all data set should be considered a partial set from an infinite set that spans an infinite period. If a wave takes infinite time to advance through one period then its frequency approaches 0. What happens to the spectrum of the frequency, as the wave passes through only the smallest fraction of its period in one second? The spectrum shows only the integral multiple of the fundamental frequency, so as the frequency approaches 0, the interval between frequencies on the spectrum grows narrower and narrower. Eventually as the period approaches infinity the spectrum becomes continuum and all frequencies are known and there is no gap between the frequencies. Since, now every sample is considered a contributing factor towards the frequency component and the number of samples is related to the number of frequencies in our spectrum, it is no longer a time domain any more. The domain is our frequency centered around $$f_0 = 0$$.

Having realized that frequency a continuum function, the fundamental and the harmonics are simply a continuous frequency variable with no gaps among them. Let’s rewrite the Equation 2.2 and use the term $$g(t)$$ instead of $$f(t)$$, as $$f$$ is being used for frequency.

$$ \\ \\ g(t) = \sum \limits_{n=-{\infty}}^{\infty}{X(\omega_n)} \times e^{jn2\pi f_0t} \cdots (2.4)\\ \\ $$

$$ \\ \\ X(\omega_n) = {1 \over T}\times\int\limits_{t=-{T\over2}}^{T\over 2}{g(t)}e^{-jn2\pi f_0t}dt \cdots (2.5)\\ \\ $$

Replacing the value of $$X(\omega_n)$$ in Equation 2.4 with the Equation 2.5.

$$ \\ \\ g(t) = \sum \limits_{n=-{\infty}}^{\infty} { \begin {pmatrix} {1 \over T} \times \int\limits_{t=-{T\over2}}^{T\over 2} g(t) e^{-jn2\pi f_0t}\times dt \end {pmatrix} } e^{jn2\pi f_0t} \cdots (2.6)\\ \\ $$

The sticky point in the Equation 2.6 is the multiplying factor $$1 \over T$$, (which is now very close to 0,) if we eliminate the multiply we may not get the actual coefficient of the sinusoid, but still the term will be unique for the sinusoid and we will leave it that way. We can also replace the quantity $$T\over 2$$with $$\infty$$ and $$-T\over{2}$$ with $$-\infty$$, since $$T$$ is infinite so $$T\over 2$$ is also infinite. The term in the bracket of Equation 2.6 (without the multiply $$1 \over T$$) is our Fourier Transform. It has the desired property that it does not depend upon the periodicity.

$$ \\ \\ X(\omega_n) = \int\limits_{t=-{\infty}}^{\infty} g(t) e^{-jn2\pi f_0t}\times dt \cdots (2.7) \\ $$

And to bring back our original function we have Inverse Fourier Transform by substituting the value back in Equation 2.6.

$$ \\ \\ g(t) = \sum \limits_{n=-{\infty}}^{\infty}{X(\omega_n)} \times e^{nj2\pi f_0t} \cdots (2.8)\\ \\ $$

The Fourier Transform simply tells us the relative magnitude of the component frequencies, and this is good enough for spectrum analysis.

Discrete Time Fourier Transform

The analog world may be continuous in time but when a computer is being used to process a signal, there is always a finite amount of time between acquisition and subsequent processing, no matter how fast the computer is. The method is essentially a sampling of event taken at discrete time. The mathematics has to be modified slightly to take into account that only at specific time the discrete time coincides with the real time.

The discrete time equivalent of Fourier Series

$$ \\ \\ f(k\Delta T) = c + \sum\limits_{n=1}^N a_n sin(n2\pi k\Delta T/T) + \sum\limits_{n=1}^N b_n cos(n2\pi k\Delta T/T)\right) \hdots \hdots \hdots(2.9) \\ \\ a_n = {2\over T} { \sum\limits_{k=0}^{K-1} {f(k\Delta T) \times sin({n2\pi \over T} k\Delta T ) \times \Delta T}} \\ \\ b_n = {2\over T} {\sum\limits_{k=0}^{K-1} {f(k\Delta T) \times cos({n2\pi \over T} k\Delta T ) \times \Delta T}} \\ \\ c = {1\over {T}}{\sum\limits_{k=0}^{K-1} f(k\Delta T) \times \Delta T} \\ \\ $$

If $$\Delta T$$ is our sampling period then the $$kth$$ sample time that coincide with the real time is $$k\Delta T$$. Each discrete sample in time of a function $$g(t)$$ is being placed in the rspective $$kth$$ position of the function and is being represented as $$g(k\Delta T)$$. If $$\omega$$ is the angular velocity in radians per second then $$\omega \Delta T$$ is the advancement of the wave in radians per sample. If $$N$$ is the number of samples per second then the nth frequency component is $$f_n=n/(\Delta t N)$$, where $$n =0,1,2,3,4…N$$. Modifying the Fourier Transform of Equation 2.7 for discrete time events we have discrete time Fourier Transform.

$$ \\ \\ G{\begin {pmatrix} n \end {pmatrix}} = \sum\limits_{k=0}^{N-1}{\Delta T} \times g(k\Delta T) \times e^{-j2\pi \Delta T k {n\over \Delta TN}} \cdots (2.10) \\ $$

Since $$\Delta T$$ is our choice, we could normalize the Equation 2.10 by choosing $$\Delta T = 1$$ and simplify the Fourier coefficient formula as,

Discrete Time Fourier Transform

$$ \\ \\ G{\begin {pmatrix} n \end {pmatrix}} = \sum\limits_{k=0}^{N-1} \times g(k) \times e^{-j2\pi k {n\over N}} \cdots (2.11) \\ $$

And to bring back our original function we have Inverse Fourier Transform by substituting the value back in Equation 2.8.

Discrete Time Inverse Fourier Transform

$$ \\ \\ g(k) = \sum \limits_{n=0}^{N - 1}{G(n)} \times e^{j2\pi k {n\over N}} \cdots (2.12)\\ \\ $$

Table 2.1

$$k$$	$$k\Delta T$$	$$2\pi fk\Delta T$$	$$f(k\Delta T)$$
0	0	0	7.5
1	0.0625	0.4	7.32
2	0.125	0.8	5.27
3	0.1875	1.2	4.04
4	0.25	1.6	4.7
5	0.3125	2	5.5
6	0.375	2.4	4.27
7	0.4375	2.8	1.38
8	0.5	3.2	-0.5
9	0.5625	3.6	0.39
10	0.625	4	2.73
11	0.6875	4.4	3.66
12	0.75	4.8	2.3
13	0.8125	5.2	0.79
14	0.875	5.6	1.73
15	0.9375	6	4.91

Notice $$e^{j2\pi k{n \over N}}$$ and $$e^{-j2\pi k{n\over N}}$$ are complex numbers that must be evaluated in terms of corresponding sin and cos values:

$$ \\ e^{j2\pi k{n\over N}} = cos (n2\pi k{n\over N}) + jsin (2\pi k{n\over N}) \\ $$

$$ \\ e^{-j2\pi k{n\over N}} = cos (n2\pi k{n\over N}) - jsin (2\pi k{n\over N}) \\ $$

The following table is a tabulation of Fourier Coefficients of sample data of Table 2.1 using Eq(2.11)

$$ k $$	$$f(k)$$	$$ cos_0 $$	$$ sin_0 $$	$$ cos_1 $$	$$ sin_1 $$	$$ cos_2 $$	$$ sin_2 $$	$$ cos_3 $$	$$ sin_3 $$	$$ cos_4 $$	$$ sin_4 $$	$$ cos_5 $$	$$ sin_5 $$	$$ cos_6 $$	$$ sin_6 $$	$$ cos_7 $$	$$ sin_7$$
0	7.5	7.5	0	7.5	0	7.5	0	7.5	0	7.5	0	7.5	0	7.5	0	7.5	0.00
1	7.32	7.32	0	6.76	2.8	5.18	5.18	2.8	6.76	0	7.32	-2.8	6.76	-5.18	5.18	-6.76	2.80
2	5.27	5.27	0	3.73	3.73	0	5.27	-3.73	3.73	-5.27	0	-3.73	-3.73	0	-5.27	3.73	-3.73
3	4.04	4.04	0	1.55	3.73	-2.86	2.86	-3.73	-1.55	0	-4.04	3.73	-1.55	2.86	2.86	-1.55	3.73
4	4.7	4.7	0	0	4.7	-4.7	0	0	-4.7	4.7	0	0	4.7	-4.7	0	0	-4.70
5	5.5	5.5	0	-2.1	5.08	-3.89	-3.89	5.08	-2.1	0	5.5	-5.08	-2.1	3.89	-3.89	2.1	5.08
6	4.27	4.27	0	-3.02	3.02	0	-4.27	3.02	3.02	-4.27	0	3.02	-3.02	0	4.27	-3.02	-3.02
7	1.38	1.38	0	-1.27	0.53	0.98	-0.98	-0.53	1.27	0	-1.38	0.53	1.27	-0.98	-0.98	1.27	0.53
8	-0.5	-0.5	0	0.5	0	-0.5	0	0.5	0	-0.5	0	0.5	0	-0.5	0	0.5	0.00
9	0.39	0.39	0	-0.36	-0.15	0.28	0.28	-0.15	-0.36	0	0.39	0.15	-0.36	-0.28	0.28	0.36	-0.15
10	2.73	2.73	0	-1.93	-1.93	0	2.73	1.93	-1.93	-2.73	0	1.93	1.93	0	-2.73	-1.93	1.93
11	3.66	3.66	0	-1.4	-3.38	-2.59	2.59	3.38	1.4	0	-3.66	-3.38	1.4	2.59	2.59	1.4	-3.38
12	2.3	2.3	0	0	-2.3	-2.3	0	0	2.3	2.3	0	0	-2.3	-2.3	0	0	2.30
13	0.79	0.79	0	0.3	-0.73	-0.56	-0.56	-0.73	0.3	0	0.79	0.73	0.3	0.56	-0.56	-0.3	-0.73
14	1.73	1.73	0	1.22	-1.22	0	-1.73	-1.22	-1.22	-1.73	0	-1.22	1.22	0	1.73	1.22	1.22
15	4.91	4.91	0	4.54	-1.88	3.47	-3.47	1.88	-4.54	0	-4.91	-1.88	-4.54	-3.47	-3.47	-4.54	-1.88
	$$\sum$$	55.99	0	16.01	12	0.01	4	16	2.39	0	0.01	0	0	-0.01	0	-0.01	0.01

Fourier Coeficients Magnitude

$$ \\ G(0) = \sqrt {(55.99)^2 + (0.00)^2} = 55.99 \\ G(1) = \sqrt {(16)^2 + (12)^2} = 20 \\ G(2) = \sqrt {(0.1)^2 + (4)^2} = 4 \\ G(3) = \sqrt {(16)^2 + (2.39)^2} = 16 \\ G(4) = \sqrt {(0.00)^2 + (0.00)^2} = 0 \\ G(5) = \sqrt {(0.00)^2 + (0.00)^2} = 0 \\ G(6) = \sqrt {(0.00)^2 + (0.00)^2} = 0 \\ G(7) = \sqrt {(0.00)^2 + (0.00)^2} = 0 \\ $$

Phase Angle

$$ \\ Phase0 = tan^{-1}{(0/55.99)} = 0 \\ Phase1 = tan^{-1}{(12/16)} = 36.9 \\ Phase2 = tan^{-1}{(4/0.1)} = 88.6 \\ Phase3 = tan^{-1}{(2.39/16)} = 8.5 \\ Phase4 = tan^{-1}{(0/0)} = 0 \\ Phase5 = tan^{-1}{(0/0)} = 0 \\ Phase6 = tan^{-1}{(0/0)} = 0 \\ Phase7 = tan^{-1}{(0/0)} = 0 $$

Fig. 1.12

Fig. 1.13

Fig. 1.14

Using inverse Fourier Transform we can bring back the function

n	X(n)real	X(n)img	g(0)	g(1)	g(2)	g(3)	g(4)	g(5)	g(6)	g(7)	g(8)	g(9)	g(10)	g(11)	g(12)	g(13)	g(14)	g(15)
0	55.99	0.00	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99	55.99
1	16.01	-12.00	16.01	19.38	19.80	17.21	12.00	4.96	-2.84	-10.20	-16.01	-19.38	-19.80	-17.21	-12.00	-4.96	2.84	10.20
2	0.01	-4.00	0.01	2.83	4.00	2.82	-0.01	-2.83	-4.00	-2.82	0.01	2.83	4.00	2.82	-0.01	-2.83	-4.00	-2.82
3	16.00	-2.39	16.00	8.33	-9.63	-15.70	-2.39	13.87	13.00	-3.92	-16.00	-8.33	9.63	15.70	2.39	-13.87	-13.00	3.92
4	0.00	-0.01	0.00	0.01	0.00	-0.01	0.00	0.01	0.00	-0.01	0.00	0.01	0.00	-0.01	0.00	0.01	0.00	-0.01
5	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
6	-0.01	0.00	-0.01	0.01	0.00	0.00	0.01	-0.01	0.00	0.00	-0.01	0.01	0.00	0.00	0.01	-0.01	0.00	0.00
7	-0.01	-0.01	-0.01	0.01	-0.01	0.01	-0.01	0.01	0.00	0.00	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.00	0.00
8	0.01	0.00	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.01	-0.01
9	-0.01	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.01	0.00	0.00	0.01	-0.01	0.01	-0.01	0.01	-0.01	0.00	0.00
10	-0.01	0.00	-0.01	0.01	0.00	0.00	0.01	-0.01	0.00	0.00	-0.01	0.01	0.00	0.00	0.01	-0.01	0.00	0.00
11	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
12	0.00	0.01	0.00	0.01	0.00	-0.01	0.00	0.01	0.00	-0.01	0.00	0.01	0.00	-0.01	0.00	0.01	0.00	-0.01
13	16.00	2.39	16.00	8.33	-9.63	-15.70	-2.39	13.87	13.00	-3.92	-16.00	-8.33	9.63	15.70	2.39	-13.87	-13.00	3.92
14	0.01	4.00	0.01	2.83	4.00	2.82	-0.01	-2.83	-4.00	-2.82	0.01	2.83	4.00	2.82	-0.01	-2.83	-4.00	-2.82
15	16.01	12.00	16.01	19.38	19.80	17.21	12.00	4.96	-2.84	-10.20	-16.01	-19.38	-19.80	-17.21	-12.00	-4.96	2.84	10.20
	$$\sum$$	120.00	117.12	84.32	64.64	75.20	88.00	68.32	22.08	-8.00	6.24	43.68	58.56	36.80	12.64	27.68	78.56
	$$\sum/16$$	7.50	7.32	5.27	4.04	4.70	5.50	4.27	1.38	-0.50	0.39	2.73	3.66	2.30	0.79	1.73	4.91

The Discrete Time Fourier Transform (DFT) is a computationaly intensive process, for sample size $$N$$ there are $$N^2$$ multiply and add add operations, not a very practical solution. J. W. Cooley and John Tukey developed an algorithm called Fast Fourier Transform (FFT) that greatly reduced the number of operations, and is being used widely in computer applications. Here is the simplification developed by Cooley and Tukey.

For a sample size 8 the total number of operations are as follows

$$ \\ X(0) = x(0)e^{-j{2\pi 0\over 8}0} + x(1)e^{-j{2\pi 0\over 8}1} + x(2)e^{-j{2\pi 0\over 8}2} + x(3)e^{-j{2\pi 0\over 8}3} + x(4)e^{-j{2\pi 0\over 8}4} + x(5)e^{-j{2\pi 0\over 8}5} + x(6)e^{-j{2\pi 0\over 8}6} + x(7)e^{-j{2\pi 0\over 8}7} \\ X(1) = x(0)e^{-j{2\pi 1\over 8}0} + x(1)e^{-j{2\pi 1\over 8}1} + x(2)e^{-j{2\pi 1\over 8}2} + x(3)e^{-j{2\pi 1\over 8}3} + x(4)e^{-j{2\pi 1\over 8}4} + x(5)e^{-j{2\pi 1\over 8}5} + x(6)e^{-j{2\pi 1\over 8}6} + x(7)e^{-j{2\pi 1\over 8}7} X(2) = x(0)e^{-j{2\pi 2\over 8}0} + x(1)e^{-j{2\pi 2\over 8}1} + x(2)e^{-j{2\pi 2\over 8}2} + x(3)e^{-j{2\pi 2\over 8}3} + x(4)e^{-j{2\pi 2\over 8}4} + x(5)e^{-j{2\pi 2\over 8}5} + x(6)e^{-j{2\pi 2\over 8}6} + x(7)e^{-j{2\pi 2\over 8}7} \\ X(3) = x(0)e^{-j{2\pi 3\over 8}0} + x(1)e^{-j{2\pi 3\over 8}1} + x(2)e^{-j{2\pi 3\over 8}2} + x(3)e^{-j{2\pi 3\over 8}3} + x(4)e^{-j{2\pi 3\over 8}4} + x(5)e^{-j{2\pi 3\over 8}5} + x(6)e^{-j{2\pi 3\over 8}6} + x(7)e^{-j{2\pi 3\over 8}7} \\ X(4) = x(0)e^{-j{2\pi 4\over 8}0} + x(1)e^{-j{2\pi 4\over 8}1} + x(2)e^{-j{2\pi 4\over 8}2} + x(3)e^{-j{2\pi 4\over 8}3} + x(4)e^{-j{2\pi 4\over 8}4} + x(5)e^{-j{2\pi 4\over 8}5} + x(6)e^{-j{2\pi 4\over 8}6} + x(7)e^{-j{2\pi 4\over 8}7} \\ X(5) = x(0)e^{-j{2\pi 5\over 8}0} + x(1)e^{-j{2\pi 5\over 8}1} + x(2)e^{-j{2\pi 5\over 8}2} + x(3)e^{-j{2\pi 5\over 8}3} + x(4)e^{-j{2\pi 5\over 8}4} + x(5)e^{-j{2\pi 5\over 8}5} + x(6)e^{-j{2\pi 5\over 8}6} + x(7)e^{-j{2\pi 5\over 8}7} \\ X(6) = x(0)e^{-j{2\pi 6\over 8}0} + x(1)e^{-j{2\pi 6\over 8}1} + x(2)e^{-j{2\pi 6\over 8}2} + x(3)e^{-j{2\pi 6\over 8}3} + x(4)e^{-j{2\pi 6\over 8}4} + x(5)e^{-j{2\pi 6\over 8}5} + x(6)e^{-j{2\pi 6\over 8}6} + x(7)e^{-j{2\pi 6\over 8}7} \\ X(7) = x(0)e^{-j{2\pi 7\over 8}0} + x(1)e^{-j{2\pi 7\over 8}1} + x(2)e^{-j{2\pi 7\over 8}2} + x(3)e^{-j{2\pi 7\over 8}3} + x(4)e^{-j{2\pi 7\over 8}4} + x(5)e^{-j{2\pi 7\over 8}5} + x(6)e^{-j{2\pi 7\over 8}6} + x(7)e^{-j{2\pi 7\over 8}7} \cdots (2.12) \\ \\ \\ $$

Let us Group Odd and Even operation for one frequency

$$ \\ \\ X(3) = { \begin {pmatrix} x(0)e^{-j{2\pi 3\over 8}0} + x(2)e^{-j{2\pi 3\over 8}2} + x(4)e^{-j{2\pi 3\over 8}4} + x(6)e^{-j{2\pi 3\over 8}6} { \end {pmatrix}} + { \begin {pmatrix} + x(1)e^{-j{2\pi 3\over 8}1} + x(3)e^{-j{2\pi 3\over 8}3} + x(5)e^{-j{2\pi 3\over 8}5} + x(7)e^{-j{2\pi 3\over 8}7} { \end {pmatrix}} \\ \\ $$

Simplify the expression using $$e^{-j{2\pi k\over N}0} = 1$$ and factroring $$e^{-j{2\pi 3\over 4}n} = e^{-j{2\pi 3\over 4}(n-1)} e^{-j{2\pi 3\over 4}1}$$

$$ \\ \\ X(3) = { \begin {pmatrix} x(0) + x(2)e^{-j{2\pi 3\over 4}1} + x(4)e^{-j{2\pi 3\over 4}2} + x(6)e^{-j{2\pi 3\over 4}3} { \end {pmatrix}} + e^{-j{2\pi 3\over 8}1} { \begin {pmatrix} x(1) + x(3)e^{-j{2\pi 3\over 4}1} + x(5)e^{-j{2\pi 3\over 4}2} + x(7)e^{-j{2\pi 3\over 4}3} { \end {pmatrix}} \\ \\ $$

Repeat the regroup Odd and Even operations until there is only one odd and one even operation in the group

$$ \\ \\ X(0) = { \begin {pmatrix} x(0) + x(4)e^{-j{2\pi 0\over 2}1} { \end {pmatrix}} + e^{-j{2\pi 0\over 4}1} { \begin {pmatrix} x(2)} + x(6)e^{-j{2\pi 0\over 2}1} { \end {pmatrix}} + e^{-j{2\pi 0\over 8}1} { \begin {pmatrix} { \begin {pmatrix} x(1) + x(5)e^{-j{2\pi 0\over 2}1} { \end {pmatrix}} + e^{-j{2\pi 0\over 4}1} { \begin {pmatrix} x(3)} + x(7)e^{-j{2\pi 0\over 2}1} { \end {pmatrix}} { \end {pmatrix}} \\ \\ \vdots \\ \\ X(7) = { \begin {pmatrix} x(0) + x(4)e^{-j{2\pi 7\over 2}1} { \end {pmatrix}} + e^{-j{2\pi 7\over 4}1} { \begin {pmatrix} x(2)} + x(6)e^{-j{2\pi 7\over 2}1} { \end {pmatrix}} + e^{-j{2\pi 3\over 8}1} { \begin {pmatrix} { \begin {pmatrix} x(1) + x(5)e^{-j{2\pi 7\over 2}1} { \end {pmatrix}} + e^{-j{2\pi 3\over 4}1} { \begin {pmatrix} x(3)} + x(7)e^{-j{2\pi 7\over 2}1} { \end {pmatrix}} { \end {pmatrix}} \\ \\ $$

Recognizing the fact that $$ e^{-j{2\pi 3\over 2}n} = -1 $$ the last operation is simply multiply by -1,

$$ \\ \\ X(0) = { \begin {pmatrix} x(0) - x(4)} { \end {pmatrix}} + e^{-j{2\pi 0\over 4}1} { \begin {pmatrix} x(2)} - x(6)} { \end {pmatrix}} + e^{-j{2\pi 0\over 8}1} { \begin {pmatrix} { \begin {pmatrix} x(1) - x(5)} { \end {pmatrix}} + e^{-j{2\pi 0\over 4}1} { \begin {pmatrix} x(3)} - x(7)} { \end {pmatrix}} { \end {pmatrix}} \\ \\ \vdots \\ \\ X(7) = { \begin {pmatrix} x(0) - x(4)} { \end {pmatrix}} + e^{-j{2\pi 7\over 4}1} { \begin {pmatrix} x(2)} - x(6)} { \end {pmatrix}} + e^{-j{2\pi 7\over 8}1} { \begin {pmatrix} { \begin {pmatrix} x(1) - x(5)} { \end {pmatrix}} + e^{-j{2\pi 7\over 4}1} { \begin {pmatrix} x(3)} - x(7)} { \end {pmatrix}} { \end {pmatrix}} \\ \\ $$

The whole operation lends itself to a nice recursve algorithm with simple multiply and add operation, The following code snippet is a C++ program implementing FFT algorithm and later on we will replace the floating point add and multiply with a fixed point arithmatic that is guaranteed to improve the performance where Floating Point CoProcessor is not avaialable: