As it known from Fourier series, any signal can be represented as a sum of different simple signals, which basically allow as to simplify a computation and makes possible FFT transformation.

To begin to dig into basics let me refresh mind with couple fundamentals equations, which help to work with Fourier.

Math
Euler equations: e^{jx} = cos(x)+jsin(x) cos(x)=\frac{e^{jx}+e^{-jx}}{2} sin(x)=\frac{e^{jx}-e^{-jx}}{2j} Trigonometry: cos(x)cos(y)=0.5(cos(x-y)+cos(x+y)) sin(x)sin(y)=0.5(cos(x-y)-cos(x+y)) cos^2(x)=0.5(1+cos(2x)) \int_{-\pi}^{\pi} cos^2(x) \,dx=\int_{-\pi}^{\pi} {(1+cos(2x))} \,dx=0.5(x+sin(x)cos(x))\Biggr|_{-\pi}^{\pi}=\pi cos(-x)=cos(x) - even\ function, sin(-x)=-sin(x) - odd\ function

OK, the theoretical basis is built, we can write down some basics of Fourier series.

FOURier series

Any signal can be represented as: f(x)=\sum_{n=0}^{\infty} a_ncos(nx) + \sum_{n=1}^{\infty} b_nsin(nx) How to calculate these coefficients? Use an ortogonal funciton! As known, if we integrate sin(x)cos(x) from -pi to pi, the result will be equal to 0, so that can be used for the calculation of coefficients. If we want to calculate the a coefficient, let’s multiply the function to cos(kx), in that case the integral will be non-zero only if k=n, so: \int_{-\pi}^{\pi} f(x)cos(kx) \,dx=\int_{-\pi}^{\pi} a_kcos^2(kx) \,dx=a_k\pi So, the coefficient is equal to: a_k=\frac{1}{\pi}\int_{-\pi}^{\pi} f(x)cos(kx) \,dx b_k=\frac{1}{\pi}\int_{-\pi}^{\pi} f(x)sin(kx) \,dx Now need to check what do we have when k equal to 0 – b_k will be equal to 0, while a_k is not: a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x) \,dx

Example 1: The delta (Dirac) function

\delta(x)=\infty\ for\ x=0\ and\ \delta(x)=0\ for\ x<>0

Integral of Dirac function is equal to 1. This is important. Let’s look closer – that function is the even function, so it doesn’t contain the sine portion, only the cosine is presented.

a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi} \delta(x) \,dx=\frac{1}{2\pi}

For other coefficient it is also should be noticed that if we multiply the delta function by any other function in the area different from 0, the result will be 0, and in 0 area cosine is equal to 1, so:

a_k=\frac{1}{\pi}\int_{-\pi}^{\pi} \delta(x)cos(kx) \,dx=\frac{1}{\pi}

Which basically means that all coefficients are equal to 1/pi, and a resulting representation is:

\delta(x)=\frac{1}{2\pi}+\frac{1}{\pi}(cosx + cos(2x) + cos(3x) + cos(4x) + ... )

More harmonic we include, closer to ideal curve we move:

What I, personally still don’t understand why we have that 1/2pi term, it changes value for non-zero x moving it away from ideal form.

Example 2: The sign function

Sign function is equal to +1 from 0 to pi, and equal to -1 from -pi to 0:

Looking at that plot we already can say that this is an ODD function. So we don’t have cosine components, and, one more useful property of the odd function – we can integrate only half of the period and multiply result by two to get the full value.

b_k=\frac{2}{\pi}\int_{0}^{\pi} f(x)sin(kx) \,dx=\frac{2}{\pi}\int_{0}^{\pi} sin(kx) \,dx=\frac{2}{\pi}(\frac{-cos(kx)}{k})\Biggr|_{0}^{\pi}=\frac{2}{\pi}(\frac{2}{1};\frac{0}{2};\frac{2}{3};\frac{0}{4};\frac{2}{5};...)

Now, lets instantiate it into Fourier series basic representation:

\frac{4}{\pi}(sin(x)+sin(3x)+sin(5x)+...)=\frac{4}{\pi}\sum_{k=1}^{\infty} \frac{sin((2k-1)x)}{2k-1}

Only odd components are survived, like for usual square wave signal frequency content.

But what if our function will be phase shifted and turned from odd into the even function?

So now the function is even and will be represented by cosine terms. Lets start from a0:

a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x) \,dx=\frac{1}{2\pi}\int_{-\pi}^{\pi} (-1) \,dx=(-x)\Biggr|_{-\pi}^{\pi}=0

I actually think that a_0 calculation is wrong there, we actually need to take integrals from all 3 areas: -pi to -pi/2, -pi/2 to pi/2, pi/2 to pi. This will lead to the same result, still 0, so we can skip that. But we also can not do integral from -pi to pi for ak coefficients, because function is changed within that period, so we also can do that division by 3 areas.

