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 δ[nk]
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

x[n] = N −1 Xk ψk[n]

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

Xk = ⟨ψk, x⟩ = ejnk/N x[n]

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);
    for i = 3 : N
        res(i,1:N) = res(2,1:N).^(i-1)

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;

Let’s choose N = 64

N = 64;

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

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

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

Published in: on 11/05/2014 at 00:10  Leave a Comment  

The URI to TrackBack this entry is:

RSS feed for comments on this post.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: