[Py] Spline function test

To interpolate data points, Scipy has several modules for this work. For what I am interested in, there are two functions seem to meet my requirement:

• scipy.interpolate.UnivariateSpline
• scipy.interpolate.InterpolatedUnivariateSpline

According to the explanations given by the above links, to use UnivariateSpline with s=1 is equivalent to use InterpolatedUnivariateSpline.

Here is the results of a simple test:

Code:

from numpy import linspace,exp
from numpy.random import randn
from scipy.interpolate import UnivariateSpline, InterpolatedUnivariateSpline
from pylab import plot, legend, show

x = linspace(-3, 3, 100)
y = exp(-x**2) + randn(100)/10

s0 = UnivariateSpline(x, y, s=0)
s1 = UnivariateSpline(x, y, s=1)
s2 = InterpolatedUnivariateSpline(x, y)

xs = linspace(-3, 3, 1000)
ys0 = s0(xs)
ys1 = s1(xs)
ys2 = s2(xs)

plot(xs,ys0,'g--',label='UnivariateSpline s=0',linewidth=4)
plot(xs,ys1,'b',label='UnivariateSpline s=1')
plot(xs,ys2,'k',label='InterpolatedUnivariateSpline')
plot(x,y,'ro',label='original signal')
legend()
show()

Producing sinusoidal signals

Sometimes I need a signal of different frequencies for testing, so I wrote a simple Python code to produce compound sinusoidal signal by adding several cosine waves of different frequencies, amplitudes, and phases.

The code is listed as the follows. (download)

from scipy import *
from pylab import *

data_size = 10 # number of sinusoidal waves

A = rand(data_size)
theta = rand(data_size)
f = [i for i in range(1,data_size+1)]
t = linspace(0,10,1024)

y = [0.0 for i in range(len(t))]

for i in range(data_size):
        y += A[i]*cos(2*pi*f[i]*t + theta[i])

plot(t,y)
show()

Trials of EMD C code

I found a simple module of EMD in C on google code: realtime-emd. The source code, however, has no main function to invoke the modules. So I wrote one as the follows.

#include <stdio>
#include "EmpiricalModeDecomposition.h"

void imfs_output(emdData* emd){
        int i, j;

        for(j = 0; j < emd->size; j++){
                for(i = 0; i < emd->order; i++) {
                        printf("%f ", emd->imfs[i][j]);
                }   
                printf("\n");
        }   
}

int main(void){
        FILE *inputfile = fopen("sinusoid.txt","r");
        int i=0;
        float signal[1024]={0.0};
        emdData data;

        while (fscanf(inputfile, "%f", &signal[i]) != EOF){
                //printf("%d %f\n",i,signal[i]);
                i++;
        }

        emdCreate(&data,i,5,10,0); // problems here
        emdDecompose(&data,signal);
        imfs_output(&data);
        emdClear(&data);

        fclose(inputfile);
        return 0;
}

The main problem is the arguments of the function emdCreate. I am not sure how to set proper parameters to obtain correct results, so please take this clip of code as just a reference.

You can download the main.c and sinusoid.txt for testing. The input file sinusoid.txt was generate by adding several sinusoidal signals.