Digital Signal Processing with Python Author: Jacek Nosal 2011 r.
Introduction : The aim of this article or should I say tutorial is to present how we can combine different modules and tools created directly for use with Python programming language. I am not going to introduce a general idea of digital signal processing, basically because of the fact that this topic is very complex and it would be difficult to cover all the aspects within a few paragraphs. If You are not familiar with DSP, I suggest You read appropriate book or watch some tutorials. Just to be clear – if I will consider some theory worth presenting, then certain paragraphs will appear. * REMEMBER ! – In order to be able to use all the stuff below, You have to install certain packages : SciPy, NumPy and MatPlotLib for Your Python interpreter – just google them, download and follow instructions.
Why DSP with Python ? Before I start, there is one simple question: what for ? or rather why ? Nowadays, there are many tools and dsp libraries for practically every programming language. So why did I choose Python ? 1. It’s an interpreted programming language, it means that there is a shell (or some gui) and the interpreter reads Your code line by line and it performs an action. It’s much slower than a compiler, but it’s pretty useful, when a problem of debugging/compiling is taken into consideration. Another example of an interpreted language is MATLAB. Why am I not using MATLAB then ? There is a license, that costs a lot of money, and despite the popularity and scale of MATLAB I think that Python is more flexible (if not right now, it will become in no time). 2. Python is a language that is still developing, I mean that there are more and more packages being released day after day and it’s pretty easy to use them. For example, combining matplotlib with scipy or using protocol buffers with Python, or having some fun with PyGame. 3. The last reason is comfort. Once You will learn how to do some DSP with Python You will be able to come up with many interesting things easily. For example, a couple of days ago I used the stuff mentioned above to create random signal generator with ploting functions, histograms and so on. I used Google Protocol Buffers and now when I’m doing some experiments I can easily generate 1000 different signals (it takes about 2 seconds to
generate them and write each signal in separate binary coded file), which are random but I have control over the generation process.
Plan : First of all I will present a NumPy – a package needed for scientific computing with Python. Next, I will cover the basics of matplotlib package and SciPy. Finally, we will combine all of them, write some Python scripts, do some DSP and then watch and discuss results.
NumPy DSP is a process that operates on discrete values of a signal. It means that continuous signals are sampled and instead of continuous set of values we get discrete set. We are doing all the A/C convertion and sampling processes just for our computer’s comfort. Microprocessors operate on binary numbers and they map infinite set of values into a finite one. Shortly, as we are operating on a finite set of numbers – a discrete one it would be nice to use a specific “notation”. So, imagine a vector of indexes that are associated with our samples: Like this: n=[0,1,2,3,...] or n=[1,2,3,....] Right now we are not going to bother ourselves with a value of the first index: should it be 0 or 1 or maybe a different one ? There’s no good answer, but I am starting to index the values from 0. And here’s a vector of our signal’s values: x=[val1,val2,val3,....] My point is: there is a necessity of a flexible, simple and easy to use representation. That’s why we are going to use NumPy. NumPy is a Python package that contains all necessary classes and functions for representing our vectors. This extension allows us to operate on multidimensional arrays and matrices as well. In conslusion: it provides arrays, matrices and also implements all the operations we are allowed to do on them. So let’s have a look on a simple session with python shell. First, let’s import numpy module : from numpy import * Now, we are going to generate some signals using numpy n=linspace(1,100,100) or n=arange(1,101) // linspace syntax is simple linspace(value_min,value_max,how
many) So we’ve created 100 elements from a interval (1;100), if we print our vector the outpul will look like this: print n [1. 2. 3. 4. 5. ... 99. 100.] Now, let’s generate our values: x1=linspace(-1,1,100) //a vector of values from a interval (-1;1), which length is 100 x2=zeros(100) // a vector of zeros, which length is 100 x3=array((x1)*2) // same as x1, but the values are multiplied by 2 => x1.*2 x4=array((x1)**2) // x1 to the power of 2 So, there are some functions that we can use to create our signals and then proceed with our experiments, but let’s have a look at this, in my opinion most realistic example: First let’s generate our sample indexes: n=arange(1,2**10) now let’s create a vector that will represent time dt=1/(sampling_frequency) t=[] for i in range(1,2**10): t.append(n[i-1]*dt) Ok, we have our discrete time, we’ve specified our sampling frequency, now we are able to sample our cosinus wave: x=AMPLITUDE*cos(2*pi*t*signal_frequency) //sometimes we may be obliged to use a loop, because there will be error while multiplying float and int (i don’t really know why it occurs) You may now think that I’ve written so much about NumPy and hardly used it, but in the example above I’ve generated a signal that consists of 2^10 samples. But we when it comes to plotting and working with images in 2-D dimension functions like array, arange, linspace and all the operations on matrices and arrays will be very useful. Consider yourself a user NumPy. But, before You move to the next chapter try typing these commands in Your python shell (remember to import certain packages): a=ones(3)
b=zeros(42) x=arange(100) x=arange(12).reshape(4,3) x=arange(24).reshape(2,3,4) x=array([[1,1],[2,2]]) x*x x=linspace(0,pi,3) y=arange(3) y**2 exp(y) sqrt(y) x=arange(100) x[-1] x[2:5] x[:6:2] n=array([1,3,4,5]) x=arange(10) x[n] xy=vstack([[1,2,3],[4,5,6]]) xy2=hstack([[1,2,3],[4,5,6]]) – don’t skip this one ! Finnaly, if You have some time, need some practice on NumPy then have a look at this website. It’s a NumPy tutorial uploaded by the creators, You can find there everything on NumPy and even more. Of course, You are able to experiment with signals without using a NumPy (just playing with Python in-built tuples,lists and lists of lists etc.) but NumPy is a part of SciPy package so it was created to improve all the operations. My recommendation to You: spend 1 hour on life on learning NumPy syntax and instructions, make some notes, practice in a shell – it’s just an hour and I am sure You will not waste it. (P.S try to generate a few sampled signals like sin,cos waves, sawtooth, square, some exponentials (maybe try
presenting some aliasing) – just a little practice)
Matplotlib Ok, we are able to store some data using python structures or numpy package, so it would be nice to plot the results. Matplotlib is a python plotting package, actually it enables You to plot in 2D in 3D, prepare some histograms, data visualisations and many many more. Just google matplotlib and You will see many tutorials and examples. But for now let’s concentrate on signals. I am going to put some commands below that enable plotting, but I’m not going to upload any images – if You want to see the results - experiment with matplotlib. import matplotlib.pyplot as plt // first of all we need to import a package from numpy import * // just for further use x1=arange(1,100) n=arange(1,2**10) fsampl=10000 // I specify sampling frequency – feel free to experiment with it ! t=[] dt=1.0/fsampl now let’s generate time vector for i in range(1,2**10): t.append(n[i-1]*dt) A1,A2,A3,A4=1,2,4,8 // our amplitudes f1,f2,f3,f4=100,200,400,800 // and frequencies x1,x2,x3,x4=[],[],[],[] for i in range(1,2**10): x1.append(A1*cos(2*pi*f1*t[i-1])) x2.append(A2*cos(2*pi*f2*t[i-1])) x3.append(A3*cos(2*pi*f3*t[i-1])) x4.append(A4*cos(2*pi*f4*t[i-1])) So it was just a signal generation, now we are going to plot sth. plt.plot(t,x1) // really simple x-values and y-values for x-y plot
plt.show() // we have to of course show it to orselves, if not it will be stored in memory until we show it, modify it or close the program plt.plot(t,x1,’r.-’) // other possibilities plt.show() plt.plot(t,x1,’ro’) plt.show() plt.plot(t,x1,’r.-’,t,x2,’b,-’,t,x3,’g.-’,t,x4,’y.-’) // it plots all the signals in one window, but the general outlook is a mess but this option is useful in case we want to compare two signals, like this : plt.plot(t,x1,’r.-’,t,x2,’g.-’) where b,g,y,r are the colors and .- is how we want our signal to be displayed Ok, so I’ve showed You how to display our signals in a really basic way – believe me it’s enough for most examples, here I want to mention that there is one option of plotting – plt.plot(x) where x is a vector, but sometimes it can be tricky. Now, I’ll show You another example, and one more way of plotting. import matplotlib.pyplot as plt from numpy import * n=arange(1,2**8) fsampl=10000 dt=1.0/fsampl x1,x2,x3=[],[],[] f1,f2,f3=100,200,1000 n1=ones(100) n2=zeros(100) wave=hstack([n1,n2,n1,n2,n1,n2,n1,n2]) // we are generating a wave similar to unipolar t=[] for i in range(1,2**8): t.append(n[i-1]*dt) x1.append(cos(2*pi*f1*t[i-1])) x2.append(cos(2*pi*f2*t[i-1])) x3.append(cos(2*pi*f3*t[i-1]))
plt.subplot(4,1,1) // subplot allows us to display many signals in one window, but seperate plots plt.plot(wave) // check out the subplot syntax in matplotlib docummentation plt.subplot(4,1,2) plt.plot(t,x1) plt.subplot(4,1,3) plt.plot(t,x2) plt.subplot(4,1,4) plt.plot(t,x3) plt.show() Summing up, matplotlib is a really useful package, the best one if it comes to plotting in Python. It’s a really handy tool and if You want to do some scientific research with Python’s help it’s a “must have”. For us, it will be helpful to for example plot fft module or angle, signals after filtering process, function families, transformations and so on. Before You move to the next chapter (yes there will be a second one, since I’ve splitted this article/tutorial into two parts) I suggest You practise a little.
Digital Signal Processing with Python chapter 2 Author: Jacek Nosal 2011 r.
Introduction Welcome to chapter 2 of this article/tutorial. I’m going to present to You, in main opinion the most important and useful functions of scipy.signal package. Also, I will combine processing with plotting to show results in a clearlier way. Both time and frequency domains will be taken into consideration. At the very beginning I’m also pleased to inform that due to some additional free time, I’ll actually prepare three chapters of this tutorial, the third and last part is going to concentrate on DIP – digital image processing.
FFT, Convolution and Signal Representation FFT – Fast Fourier Transform, which definition I’m not going attach here, gives us the information about the frequencies associated with our signals. Shortly, transform one function into another, and this transformation gives us frequency domain representation. What is really great about FFT is the efficiency of the algorithm, compared to DFT (Discrete Fourier Transform) – it’s much more faster and we receive exactly the same results, no matter which transformation we use. So, if the input signal will be defined like this x=2*cos(2*pi*1000*t)+sin(2*pi*700*t) we expect FFT to (after proper calculations) mention frequencies 700 Hz and 1000 Hz somehow. But first, let’s see how can we perform FFT on vectors, and how can we represent information it gives us. At the start, I import appropriate packages and generate several signals : from numpy import * from pylab import * from scipy.signal import * n=arange(1,2**10+1) fsamp=10000 dt=1.0/fsamp
df=1.0/((2**10+1)*dt) t,f=[],[] for i in range(1,2**10+1): t.append(n[i-1]*dt) f.append(n[i-1]*df) x1=[1,2,3,4,-1,-2,-3,-4] n1=[1,2,3,4,5,6,7,8] x2,x3,x4=[],[],[] for i in range(1,2**10+1): x2.append(0.5*sin(2*pi*3000*t[i-1])) x3.append(5*cos(2*pi*2000*t[i-1])) x4.append(2*cos(2*pi*4000*t[i-1])) X1=fft(x1) subplot(3,1,1) stem(n1,x1) subplot(3,1,2) stem(n1,abs(X1) subplot(3,1,3) stem(n1,angle(X1) show()
Image shows us the results only for x1, first plot – input signal, second plot – abs(fft(x1)), third plot – angle(fft(x1)). As FFT gives us complex numbers, there is a trend to plot the modulus and the argument of all components. We see that there are several lines on a second plot, which represents modulus. First conclusion – there is a certain symmetry modulus is an even function and argument odd. However, right now we are not able to tell which frequency is associated with our signal, because of the fact that we are plotting our signal like this index-value, instead of frequency-value. Next example will show how to properly plot a signal in frequency demain. Summing up, it’s pretty easy to calculate fft and present some values with SciPy and Matplotlib. subplot(3,1,1) plot(f,x2,’.-’) subplot(3,1,2) plot(f,abs(fft(x2)),’.-’) subplot(3,1,3) plot(f,angle(fft(x2)),’.-’)
Ok, now we are dealing with x2 – a cos wave, which frequency is equal to 3000 Hz. Because, of this high frequency plot of the input signal is not very clear. In this example, we are plotting the signal in frequency domain, so we can see two peaks – one near 3000 Hz and the second one near 7000 Hz. This is very common problem, and I’m not going to describe it’s background – shortly, first peak represents main frequency let’s call it f0, and the second peak represents fp-f0 = 10000-3000=7000 Hz, actually it’s due to symmetry of fft, the fact that we are dealing with complex numbers, sampling process itself and PERIODICITY of the signal. So, our plot’s centre of symmetry is in 5000 Hz (fsampling/2). It would be nicer if it would be around 0, then we will be just concerned by the part on the right. x=[] for i in range(1,2**10+1):
x.append(x2[i-1]+x3[i-1]+x4[i-1]) # I’m going to perform 1024-FFT and use fftshift like this: w=fftfreq(1024)*fsamp #generation of frequency vector, and also I normalize it using sampling frequency X=fft(x,1024) # 1024-FFT – FFT that computes for us 1024 values X=fftshift(X) # shifting operation w=fftshift(w) # shifting operation subplot(3,1,1) plot(n,x,’r.-’) subplot(3,1,2) plot(w,abs(X)) subplot(3,1,3) plot(w,angle(X)) show()
(Sorry for poor image quality). Now as we added all three signals, we were dealing with 3 frequencies 2000,3000 and 4000 Hz. After performing 1024-FFT, and shifting operations, we are able to tell that those 3 frequencies actually represent our signal in frequency domain. So, we are able to perform spectral analysis of practically any signal (try experimenting with signal + AWGN (additive withe gaussian noise)), and present results in a very clear way. Let’s talk about a convolution for a while. Convolotion is the basic process in case of filtering audio signals as well as images. But engineers don’t actually find this operation “user-friendly”, it’s very confusing and tricky. But ! We have our FFT and one of it’s (actually Fourier Transform’s) properties is that we can replace x1(t)[convolution]x2(t) X1(2pif)*X2(2*pi*f). So convolution in time domain, equals multiplication in frequency domain. There is one point that needs to be told, there are two kinds of convolution – linear convolution and circular convolution. It’s easy to use fft and ifft to calculate
circular convolution, but it becomes harder with the linear one. Below You can see a list of functions dealing with these topics : convolve(x1,x2) #linear convolution of two signals {numpy package } fftconvolve(x1,x2) #where x1,x2 – arrays, fftconvolve can be found in scipy.signal.signaltools x1=[ ... ], x2 =[ ... ], x=fft(x1)*fft(x2) x=ifft(x) #this is how we can perform circular convolution it’s similar with correlation of two signals
Window Functions Usually at this point the author starts with difference equations, that describe LTI systems. But, as the windows topic has not been covered yet I think that it’s a good time to mention window functions. The main concept of window function is to “cut” a certain part of a signal and put an emphasis on it. Windows are often used with FFT, I don’t know if You remember of not but due to fft symmetry and signal periodicty we were observing peaks in places that corresponded following frequencies: f0,fsamp-f0,2fsamp-f0 and so on… .We have coped with the problem using FFT-1024 and fftshift, another idea is to apply a window on a signal that will “cut” desired part. Below You can find a plot that shows most commonly used windows. Shall I mention that they’re also used in spectrogram functions, and filtering itself. To use a window just perform something like that : import scipy.signal as tool x=tool.window_name(window_length) YELLOW – Hamming GREEN – Hann BLUE – Bohmann The length of all windows is 1024, we can see the difference – not a big one, but sometimes the window parameters are really important. Of course, there are a few more window functions e.g boxcar, gaussian, kaiser or triang.
As another example of windowing functions practical use we will modify previous FFT computing (actually plotting and output modifying, cause calculation process remains the same), below You can find set of instructions and certain image (remember that our x2 signal remains unchanged) : import scipy.signal as tool X2=fft(x2) window=tool.hann(512) zeros=zeros(512) fun=hstack([window,zeros]) subplot(2,1,1) plot(f,abs(X2)),’.-’) subplot(2,1,2) plot(fun) X2=X2*fun plot(f,abs(X2) show() Before I apply window function
After windowing process
Note that in the first plot one signal is represented in a frequency domain and a second one just uses ‘index – value’ representation. I had to apply zeros to my window just to obtain the same signal length as our input’s signal, of course I could have modified my window function in completely another way (I may have used another window function actually – feel free to experiment with that).
