Sunday, October 17, 2010

Signal processing in sagemath

I was experimenting with some tools for signal processing. The first tools that come to mind are opensource tools scilab and octave. I spent a day looking at the tools, the UI. I figured I could do definite integrals quite easily. Having spent the day and not having the made the progress I would have liked to, I turned to sagemath.

My first delight was to learn about piecewise functions. I decided to play with a square wave of time period 20.

I started with a function that looked like

f1(x) = 1
f2(x) = 0
f = Piecewise([[(0,10),f1],[(10,20),f2]])

Plotting the function with plot(f) showed






Which is exactly what I wanted. Piecewise class also supports a number of very useful functions related to fourier series.

One particularly useful one is a plot fourier series partial sum. The function shows how as we add more frequencies the fourier approximation to the original wave gets better

I started with

f.plot_fourier_series_partial_sum(5,10,-20, 20)

f.plot_fourier_series_partial_sum(15,10,-20, 20) gave me



f.plot_fourier_series_partial_sum(150,10,-20, 20) gave me



This looks like a good approximation to the wave we started with and almost looks like the signals I got on my oscilloscope during my engineering days :)

The next important thing was to get the sine and cosine coefficients

print "sine terms"
for j in range(0,21):
    pretty_print(f.fourier_series_sine_coefficient(j, 10))
print "cosine terms"
for j in range(0,21):
    pretty_print(f.fourier_series_cosine_coefficient(j, 10))



sine terms
0
2/pi
0
2/3/pi
0
2/5/pi
0
2/7/pi
0
2/9/pi
0
2/11/pi
0
2/13/pi
0
2/15/pi
0
2/17/pi
0
2/19/pi
0
cosine terms
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
 
 
Overall, I had good day with sagemath. I need to experiment with some of the discrete functions. Keep tuned in, I'll try and keep you updated on how it goes.
Post a Comment