a_{k*}=\frac{1}{\pi}\int_{-\pi}^{\frac{-\pi}{2}} f(x)cos(kx) \,dx=\frac{1}{\pi}(-\frac{sin(kx)}{k}\Biggr|_{-\pi}^{\frac{-\pi}{2}})=\frac{1}{\pi}(1;\frac{-1}{3};\frac{1}{5};...)

a_{k**}=\frac{1}{\pi}\int_{\frac{\pi}{2}}^{\pi} f(x)cos(kx) \,dx=\frac{1}{\pi}(-\frac{sin(kx)}{k}\Biggr|_{\frac{\pi}{2}}^{\pi})=\frac{1}{\pi}(1;\frac{-1}{3};\frac{1}{5};...)

a_{k***}=\frac{1}{\pi}\int_{\frac{-\pi}{2}}^{\frac{\pi}{2}} f(x)cos(kx) \,dx=\frac{1}{\pi}(\frac{sin(kx)}{k}\Biggr|_{\frac{-\pi}{2}}^{\frac{\pi}{2}})=\frac{2}{\pi}(1;\frac{-1}{3};\frac{1}{5};...)

Now, all we need to do just to sum all 3 parts:

a_{k*}+a_{k**}+a_{k***}=\frac{4}{\pi} \sum_{k=1}^{\infty}\frac{cos((2k-1)x)}{(2k-1)(-1)^(k+1)}

Let’s plot that function:

Looks similar to our initial function, so seems I did my calculation right.

Complex form of Fourier series.

Looking at coefficient derived in the first part it can be seen that a_0 is divided by two in comparison to usual (why?)

f(x)=\frac{a_0}{2} + \sum_{n=1}^{\infty} (a_ncos(nx) + b_nsin(nx))=\frac{a_0}{2} + \sum_{n=1}^{\infty}\frac{a_n-jb_n}{2}e^{jnx} + \sum_{n=1}^{\infty}\frac{a_n+jb_n}{2}e^{-jnx} = \sum_{-\infty}^{\infty}c_ne^{jnx}

c_0=\frac{a_0}{2}, c_n=\frac{a_n-jb_n}{2}, c_{-n}=\frac{a_n+jb_n}{2}

So we have three coefficient, the zero one and two complex conjugates.

c_0=\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x) \,dx

Coefficients would be:

c_n=\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)e^{-jnx} \,dx, n=0,\pm1,\pm2,\pm3...

Let’s try to calculate Fourier series in complex form for sign function. C_0 component will still be 0, it is easy to check by integration from -pi to 0 plus integral from 0 to pi.

The same sum of integrals should be taken for c_n as well:

[katex]c_n=\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)e^{-jnx} \,dx=\frac{1}{2\pi}(\int_{-\pi}^{0} (-1)e^{-jnx} \,dx+\int_{0}^{\pi} e^{-jnx} \,dx)=\frac{1}{n\pi}(\frac{e^{jn\pi}+e^{-jn\pi}}{2}-1)=\frac{1}{n\pi}(cos(n\pi)-1)=\frac{1}{n\pi}((-1)^n-1)

If n is even, the coefficient is equal to 0, only when it is odd (equal to 2k-1):

c_{2k-1}=-\frac{2j}{(2k-1)\pi}

And the function can be represented as:

f(x)=sign(x)=-\frac{2j}{\pi}\sum_{k=-\infty}^{\infty}\frac{e^{j(2k-1)x}}{2k-1}

Lets try to check in Octave, if the function we derived still equal to sign shape:

Looks similar, let's make something else I think.

f(x)=sin(8x)+sin(9x)+cos(x)

Now the tough part - that function requires solving of the complex integral, which appeared not quite straightforward for Octave. Python looks like using the same quad algorithm came from Fortran. Matlab has that ability, but I dont have an ability to buy Matlab license 🙂

So the story is that I ended up with using of a webservice wolfram alpha, and hand calculated coefficients up to 9, still easier than hand calculate on the paper huh?

The code for octave looks like:

clf;
clear;
st = pi/256;
step = -pi:st:pi;
y = (sin(8*step)+sin(9*step)+cos(step));

plot(step,y)
hold on;


n = 9;

c_p(1) = 0.5;
c_p(2) = 0;
c_p(3) = 0;
c_p(4) = 0;
c_p(5) = 0;
c_p(6) = 0;
c_p(7) = 0;
c_p(8) = -0.5*j;
c_p(9) = -0.5*j;


c_n(1) = 0.5;
c_n(2) = 0;
c_n(3) = 0;
c_n(4) = 0;
c_n(5) = 0;
c_n(6) = 0;
c_n(7) = 0;
c_n(8) = 0.5*j;
c_n(9) = 0.5*j;

y_aprox = 0;

for i = 1:n
  y_aprox = y_aprox + c_p(i)*exp(j*i*step) + c_n(i)*exp(-j*i*step);
endfor

plot(step, y_aprox)
legend("original", "approximated")
hold on
grid on

It is enough to fully represent the curve:

Looks right, but it is a pity that octave doesn't allow to make a complex integration. Or maybe not, if it will not impact the FFT by itself, will see.