How to Use Scilab: fft

Showing posts with label fft. Show all posts

Thursday, October 18, 2012

Scilab as software for engineers

I'm an engineer (of telecommunications) and have used Scilab for several years (since 2005). There are many useful tools for image and signal processing, optimization, linear systems and transformations, and general simulations.

Scilab is a software for numerical simulation supported by Scilab Consortium. This group of companies provides resources and services for users and other companies those intend to migrate to an open and stable software platform.

Beyond formal support, Scilab Consortium is an open environment for doubts and discussions where everyone can help and be helped.

About technical resources, Scilab has two ways for developing and application or simulation: script language (Scinotes) or drag-and-drop interface (Scicos), both are powerful but I believe script language is supported by more external toolboxes.

Scilab external toolboxes are freeware and can be downloaded and installed through ATOMS tool (see next picture, click in the figure to enlarge it).

Monday, March 15, 2010

Discrete Fourier Transform - DFT

I made some posts about Discrete Fourier Transform (DFT), they're here.

But, I received a comment asking me how to do a function that implements DFT.

Let's only write a script, not optimized, that calculates DFT:

function X = DFT(x)
n = length(x); //number of elements in 'x'

omega = exp(-2*%pi*%i/n);
j=0:(n - 1);
F=omega.^(j'*j); //Fourier matrix

X = F*x(:); //X = DFT(x)
endfunction;

If anyone wants to try the showed function, you could do like this:

N = 10;
x = rand(N, 1);

X1 = DFT(x)
X2 = dft(x, -1) // dft(.) function uses flag = -1 for direct transform and flag = 1 for inverse transform
X3 = fft(x)

Now, look the values of X1, X2 and X3.

Did I help?

Thursday, November 5, 2009

More frequencies in FFT

Hi Transmogrifox, thank you very much for your help. We have this blog to interact and to help ourselves.

The Transmogrifox's code is given following:

//Frequency components of a signal
//----------------------------------
// build a noides signal sampled at 1000hz containing to pure frequencies
// at 50 and 70 Hz
sample_rate=1000;
t = 0:1/sample_rate:0.6;
N=size(t,'*'); //number of samples
s=sin(2*%pi*50*t)+ 0.3*sin(2*%pi*70*t+%pi/4)+ 0.2*grand(1,N,'nor',0,1);

y=fft(s);
//the fft response is symetric we retain only the first N/2 points
f=sample_rate*(0:(N/2))/N; //associated frequency vector
n=size(f,'*')
yy = y*2/N;
clf()
plot2d(f,abs(yy(1:n)))

And the following picture is the result.

So, if anyone has any question, I think I and Transmogrifox are ready for help.

Friday, October 23, 2009

Peaks in FFT transformation

Hi Kaustubh, excuse me for the time of the answer.

Do you want to know about the peaks in the FFT transformation (here), right?

Well, this blog is about Scilab, and you can find what you want in any book of Digital Signal Processing.

But God loves you and I'll try to explain what the peaks mean.

Think in a pure frequency signal x(t) like a sine or cosine, it has only the own frequency (f0). The FFT transform is two impulses (X(f) = d(-f0) + d(f0)) (similar to this figure)

If you plot the FFT transform of x(t), then you obtain the peaks in -f0 and f0.

Now, think in a signal composed by two pure frequency signals, y(t) = x1(t) + x2(t), of different frequencies (f1 and f2).

The FFT transform of y(t) is four impulses, two for x1(t) and two for x2(t).

If x1(t) and x2(t) have different amplitudes, for example: x1(t) = cos(2t) and x2(t) = 3 cos(5t), then the peaks of the FFT transform of y(t) are weighted,

Following the example, the FFT transform of x1(t) is X1(f) = d(-f1) + d(f1) and the FFT transformation of x2(t) is X2(t) = 3(d(-f2) + d(f2)). So, the FFT transform of y(t) = x1(t) + x2(t) is Y(f) = X1(f) + X2(f) = (d(-f1) + d(f1)) + 3(d(-f2) + d(f2)).

Now, if you have a real signal xr(t), it has many frequencies and it FFT transform XR(f) presents many peaks, one by each frequency. The peaks' amplitude represents the amplitude of each frequency in xr(t).

My friend, I said this not a blog about Signal Processing, but I tried to help you.
If you have more question, I'll try to help you again.

Monday, September 7, 2009

Discrete convolution properties - 2

Anyone called Younes asked something about convolution via Fourier Transform.
I made the following deduction for you Younes (click on the image for see it in the real size):

Now, I ask: try to confirm the property using Scilab.

Don't forget computers don't make infinity operations.

Monday, August 10, 2009

FFT - specifying the frequencies

I'd like to know my readers, but the "no name" readers deserve respect, as everybody.

This post is about the FFT function, and anyone wants to know how to specify the frequencies for plot the values.

A few of theory:

If your signal has N values [0, N - 1], then the Fourier Transform has N distinct values.

The function fft(.) returns the signal in the interval [0, N - 1], so you can use the function fftshift(.), over the function fft(.) like this: X = fftshift(fft(x)), that it returns the Fourier Transform in the intervals:

[-(N - 1)/2, (N - 1)/2], if N is odd.
[-N/2, N/2 - 1], if N is even.

About the frequencies, if your signal was sampled with a rate T (T samples by second), the indexes are given multiplying the interval by: 2*%pi*T.

Look the example:

N = 100;

n = 1:N;

T = 0.1; // one sample by each 0.1s

w = 0.5; // frequency of the sampled signal (in radians)

x = cos(w*n);

plot(n*T, x); // plot the signal indexed by seconds

f = [-N/2:N/2 - 1];

X = fftshift(fft(x));

scf();
plot(f*2*%pi*T, X); // plot the signal indexed by radians

The result is:

Observing that w = 0.5 is the frequency of the sampled signal, the true frequency (of the analog source signal) is given by w_source = w/T = 0.5/0.1 = 5 rad/s. Now, click on the image and look where the peaks are.

Saturday, August 1, 2009

Fast Fourier Transform - FFT

This post is about a good subject in many areas of engineering and informatics: the Fourier Transform. The continuous Fourier Transform is defined as:

f(t) is a continuous function and F(w) is the Fourier Transform of f(t).

But, the computers don't work with continuous functions, so we should use the discrete form of the Fourier Transform:

f[n] is a discrete function of N elements, F[p] is a discrete and periodic function of period N, so we calculate just N (0 to N - 1) elements for F[p].

Who studies digital signal processing or instrumentation and control knows the utilities of this equation.

Now, how to use the Fourier Transform in Scilab?

If we are using large signals, like audio files, the discrete Fourier Transform is not a good idea, then we can use the fast Fourier Transform (used with discrete signals), look the script:

-->N = 100; // number of elements of the signal

-->n = 0:N - 1;

-->w1 = %pi/5; // 1st frequency

-->w2 = %pi/10; // 2nd frequency

-->s1 = cos(w1*n); // 1st component of the signal

-->s2 = cos(w2*n); // 2nd component of the signal

-->f = s1 + s2; // signal

-->plot(n, f);

The result is:

Now, let's study the Fourier Transform of our signal.

-->F = fft(f); // it calculates the Fourier Transform

-->F_abs = abs(F); // F_abs is the absolute value of each element of F

-->plot(n, F_abs);

The result is:

Look at the two graph's peaks, one for the component s1 and the other for the component s2.

The Fourier Transform is a linear transformation, thus it has a inverse transformation: the Inverse Fourier Transform.

Scilab has the function ifft(.) for obtain the original signal from it Fourier Transform.

If anyone wants to know, I can make a new post about how to identify the frequencies of the original signal in the Fourier Transform.

How to Use Scilab

Pages