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
|x[n] =||∑||xk δk[n] =||∑||xk δ[n−k]|
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
|x[n] = N −1||∑||Xk ψk[n]|
where ψk[n] = exp (j 2πnk/N) = (ψ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
|Xk = 〈ψk, x〉 =||∑||e −j 2πnk/N x[n]|
In matrix form, X = Ψ†x and x = ΨX/N, where Ψnk = (ψ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)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 —
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).