Lagrange polynomial interpolation in Scilab

Assume we have data (xi, yi), i = 1, …, n.

The data don’t have to be equally spaced.

We can pass a Lagrange polynomial P(x) of degree n−1 through these data points.


function y = lagrange_interp(x,data)
    
    for i = 1:length(x)
        y(i) = P(x(i),data);
    end
    
endfunction

The polynomial P(x) is a linear combination of polynomials Li(x), where each Li(x) is of degree n−1

P(x) = y1L1(x) + y2L2(x) + … + ynLn(x)

The polynomials Li(x) are defined as

Li(x) = (xx1) … (xxi-1)(xxi+1) … (xxn) / αi

where

αi = (xix1) … (xixi-1)(xixi+1) … (xixn)

It is obvious that Li(xj) = δij and therefore P(xi) = yi meaning that the polynomial P passes through all the data (xi, yi), i = 1, …, n.

If the the data (xi, yi) are organized in an n×2 matrix, then the polynomial P can be computed as follows


function y = P(x,data)

    n = size(data,1);

    xi = data(:,1);
    yi = data(:,2);

    L = cprod_e(xi,x) ./ cprod_i(xi); 
    y = yi' * L;

endfunction

where the function cprod_e calculates the numerators of Li(x) for each i, namely all the products

(xx1) … (xxi-1)(xxi+1) … (xxn)


function y = cprod_e(x,a)
    
    n = length(x);
    y(1) = prod(a-x(2:$));
    for i = 2:n-1
        y(i) = prod(a-x(1:i-1))*prod(a-x(i+1:$));
    end
    y(n) = prod(a-x(1:$-1));
    
endfunction

The function cprod_i calculates all αi for i = 1, …, n.


function y=cprod_i(x)
    
    n = length(x);
    y(1) = prod(x(1)-x(2:$));
    for i = 2:n-1
        y(i) = prod(x(i)-x(1:i-1))*prod(x(i)-x(i+1:$));
    end
    y(n) = prod(x(n)-x(1:$-1));
    
endfunction

Now, let’s test our code.

I take two examples from the book “Fundamentals of Engineering Numerical Analysis” by Prof. Parviz Moin.

Both examples use data obtained from the Runge’s function

y = 1/(1+25 x2)

The data in the first example are equally spaced:


data1=[...
-1.0 0.038;
-0.8 0.058;
-0.6 0.1;
-0.4 0.2;
-0.2 0.5;
0.0 1.0;
0.2 0.5;
0.4 0.2;
0.6 0.1;
0.8 0.058;
1.0 0.038];

The commands to perform Lagrange interpolation are:


x = linspace(-1,1,1000);
y = lagrange_interp(x,data1);

Here’s the plot of the data (red circles), the Lagrange polynomial (solid line) and the Runge’s function (dashed line):


plot(x,y,'-b',...
x, 1.0 ./ (1+25.*x.^2),'--g',...
data1(:,1),data1(:,2),'or')

L_eg_1

The second example takes non-equally spaced data:


data2 = [...
-1.0 0.038;
-0.95 0.042;
-0.81 0.058;
-0.59 0.104;
-0.31 0.295;
0.0 1.0;
0.31 0.295;
0.59 0.104;
0.81 0.058;
0.95 0.042;
1.0 0.038];

Interpolating and plotting them as we did in the previous example produces the following picture.


y = lagrange_interp(x,data2);

plot(x,y,'-b',...
x, 1.0 ./ (1+25.*x.^2),'--g',...
data2(:,1),data2(:,2),'or')

L_eg_2

Published in: on 04/07/2014 at 08:18  Leave a Comment  
Tags:

Measuring banana radioactivity

Banana equivalent dose

Banana equivalent dose: link

We did it with this kit: link

A video from the kit manufacturer’s site:

Published in: on 02/07/2014 at 14:58  Leave a Comment  

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

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⟩ = ejnk/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):

psi

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

delta

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.

fin

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  
Tags:

Digital signal processing in Scilab and Xcos. Part 2: Playing sound files

In my previous post on DSP in Scilab, I gave an example on how to open a wav file and make it available for further processing.

Today we’ll discuss the opposite task: how to generate signals and store them into wav files.

Open the Scilab’s editor


editor

Choose the sampling rate, the number of bits per sample, and the time duration for your signals


Fs = 11025; // samples per second
bits = 16; // bits per sample
t_total = 10; // seconds

The total number of samples in the signal:


n_samples = Fs * t_total;

Create an array of time points at which the samples will be synthesized:


t = linspace(0, t_total, n_samples);

The following code will generate a 440 Hz sine wave, and save it into a wav file sin440.wav


f=440; // sound frequency

// Sine wave
sin_wave = sin(2*%pi*f*t);
sin_file = "C:\Users\Vlad\wav\sin"+string(f)+".wav";
wavwrite(sin_wave, Fs, bits, sin_file);

Let’s now generate a sawtooth wave of the same frequency and output it into another wav file:


// Sawtooth wave
saw_wave=2*(f*t-floor(0.5+f*t));
saw_file = "C:\Users\Vlad\wav\saw"+string(f)+".wav";
wavwrite(saw_wave, Fs, bits, saw_file);

Similarly, for a triangle and a square waves:


// Triangle wave
tri_wave=(2/%pi)*asin(sin(2*%pi*f*t));
tri_file = "C:\Users\Vlad\wav\tri"+string(f)+".wav";
wavwrite(tri_wave, Fs, bits, tri_file); 

// Square wave
sq_wave=sign(sin(2*%pi*f*t));
sq_file = "C:\Users\Vlad\wav\sq"+string(f)+".wav";
wavwrite(sq_wave, Fs, bits, sq_file); 

Let’s simulate a guitar sound using the Karplus–Strong string synthesis method.


// Karplus-Strong
n_width=100;
ks=-1+2*rand(1,n_width,"uniform");
alpha=0.96;

while ( length(ks) < n_samples )
    ks=[ks,alpha*ks($-n_width+1:$)];
end
ks=ks(1:n_samples);

We can play the generated sound in Scilab with


sound(ks)

and plot its waveform


plot(t,ks)

ks

Next: part 3

Published in: on 08/05/2014 at 19:36  Leave a Comment  
Tags:
Follow

Get every new post delivered to your Inbox.

Join 73 other followers