Previous parts: part 1, part 2.

Scilab is optimized for matrix calculations, so I will perform Fourier analysis here in the language of matrices and vectors.

For theoretical details about the matrix treatment of Fourier transform, please check Chapter 4 of an excellent book, Signal Processing for Communications, by P.Prandoni and M. Vetterli, 2008, EPFL Press.

Consider finite length complex-valued signals

_{N−1} |
_{N−1} |
|||

x[n] = |
∑ | x_{k} δ_{k}[n] = |
∑ | x_{k} δ[n−k] |

^{k=0} |
^{k=0} |

where δ[*n*−*k*] are the discrete-time impulses:

δ[*n*] = 1 when *n* = 0 and δ[*n*] = 0 when *n* ≠ 0.

The set of discrete-time impulses δ_{k}[*n*] = δ[*n-k*], *k* = 0, …, *N*−1 forms an orthonormal basis in the space of discrete signals of length *N*.

The Fourier transform represents a signal as a superposition of harmonic oscillations of different frequencies.

In other words, Fourier transform of the signal *x*[*n*] is its representation in another basis, namely, in the basis formed by the set of all harmonic oscillators of length *N*

_{N−1} |
||

x[n] = N^{ −1} |
∑ | X_{k} ψ_{k}[n] |

^{k=0} |

where ψ_{k}[*n*] = exp (*j* 2π*nk*/*N*) = (ψ_{1}[1])^{nk}

The vectors ψ_{k}[*n*] form an orthogonal basis in *C* ^{N}

⟨ψ_{k}, ψ_{h}⟩ = *N* when *k* = *h* and ⟨ψ_{k}, ψ_{h}⟩ = 0 when *k* ≠ *h*

Any signal that neither explodes nor disappears has to be composed of periodic signals.

The Fourier coefficients are

_{N−1} |
||

X_{k} = ⟨ψ_{k}, x⟩ = |
∑ | e^{ −j 2πnk/N} x[n] |

^{n=0} |

In matrix form, *X* = Ψ^{†}*x* and *x* = Ψ*X*/*N*, where Ψ_{nk} = (ψ_{1}[1])^{nk}

*X*_{0}/*N* is the average value of the signal *x*[*n*].

The orthonormal basis vectors ψ′_{k}[*n*] = *N*^{ −1/2} ψ_{k}[*n*] are less convenient here because

ψ′_{k}[*n*] ≠ (ψ′_{1}[1])^{nk}

therefore, if we used them, we wouldn’t be able to construct the matrix Ψ in such a beautiful way as we did with the non-orthonormal basis vectors ψ_{k}

A Scilab implementation of the matrix Ψ:

```
function [res] = psi(N)
psi11 = exp(%i*2*%pi/N);
res(1,1:N) = ones(1,N);
for j = 1 : N
res(2,j) = psi11^(j-1);
end
for i = 3 : N
res(i,1:N) = res(2,1:N).^(i-1)
end
endfunction
```

The *k*-th row (or the *k*-th column) of Ψ is the basis vector ψ_{k}

The following picture shows the real and imaginary parts of the basis vectors ψ_{k} with *N* = 64 (click on the picture to enlarge):

The Scilab code for a function that finds the Fourier coefficients:

```
function [res] = F(x)
N = length(x);
res = psi(N)' * x;
endfunction
```

Let’s choose *N* = 64

```
N = 64;
n=[0:1:N-1]';
```

Then, for any signal `x`

, it’s Fourier coefficients are found as

```
X = F(x);
```

I also set those coefficients that are very small to zero, otherwise they will look ugly in pictures

```
eps = 1e-8;
realzero = find(abs(real(X))<eps);
X(realzero) = %i*imag(X(realzero));
imagzero = find(abs(imag(X))<eps);
X(imagzero) = real(X(imagzero));
abszero = find(abs(X)<eps);
X(abszero) = 0;
```

The real parts of the Fourier coefficients are obtained with the function `real(X)`

, their imaginary parts — `imag(X)`

, and their amplitudes — `abs(X)`

.

The phases of the Fourier coefficients are obtained as follows

```
// phase
ph = zeros(X);
nonzero = find(abs(real(X))>eps & abs(imag(X))>eps);
if length(nonzero) > 0 then
ph(nonzero) = atan(imag(X(nonzero)),real(X(nonzero)));
end
```

I put the phase to zero if either the real or imaginary part of *X* is very small.

The following animation shows the Fourier transforms of the basis vectors δ[*n*−*k*] (click on the picture to enlarge):

We see that the Fourier coefficients of δ[*n*−*k*] are equal to the components of the vectors ψ^{*}_{k}

Conversely, the Fourier coefficients of ψ_{k} will be the components of the vectors *N* δ[*n*−*k*]

Now, let’s study Fourier transform of a finite length input.

A finite length signal can be implemented by the following Scilab function:

```
function [res] = fin(n,N)
res = zeros(N,1);
N1=N/2-n/2+1;
res(N1:N1+n-1)=1;
endfunction
```

The animation below shows how the Fourier coefficients change as the length of the signal increases from 1 to 64.

From this animation we see that the longer the signal is in the time domain, the narrower it is in the frequency domain; and vice versa — the shorter the signal is in the time domain, the wider it is in the frequency domain.

A single impulse (of length 1) contains all frequencies in its spectrum while a constant signal *x*[*n*] = 1 has only one frequency (equal to zero).

## Leave a Reply