## Digital signal processing in Scilab and Xcos. Part 3: Fourier Analysis

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] = ∑ xk δk[n] = ∑ xk δ[n−k] k=0 k=0

where δ[nk] 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 ∑ Xk ψk[n] k=0

where ψk[n] = exp (jnk/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 kh

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

The Fourier coefficients are

 N−1 Xk = ⟨ψ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

X0/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 δ[nk] (click on the picture to enlarge):

We see that the Fourier coefficients of δ[nk] are equal to the components of the vectors ψ*k

Conversely, the Fourier coefficients of ψk will be the components of the vectors N δ[nk]

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