cresspahl
science, computing and nice little unicorns
May 7, 2015
RehLivebot
... als Ersatz für den ReLiveBot. Mitschnitte von xenim zum runterladen, inkl. Pre und Postshow. Zum RehLiveBot.
Labels:
download,
mitschnitt,
podcast,
rehlivebot,
relive,
relivebot,
streamrip,
xenim
April 17, 2015
Howto: Access/Convert HAMEG *.mes file from Linux
We have a Hameg HM50142 spectrum analyzer in our lab, which can be controlled via RS232 from a windows box.
For this purpose, the software "AS100E V305 spectrum analyzer software" from www.hameg.com is installed on a Win XP laptop and analyzer is accessed via a USBtoRS232 dongle.
So far, so good. But the stupid thing is that you can not save/export your data in standard formats (e.g. *.dat, *.csv, *.txt), which can be opened with any text editor. Instead, you only can save a *.mes file.
This file turns out to be a Microsoft JET Database (http://en.wikipedia.org/wiki/Microsoft_Jet_Database_Engine), which is proprietary and also flagged as depricated by MS. To access such an database, you could use "MS Access", the database program from Microsoft (office package).
Fortunately, there is a way to dump the data using free software: MDBTOOLS (https://github.com/brianb/mdbtools). You can compile it from github source, but it also should be available via the package managing source of your distro (aptget, pacman, etc).
I have a test measurement called "hameg.mes".
To see which tables are stored inside this file, one can use the tool mdbtables:
PS: the our version of the HM50142 creates databases in the JET3format.
For this purpose, the software "AS100E V305 spectrum analyzer software" from www.hameg.com is installed on a Win XP laptop and analyzer is accessed via a USBtoRS232 dongle.
So far, so good. But the stupid thing is that you can not save/export your data in standard formats (e.g. *.dat, *.csv, *.txt), which can be opened with any text editor. Instead, you only can save a *.mes file.
This file turns out to be a Microsoft JET Database (http://en.wikipedia.org/wiki/Microsoft_Jet_Database_Engine), which is proprietary and also flagged as depricated by MS. To access such an database, you could use "MS Access", the database program from Microsoft (office package).
Fortunately, there is a way to dump the data using free software: MDBTOOLS (https://github.com/brianb/mdbtools). You can compile it from github source, but it also should be available via the package managing source of your distro (aptget, pacman, etc).
I have a test measurement called "hameg.mes".
To see which tables are stored inside this file, one can use the tool mdbtables:
>mdbtables hameg.mes
Measurement Limitlines SettingsNow one can dump one specific table using mdbdump:
> mdbexport hameg.mes Measurement
Frequency,Sample_dBmW,Average_dBmW,Maxhold_dBmW,QPeak_dBmW,memory_dBmW,Ambient_H_dBmW,Ambient_V_dBmW,V_memory_dBmW,H_memory_dBmW,H_Average_dBmW,H_Maxhold_dBmW,H_QPeak_dBmW,V_Average_dBmW,V_Maxhold_dBmW,V_QPeak_dBmW
2.50000000e+01,5.04000015e+01,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,,,,,,,,,,
2.51000004e+01,5.04000015e+01,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,,,,,,,,,,
2.52000008e+01,
....As visible, the dump is in the CSV format which can be easily piped into a csvfile:
> mdbexport hameg.mes Measurement > hameg_measurement.csv
PS: the our version of the HM50142 creates databases in the JET3format.
>mdbver hameg.mesI was able to view/dump/edit this version with mdbtools version 0.7.1 (clone from github)
JET3
June 7, 2014
Numpy / Scipy / Matplotlib  How did I do ... again?
I almost completely switched my computing/plotting tasks from octave to numpy / scipy and matplotlib.
Here's a little list about things i constantly keep forgetting
Matplotlib
 manually set the ticks on colorbars:
from matplotlib import pyplot as plt
from matplotlib.ticker import FixedLocator
...
cbar = plt.colorbar()
locs = [0,10,20,30]
cbar.locator = FixedLocator(locs)
cbar.update_ticks()
 define position and size of the colorbar (place it freely somewhere)
fig = plt.figure(1)
...
is1 = plt.imshow( .... )
cbaxes = fig.add_axes([0.9, 0.05, 0.02, 0.3])
cb = fig.colorbar(is1, cax = cbaxes))
cb.ax.tick_params(labelsize=7)
sources: http://python4mpia.github.io/plotting/advanced.html
 set position of ticks and labels from left to right axis:
ax1 = plt.subplot(111)
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position("right")
 reduce the size of the legend
plt.rcParams['legend.fontsize']=8
plt.rcParams['legend.borderpad']=0.3
plt.rcParams['legend.handlelength']=0.1  get a second yaxis: ax1a = ax1.twinx()
 create a plot on a machine that has no X (running):

import matplotlib # Force matplotlib to not use any Xwindows backend. matplotlib.use('Agg')
 Legend location codes:
more legend thingies: http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.legend
‘best’ 0 ‘upper right’ 1 ‘upper left’ 2 ‘lower left’ 3 ‘lower right’ 4 ‘right’ 5 ‘center left’ 6 ‘center right’ 7 ‘lower center’ 8 ‘upper center’ 9 ‘center’ 10
 3d plot minimal example with mplot3d #from plottingtools import *
from matplotlib import pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize']=(5,4)
x = np.linspace(1.5,1.5,501)
[xx,yy]=np.meshgrid(x,x)
# some function
f = np.exp( (xx**2+yy**2)/0.5**2) \
+ 0.4*np.exp( ((xx1)**2+(yy+0.2)**2)/0.5**2) \
+ 0.8*np.exp( ((xx+1.2)**2+(yy0.5)**2)/0.5**2)
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)#fig.gca(projection='3d')
surf = ax.plot_surface(xx,yy,f,
#cmap=sunrise(),
antialiased=True,
linewidth=1,
rstride=10,
cstride=10,
shade=False)
ax.view_init(75,80)
fig2 = plt.figure()
ax2 = Axes3D(fig2)
ax2.plot_wireframe(xx,yy,f,
rstride=10,
cstride=10)
ax2.view_init(75,80)
Numpy
 Sort an 2d array like a table in a spreadsheet: with respect to some column/row but keeping the respective rows/columns together
import numpy as np
myarr= np.array( [[4,.0004],[3,.003],[2,.02],[1,.1]])
indl=np.lexsort((myarr[:,1],myarr[:,0]))
myarrsort = myarr[indl]
print indl
print myarrsort
source / other methods: http://stackoverflow.com/questions/8153540/sortanumpyarraylikeatable  load CSV files, beeing able to skip rows etc.
import numpy as np
np.genfromtxt("data.csv", skip_header=1,delimiter=',')
source: http://stackoverflow.com/questions/3518778/howtoreadcsvintorecordarrayinnumpy
December 11, 2013
List of useful linux tools whose names i keep forgetting
 arping  ping IP and get the MAC adress
 engauge (aka engaugedigitizer) get data points from images of plots
 gpicview  an image viewer
 pdfposter  a tool to scale pdfs and print (tiled) poster
 ripit most easypeasy cd ripper and mp3 encoding script
 sitecopy a tool to sync folders via ftp
 tr  can be used to translate characters of a text file  or delete them
December 1, 2013
Xenim streamrecorder.
Als Podcasthörer stellte sich mir folgendes firstworldproblem: Es gibt livePodcasts, die interessante Post/Preshows haben, die nachher nicht im Podcast sind. Weiterhin dauert es mir manchmal zu lange bis zur Veröffentlichung.
Ich habe ein wenig mit Python rumgespielt: Herausgekommen ist ein Skript, welches automatisiert die API von xenim abfragt und gegebenfalls mit VLC einen Stream mitschneidet. Das Skript sollte auch mit Streamabbrüchen umgehen können (Wiederaufnahme).
Vielleicht kann jemand was damit anfangen.
https://github.com/xmhk/xstreamrecorder
Das Ding läuft unter Linux (bei mir auf dem Raspberry Pi), braucht Python 2.(6?), und das Pythonpaket 'requests' (Debian, Ubuntu: sudo aptget install pythonrequests, Arch Linux: pacman S python2requests).
Mit dem Programm 'screen' kann das Skript im Hintergrund laufen.
Anpassung auf den Mac sollte relativ einfach machbar sein, Windows eher: mehhh.
Ich habe ein wenig mit Python rumgespielt: Herausgekommen ist ein Skript, welches automatisiert die API von xenim abfragt und gegebenfalls mit VLC einen Stream mitschneidet. Das Skript sollte auch mit Streamabbrüchen umgehen können (Wiederaufnahme).
Vielleicht kann jemand was damit anfangen.
https://github.com/xmhk/xstreamrecorder
Das Ding läuft unter Linux (bei mir auf dem Raspberry Pi), braucht Python 2.(6?), und das Pythonpaket 'requests' (Debian, Ubuntu: sudo aptget install pythonrequests, Arch Linux: pacman S python2requests).
Mit dem Programm 'screen' kann das Skript im Hintergrund laufen.
Anpassung auf den Mac sollte relativ einfach machbar sein, Windows eher: mehhh.
November 2, 2013
Run Docketport 467 Scanner under Linux
The Problem
But it didn't work out of the box under Linux (Ubuntu 12.04).
Using the sane software tools, i was able to see that the scanner at least is detected:
>sanefindscanner
found USB scanner (vendor=0x1dcc [Document Capture Technologies Inc.], product=0x4812 [DocketPORT 467], chip=GL842?) at libusb:001:003
# Your USB scanner was (probably) detected. It may or may not be supported by
# SANE. Try scanimage L and read the backend's manpage.
When i try
scanimage L
i get
No scanners were identified. If you were expecting something different,
check that the scanner is plugged in, turned on and detected by the
sanefindscanner tool (if appropriate). Please read the documentation
which came with this software (README, FAQ, manpages)
basically the same happens when i start xsane (a graphical scanning programm).
What is wrong: it is the permissions! When i tried
sudo xsane
the scanner was found and worked fine.
Because the permission situation, one is not allowed to use all usb devices available.
Solutions
1. lazy bum version
just use "sudo xsane" to access the scanner. When you save the images, maybe you have to adjust the file user and group of scanned images afterwards.2. make a udev rule
the permissions for usb devices can be adjusted writing rules for the udev (device daemon). You can place them in the folder /etc/udev/rules.d.The file name must have the ending .rules and can be for example 90scanner.rules. It can be important that the number at the beginning is higher that the number of the other filenames, because the system reads the files according to the alphanumeric ordering.
Our file contains the following:
#usb scanner readable by me#
SUBSYSTEMS=="usb", ATTRS{idVendor}=="1dcc", ATTRS{idProduct}=="4812", OWNER="THISisYOU"
(you have to change the owner (here THISisYOU) to your own username). Hint: Everything from "Subsystem" up to the end has to be put in one single line.
Now i was able to access and scan without sudo.
UPDATE:
another variation of udev rules is described on the bottom of this page (paragraph "Permission problem"): https://wiki.archlinux.org/index.php/Scanner
August 27, 2012
How to get rid of old google chrome versions on your Mac
Google chrome is a great browser, which automatically keeps itself up to date. For some reasons, it keeps all older versions which eat up disk space. This becomes obvious when your hard drive storage is not *that* big (as on my Macbook Air with a 64 GB SSD).
I wondered why my free disk space constantly shrinks. On linux, baobab is a great tool to visualize the file sizes on you system. I googled for alternatives on mac and found a java program called jdiskreport, which does the job.
In my case, the google chrome app folder was about 2 gigs big, growing since about 1.5 years of usage. It turned out that over 90 percent of that usage was taken by old versions of chrome which aren't useful anymore.
Where are the versions located?
the standard path is
/Applications/Google Chrome.app/Contents/Versions
How to get rid of the older versions?
You can simply delete them.
Just type
cd /Applications/Google Chrome.app/Contents/Versions
into a terminal and use ls to see which versions are installed.
now you should delete all the versions but the actual version (with the highest number), e.g.
sudo rm fr 21.0.1180.79
to delete the version 21.0.1180.79.
The command remove (rm) has to be used with sudo to gain administrative rights.
Another, even simpler possibility ist to just delete chrome from the /Application folder and download and reinstall chrome.
I wondered why my free disk space constantly shrinks. On linux, baobab is a great tool to visualize the file sizes on you system. I googled for alternatives on mac and found a java program called jdiskreport, which does the job.
In my case, the google chrome app folder was about 2 gigs big, growing since about 1.5 years of usage. It turned out that over 90 percent of that usage was taken by old versions of chrome which aren't useful anymore.
Where are the versions located?
the standard path is
/Applications/Google Chrome.app/Contents/Versions
How to get rid of the older versions?
You can simply delete them.
Just type
cd /Applications/Google Chrome.app/Contents/Versions
into a terminal and use ls to see which versions are installed.
now you should delete all the versions but the actual version (with the highest number), e.g.
sudo rm fr 21.0.1180.79
to delete the version 21.0.1180.79.
The command remove (rm) has to be used with sudo to gain administrative rights.
Another, even simpler possibility ist to just delete chrome from the /Application folder and download and reinstall chrome.
Labels:
baobab,
chrome,
delete,
disk space,
jdiskreport,
mac,
old,
usage,
versions
July 23, 2012
octave quickies: change the position of the yaxis from left to right
Sometimes (e.g. when you have multible plots on one page) you want the yticks and ylabeling on the right side of the plot (default is left).
This can be done by changing the current axis properties:
set ( gca(), "yaxislocation","right")
This can be done by changing the current axis properties:
set ( gca(), "yaxislocation","right")
May 20, 2012
high quality plots in octave
When you're about to publish in a scientific journal, one problem appears: how do I create high quality plots for my paper?
You can use all kinds of expensive programs (Origin, for example), but free software is as well capable of doing this. For print, you definitely want to use some vector based image formats, of which EPS (encapsulated postscript) is the most common.
Well let's look at a simple eps, created with octave:
This looks quite poor:
A good idea is to look whether the journal you're intending to publish in has any style guidelines. Often, detailed information about the size, font, and axis labels can be found there.
A list of task that can be done to improve the above plot:
I hope this is a good starting point. The above pictures where created on an Ubuntu 10.04 machine with octave 3.6.x. Now the code:
You can use all kinds of expensive programs (Origin, for example), but free software is as well capable of doing this. For print, you definitely want to use some vector based image formats, of which EPS (encapsulated postscript) is the most common.
Well let's look at a simple eps, created with octave:
generic EPS plot. we have to work on this. code for this example is below 
 the proportions of data points, axis, ticks are not good
 the fonts look messy
 there is no color
A good idea is to look whether the journal you're intending to publish in has any style guidelines. Often, detailed information about the size, font, and axis labels can be found there.
A list of task that can be done to improve the above plot:
 we'll use colors
 we'll use another font: a sansserif is standard in scientific plots (e.g. Helvetica or Arial). The font size will be changed.
 we'll change the proportions of the image. The data points and lines are the most important issue on our plot, so we'll make them quite big and fat
 the outer proportions of the image are changed. In a twocolumn layout which is common in many publications, often a widthtoheight ratio of 2:1 is demanded
 we'll make the axis a little thicker. In some cases, it is desirable to have the ticks outside
 Often, the axis are heavily overloaded, there is too much information distracting from the plot content. This can be reduced by leaving some ticklabels out.
The result looks more nicely:
an improved version. 
I hope this is a good starting point. The above pictures where created on an Ubuntu 10.04 machine with octave 3.6.x. Now the code:
clear all;
more off;
%
% create some data: polynoms with random noise
%
points = 15
x = linspace(5,5,points);
y1 = 0.3*x.^2 + 5 *(rand(1,points).5)+15;
y2 = 0.1*x.^3 + 5 *(rand(1,points).5);
%
% fit some curves to the noisy data
%
[p1,s1] = polyfit(x,y1,3);
[p2,s2] = polyfit(x,y2,3);
tx = linspace(min(x),max(x),500)
%
% create a generic plot
%
figure(1)
plot( x,y1,'x',
x,y2,'o',
tx,polyval(p1,tx),
tx,polyval(p2,tx))
xlabel 'X'
ylabel 'Y'
legend('y1',
'y2',
'fit1',
'fit2')
print('generic.eps')
%
% create a more advanced version of the plot
%
figure(2)
%
% in the plot command
%  define markersize
%  define linewidth
%  assign colors
%
plot( x,y1,'x'
,'linewidth',6
,'color',[1,0,0]
,'markersize',6
,x,y2,'o'
,'color',[0,0,1]
,'markerfacecolor',[0,0,1]
,'markersize',6,
tx,polyval(p1,tx),
'linewidth',3,
'color',[1,0,0],
tx,polyval(p2,tx),
'linewidth',3,
'color',[0,0,1])
%
% now let's adjust the axis:
%  make the line thicker
%  make the ticks outside
%  adjust ticklength
%  use custom ticks
%  apply custom ticklabels
%
set(gca(),
'linewidth',2,
'tickdir','out',
'ticklength',[0.005,0.005],
'xtick',[6:6],
'xticklabel',{'','5','','','','','0','','','','','5',''},
'ytick',[20:10:30],
'ylim',[20,30],
'yticklabel',{'','20','','0','','20',''})
xlabel 'X'
ylabel 'Y'
%
% the legend shall appera on the bottom right of the plot
%
legend('y1','y2','fit1','fit2','location','southeast')
% let's print this thing to a eps.
% we use a custom size of 900x450,
% a defined Font and Size
% and the depsc2 driver (EPS 2 color) driver
print('test2.eps','S900,450',
'depsc2',
'F:Helvetica:12',
'tight'
)
I guess this is a good starting point for creating printable plots that look good. When you're not sure what can be done: Have a look at the attributes of your graphic and axis objects:
get( gcf()) %get current figure
get( gca()) %get current axis
you can set it with the set command
set( gca(), 'property','value')
Have a look at the octave manual. The mathworks website can also be useful, sometimes ;)
Labels:
axis,
color,
custom,
customize,
depsc2,
eps,
figure,
font,
high quality,
labels,
legend,
linewidth,
markerfacecolor,
markersize,
matlab,
octave,
origin,
plot,
publication,
tick
April 21, 2012
octave quickies: customizing the colorbar
When you create colormap or 3Dplots, you may wish to customize the axis of your plot (e.g. change the linewidth, etc). One more difficult task is to change the colorbar of the plot.
As in matlab, in octave things like lines, plot points, axis and so on are basically graphical objects (see the octave manual). To change the properties of some object, you have to know its handle (although for some, e.g. xlim, there are special function to do so).
One thing I stumble across regularly is the colorbar, where i want to change the linewidth and the ticks. The octave manual is not very verbose at this point, but you basically can use the matlab doc.
Example: we create a blank plot with a colorbar. Then we aquire the colorbar handle and play around with the colorbar properties:
figure(1)
caxis([0,50]) %change the range
colorbar()
colormap( summer(250)) %use some different map than jet
%get the colorbar handle
cbh = findobj( gcf(), 'tag', 'colorbar')
set( cbh, %now change some colorbar properties
'linewidth', 2,
'tickdir', 'out',
'ticklength',[0.005,0.005],
'ytick', [0, 10, 23, 42, 50],
'yticklabel',{'zero','ten','23','42','fifty'})
As for all objects, you get get a complete list of the actual colorbars properties with the get function:
get( cbh)
Example : you can change the dimensions (position, height, width) of the colorbar using the position property:
%get the standard value
octave3.4.2:8> get(cbh,'position')
ans = 0.827500 0.110000 0.046500 0.815000
%set new value
set(cbh,'position',[0.8275 0.11, 0.02325 0.81500])
%(makes a slicker colorbar)
you should customize the colorbar to your needs! 
As in matlab, in octave things like lines, plot points, axis and so on are basically graphical objects (see the octave manual). To change the properties of some object, you have to know its handle (although for some, e.g. xlim, there are special function to do so).
One thing I stumble across regularly is the colorbar, where i want to change the linewidth and the ticks. The octave manual is not very verbose at this point, but you basically can use the matlab doc.
Example: we create a blank plot with a colorbar. Then we aquire the colorbar handle and play around with the colorbar properties:
figure(1)
caxis([0,50]) %change the range
colorbar()
colormap( summer(250)) %use some different map than jet
%get the colorbar handle
cbh = findobj( gcf(), 'tag', 'colorbar')
set( cbh, %now change some colorbar properties
'linewidth', 2,
'tickdir', 'out',
'ticklength',[0.005,0.005],
'ytick', [0, 10, 23, 42, 50],
'yticklabel',{'zero','ten','23','42','fifty'})
As for all objects, you get get a complete list of the actual colorbars properties with the get function:
get( cbh)
Example : you can change the dimensions (position, height, width) of the colorbar using the position property:
%get the standard value
octave3.4.2:8> get(cbh,'position')
ans = 0.827500 0.110000 0.046500 0.815000
%set new value
set(cbh,'position',[0.8275 0.11, 0.02325 0.81500])
%(makes a slicker colorbar)
March 29, 2012
simple algorithm for 2DPeakfinding
When it comes to analyzing data, a common task is to identify peaks in an xy dataset. There are several ways to do that. Maybe the first idea that would come into one's mind is to look for zeros in the derivatives. This works fine for smooth theoretical curves but fails when noise is present (as it is certainly all kinds of experimental data).
To get rid of noise you can apply a lowpass filter on your data, but in some cases this is not desired. A simple approach is to use an algorithm that can does
I've implemented such piece of code for octave that does the job. The code below created the picture shown. The function expects three input parameters: the environment for maxima, the data list and a threshold value.
function peaklist = SimplePeakFind ( environment, data, thresh)
%
% 
%
% function peaklist = SimplePeakFind ( environment, data, thresh)
%
% finds the positions of peaks in a
% given data list. valid peaks
% are searched to be greater than
% the threshold value.
% peaks are searched to be the maximum in a certain environment
% of values in the list.
%
% 
%
listlength = length( data );
peaklist = zeros(listlength,1); % create blank output
SearchEnvHalf = max(1,floor( environment/2));
%
% we only have to consider date above the threshold
%
dataAboveThreshIndx = find (data >= thresh);
for CandidateIndx = 1 : length(dataAboveThreshIndx)
Indx = dataAboveThreshIndx(CandidateIndx);
%
% consider list boundaries
%
minindx = Indx  SearchEnvHalf;
maxindx = Indx + SearchEnvHalf;
if (minindx < 1)
minindx = 1;
end
if ( maxindx>listlength)
maxindx = listlength;
end
%
% data( CandidateIndx ) == maximum in environment?
%
if (data(Indx) == max( data( minindx : maxindx)))
peaklist(Indx) = Indx;
end
end
% finally shrink list to nonzerovalues
peaklist = peaklist( find( peaklist));
end
In many cases, this works fine for me. Let's look at a simple test case: we numerically create a number of gaussian pulses with random widths, positions and heights. Then we add noise.
clear all;
more off;
N = 1500;
T = linspace(20,20,N); %xvector
% create some random gaussian peaks
Npeaks = 9
ppos = T(round(rand(1,Npeaks)*N))
pheight = rand(1,15)*1;
pdur = rand(1,15)*1;
data = zeros(1,N);
for i = 1:Npeaks
data = data + pheight(i) * exp( ((Tppos(i))/pdur(i)).^2);
endfor
% apply noisenoise = rand(1,N)*0.1;
datanoise = data+noise;
% Find peaks that are maxima in an
% environment of 25 data points and that
% bigger than 0.15 (noise level)
peaklist = SimplePeakFind(25, datanoise, 0.2);
%Plot result
figure(1)
plot(T,datanoise,"linewidth",2,
T(peaklist), datanoise(peaklist),
"o","markersize",5,"markerfacecolor",[1,0,0])
legend("noisy data","found peaks")
xlabel "x"
ylabel "y"
set(gca(),"tickdir","out","linewidth",2,"ticklength",[0.005,0.005])
To get rid of noise you can apply a lowpass filter on your data, but in some cases this is not desired. A simple approach is to use an algorithm that can does
 find peak candidates that are above some userdefined threshold level
 find out for each candidate whether it is the local maximum of a user defined environment
 return the valid peaks
the algorithm is able to find peaks even in noisy data 
I've implemented such piece of code for octave that does the job. The code below created the picture shown. The function expects three input parameters: the environment for maxima, the data list and a threshold value.
function peaklist = SimplePeakFind ( environment, data, thresh)
%
% 
%
% function peaklist = SimplePeakFind ( environment, data, thresh)
%
% finds the positions of peaks in a
% given data list. valid peaks
% are searched to be greater than
% the threshold value.
% peaks are searched to be the maximum in a certain environment
% of values in the list.
%
% 
%
listlength = length( data );
peaklist = zeros(listlength,1); % create blank output
SearchEnvHalf = max(1,floor( environment/2));
%
% we only have to consider date above the threshold
%
dataAboveThreshIndx = find (data >= thresh);
for CandidateIndx = 1 : length(dataAboveThreshIndx)
Indx = dataAboveThreshIndx(CandidateIndx);
%
% consider list boundaries
%
minindx = Indx  SearchEnvHalf;
maxindx = Indx + SearchEnvHalf;
if (minindx < 1)
minindx = 1;
end
if ( maxindx>listlength)
maxindx = listlength;
end
%
% data( CandidateIndx ) == maximum in environment?
%
if (data(Indx) == max( data( minindx : maxindx)))
peaklist(Indx) = Indx;
end
end
% finally shrink list to nonzerovalues
peaklist = peaklist( find( peaklist));
end
In many cases, this works fine for me. Let's look at a simple test case: we numerically create a number of gaussian pulses with random widths, positions and heights. Then we add noise.
clear all;
more off;
N = 1500;
T = linspace(20,20,N); %xvector
% create some random gaussian peaks
Npeaks = 9
ppos = T(round(rand(1,Npeaks)*N))
pheight = rand(1,15)*1;
pdur = rand(1,15)*1;
data = zeros(1,N);
for i = 1:Npeaks
data = data + pheight(i) * exp( ((Tppos(i))/pdur(i)).^2);
endfor
% apply noisenoise = rand(1,N)*0.1;
datanoise = data+noise;
% Find peaks that are maxima in an
% environment of 25 data points and that
% bigger than 0.15 (noise level)
peaklist = SimplePeakFind(25, datanoise, 0.2);
%Plot result
figure(1)
plot(T,datanoise,"linewidth",2,
T(peaklist), datanoise(peaklist),
"o","markersize",5,"markerfacecolor",[1,0,0])
legend("noisy data","found peaks")
xlabel "x"
ylabel "y"
set(gca(),"tickdir","out","linewidth",2,"ticklength",[0.005,0.005])
The result is shown in the picture above: all 9 peaks are identified correctly. Depending on your data, you'll have to play around with the parameters environment and thresh. Having the the peak positions, further data analysis (like curve fitting etc) is much easier.
Links: of course people already thought of this problem:
Links: of course people already thought of this problem:
Labels:
algorithm,
fitting,
matlab,
noise,
noisy data,
octave,
peak finding,
peak identifying algorithm
March 26, 2012
simulating nonlinear oscillators using the julia language
In my post 'simulating nonlinear oscillators using octave' i showed you how to ...  i guess you can read. Today i discovered, that a package similar to the odepkg also is included in julia. You just have to find it.
Basically it works quite similar, with some little tweaks.
There is a file called 'ode.jl' in the 'extras' folder located in the julia source folder. We just have to load this thing and have some ode solvers. The drawback is that up to now there the only way (i found) to adapt the accuracy is to change the lines
rtol = 1.e5
atol = 1.e8
in the function ode23
and
tol = 1.0e7
in the definition of oderkf{T}, which is a quite dirty method.
Another issue is that the only arguments of the ode solvers are up to now
(F::Function, tspan::AbstractVector, y_0::AbstractVector)
which means we have to hide additional parameters in our solution vector (and make its derivatives zero).
A code example, that works for me (asymmetric nonlinear potential)
# add the path to your juliasource/extras folder
load("/Users/user/juliadir/extras/ode.jl")
function rhs_asym(t, x)
#
# position and velocity are stored in x[1] and x[2]
#
xp = x[1]
vx = x[2]
#
# the parameters are stored in x[3] and x[4]
#
k=x[3]
alpha =x[4]
dxp = vx
dv =  (k * xp  alpha * xp^2)
dx = [dxp; dv;0;0] # the two zeros ensure that k and alpha
# will not change
end
# problems can occur when you forget to give these values
# explicitly as Floats ...
k=1.
alpha =0.2
tend=42.
initialDisplacement = 0.0
initialVelocity = 1.
xini = [initialDisplacement; initialVelocity; k; alpha]
#
# to get higher accuracy, change the values of rtol atol or tol
# in the ode.jlfile (or better in a copy of it).
#
# now let's integrate
(T, X) = ode45( rhs_asym, [0; tend], xini)
#
# we store the whole data in one matrix and save it,
# so we can load and plot with e.g. octave
#
M=[T1 Xlin]
csvwrite("oscil.csv",M)
Basically it works quite similar, with some little tweaks.
There is a file called 'ode.jl' in the 'extras' folder located in the julia source folder. We just have to load this thing and have some ode solvers. The drawback is that up to now there the only way (i found) to adapt the accuracy is to change the lines
rtol = 1.e5
atol = 1.e8
in the function ode23
and
tol = 1.0e7
in the definition of oderkf{T}, which is a quite dirty method.
Another issue is that the only arguments of the ode solvers are up to now
(F::Function, tspan::AbstractVector, y_0::AbstractVector)
which means we have to hide additional parameters in our solution vector (and make its derivatives zero).
A code example, that works for me (asymmetric nonlinear potential)
# add the path to your juliasource/extras folder
load("/Users/user/juliadir/extras/ode.jl")
function rhs_asym(t, x)
#
# position and velocity are stored in x[1] and x[2]
#
xp = x[1]
vx = x[2]
#
# the parameters are stored in x[3] and x[4]
#
k=x[3]
alpha =x[4]
dxp = vx
dv =  (k * xp  alpha * xp^2)
dx = [dxp; dv;0;0] # the two zeros ensure that k and alpha
# will not change
end
# problems can occur when you forget to give these values
# explicitly as Floats ...
k=1.
alpha =0.2
tend=42.
initialDisplacement = 0.0
initialVelocity = 1.
xini = [initialDisplacement; initialVelocity; k; alpha]
#
# to get higher accuracy, change the values of rtol atol or tol
# in the ode.jlfile (or better in a copy of it).
#
# now let's integrate
(T, X) = ode45( rhs_asym, [0; tend], xini)
#
# we store the whole data in one matrix and save it,
# so we can load and plot with e.g. octave
#
M=[T1 Xlin]
csvwrite("oscil.csv",M)
Labels:
code,
example,
julia,
matlab,
nonlinear oscillator,
octave,
ode23,
ode45,
simulation
March 22, 2012
Simulating nonlinear oscillators using octave
Hi there,
numerics are usefull especially when nonlinear problems are considered. Here, I want to show how to calculate the trajectory for some simple onedimensional oscillators.
You can imagine such an oscillator as a little mass attached to a spring. There is a point of rest, where all forces vanish. When the oscillator is displaced, a force appears originating from the spring.
When energy is put in the system, the mass oscillates. Depending on the form of the potential (which basically is the integral of the force with respect to the displacement) we can distinguish harmonic and anharmonic oscillators.
From textbooks (or wikipedia) we know the following items, which we want to check using numerics:
 for the harmonic (linear) oscillator, the trajectory is a pure sine wave of some frequency
 for the anharmonic (nonlinear) oscillator, more frequency components appear:
 there are only odd components for the symmetric nonlinear potential
 there are odd and even components for the asymmetric nonlinear potential
1. The linear (harmonic) oscillator
This is a well known textbook example. The potential energy of the harmonic oscillator is simply a parabula:
Epot = 1/2 k x^2
The resulting force is the negative derivative of the potential:
F =  dE/ds =  k x
This is Hooke's law which is a simple model for springs.
It is also the reason why the harmonic oscillator is called the linear oscillator: the force is directly proportional to the displacement (linear dependency).
Epot = 1/2 k x^2
harmonic potential. The potential energy grows quadratically with the displacement. 
F =  dE/ds =  k x
This is Hooke's law which is a simple model for springs.
Force acting in the harmonic potential 
The equation of motion is can be derived from the force
d^2 x /dt^2 = F/m = k/m x
where x is the displacement and k the spring constant. We can easily give the solution for the equation:
x(t) = A1 cos( k^2/m^2 t + A2)
where A1 and A2 are two constants that have to be adapted to the initial values of the problem. The numerical solution using ODEs can be done following the steps described in one earlier post. As the linear oscillator is a special case of the more general problem i'll give the code later in this text.
After integration, we get a trajectory looking as the following (k=1, initial position =0, initalspeed = 2)
This looks like a nice, pure sine wave (as expected). To get more insight, we use the fourier fransform
trajectory of the harmonic oscillator 
of the trajectory x(t). You have to be careful using results from octaves ode solvers, because in most cases the output points are not equally spaced in time. Using the interpolation function interp1 can help to linearize the data.
The fourier spectrum looks like the following:
As expected, we only see one single frequency component, which means we indeed have a harmonic oscillation here.
The fourier spectrum looks like the following:
fourier components for the harmonic oscillator: there is only one frequency present 
2. The anharmonic / nonlinear oscillator
Considering the example with the spring, Hooke's law is valid for small displacements. When the displacement gets bigger, additional terms can appear. We want to consider two cases: a symmetric and an antisymmetric nonlinear potential. We give the potential energy:
Epot,sym = 1/2 k x^2 + 1/3 alpha x^3
Epot,asym = 1/2 k x^2 + 1/3 alpha x^3
here alpha is a nonlinearity constant. The curves of the potential are shown below:
Again, the forces occuring are the derivatives of the potential with respect to the displacement:
Considering the example with the spring, Hooke's law is valid for small displacements. When the displacement gets bigger, additional terms can appear. We want to consider two cases: a symmetric and an antisymmetric nonlinear potential. We give the potential energy:
Epot,sym = 1/2 k x^2 + 1/3 alpha x^3
Epot,asym = 1/2 k x^2 + 1/3 alpha x^3
here alpha is a nonlinearity constant. The curves of the potential are shown below:
Two aharmonic potentials, compared with the harmonic one. They differ in their symmetry 
Fsym =  dEpot,sym/ds =  (k x + alpha x x)
Fasym=  dEpot,asym/ds =  (k x + alpha x^2)
We see additional dependencies between the force and the displacement which are of the order of x^2. This is why these oscillators are called nonlinear oscillators. The difference to the linear case becomes clear when we look at the force plot:
Forces derived from the potentials above. Note that the asymmetric force is equal to the symmetric force for positive displacements, but differs for negative displacements. 
For small displacements, the forces are almost identical. They begin to differ considerably, when the displacement grows. For the symmetric potential, the absolute value of the force is equal for x and x. For the asymmetric potential, the force for negative displacements is smaller than for positive values.
The equations of movements can be written as
(symmetric)
d^2 x /dt^2 = F/m =  (k/m x + alpha/m x x)
(asymmetric)
d^2 x /dt^2 = F/m =  (k/m x + alpha/m x^2)
Again, we can easily solve this with the odepkg (see code sample below). I used the same input energy as for the harmonic example resulting in the following trajectories:
The two new trajectories basically also show sine osciallations. As the force for the symmetric potential is stronger for big displacements, the maximum amplitude is smaller. This also results in a higher frequency for the oscillation. For the asymmetric potential, the oscillation is shifted slightly towards negative values in the mean.
We can get the fourier components of this oscillations when analyzing the time traces x(t) of the signals.
We have a closer look on the symmetric harmonic potential:
Besides the fact, that the base frequency is higher than for the harmonic oscillator, we see an additional frequency component! Its value is three times the one of the base frequency (third harmonic). Also, energy is converted into the fifth harmonic.
Trajectories for the nonlinear oscilators (two kinds), compared with the one of the harmonic oscillator 
The two new trajectories basically also show sine osciallations. As the force for the symmetric potential is stronger for big displacements, the maximum amplitude is smaller. This also results in a higher frequency for the oscillation. For the asymmetric potential, the oscillation is shifted slightly towards negative values in the mean.
We can get the fourier components of this oscillations when analyzing the time traces x(t) of the signals.
We have a closer look on the symmetric harmonic potential:
spectrum for the symmetric nonlinear oscillator. Only odd multibles of the base frequency are generated 
A similar behavior can be seen for the asymmetric potential:
Here we have generated additional frequencies at two, three and four times the base frequency. The conversion efficiency is higher than for the symmetric potential.
spectrum for the asymmetric nonlinear oscillator. Here, even and odd multibles of the base frequency are generated 
As a result, the simulations show the expected behavior. We can see that nonlinear oscillators are capable of generating multiples of the base frequency. This is extensively used in nonlinear optics, where light is sent through nonlinear crystals. The crystal structure is responsible for the potential, in most cases they have asymmetric potentials.
When e.g. light from a laser (at a wavelength of 1064 nm) of high intensity passes such a crystal, green light (532 nm) can be generated.
see
for further details.
3. Code to solve the linear and nonlinear oscillators using octave
clear all;
more off;
%
% right hand side of the ODE for the symmetric harmonic potential
%
function dx = rhs_sym(t, x, k, alpha,tend)
t/tend %just prints the progress
xp = x(1); % position
vx = x(2); % velocity
dxp = vx; % position change
dv = (k * xp + alpha * xp * abs(xp) ); %velocity change
dx = [dxp, dv]; %give the derivative
endfunction
%
% right hand side of the ODE for the asymmetric harmonic potential
%
function dx = rhs_asym(t, x, k, alpha,tend)
t/tend
xp = x(1);
vx = x(2);
dxp = vx;
dv = (k * xp + alpha * xp^2);
dx = [dxp, dv];
endfunction
% program parameters
k =1;
alpha =0.2;
initialDisplacement = 0;
initialVelocity = 2;
xini = [initialDisplacement, initialVelocity];
tend = 15; %integrate up to that time
% set integration options
options = odeset( 'InitialStep',tend/1e3,
'MaxStep',tend/1e3,
'RelTol',1e6,
'AbsTol',1e6,
'NormControl','on')
% linear case, using alpha = 0
[T1, Xlin] = ode45( @rhs_sym, [0, tend], xini, options, k, 0, tend);
%NL case, symetric potential:
[T2, XnlS] = ode45( @rhs_sym, [0, tend], xini, options, k, alpha,\
tend);
%NL case, symetric potential:
[T3, XnlA] = ode45( @rhs_asym, [0, tend], xini, options, k, alpha,\ tend);
% plot the result
plot (T1, Xlin(:,1),'color',[0.5,0.5,0.5],
T2, XnlS(:,1),'.','color',[0,0,1],
T3, XnlA(:,1),'.','color',[1,0,0]);
xlabel 'time'
ylabel 'displacement'
Hint: for the analysis of the frequency components shown above, the integration time has to be set to a higher value and the odeoptions have to be changed to higher accuracy...
%
% right hand side of the ODE for the symmetric harmonic potential
%
function dx = rhs_sym(t, x, k, alpha,tend)
t/tend %just prints the progress
xp = x(1); % position
vx = x(2); % velocity
dxp = vx; % position change
dv = (k * xp + alpha * xp * abs(xp) ); %velocity change
dx = [dxp, dv]; %give the derivative
endfunction
%
% right hand side of the ODE for the asymmetric harmonic potential
%
function dx = rhs_asym(t, x, k, alpha,tend)
t/tend
xp = x(1);
vx = x(2);
dxp = vx;
dv = (k * xp + alpha * xp^2);
dx = [dxp, dv];
endfunction
% program parameters
k =1;
alpha =0.2;
initialDisplacement = 0;
initialVelocity = 2;
xini = [initialDisplacement, initialVelocity];
tend = 15; %integrate up to that time
% set integration options
options = odeset( 'InitialStep',tend/1e3,
'MaxStep',tend/1e3,
'RelTol',1e6,
'AbsTol',1e6,
'NormControl','on')
% linear case, using alpha = 0
[T1, Xlin] = ode45( @rhs_sym, [0, tend], xini, options, k, 0, tend);
%NL case, symetric potential:
[T2, XnlS] = ode45( @rhs_sym, [0, tend], xini, options, k, alpha,\
tend);
%NL case, symetric potential:
[T3, XnlA] = ode45( @rhs_asym, [0, tend], xini, options, k, alpha,\ tend);
% plot the result
plot (T1, Xlin(:,1),'color',[0.5,0.5,0.5],
T2, XnlS(:,1),'.','color',[0,0,1],
T3, XnlA(:,1),'.','color',[1,0,0]);
xlabel 'time'
ylabel 'displacement'
Hint: for the analysis of the frequency components shown above, the integration time has to be set to a higher value and the odeoptions have to be changed to higher accuracy...
March 20, 2012
How to solve ODEs using octave.
The odepkg package for octave is a great tool to solve ordinary differential equations (ODEs). To solve some problem we only have to bring it into a form octave understands. We can do that following four steps:
 determine the variables of the problem
 figure out which variables change and how they do it
 if higher order derivatives do appear, reduce the order of the system
 write a function that can gives the first derivatives for all values considered and solve it with octave's ode solvers
First example: Shooting a paperball using a slingshot (no friction).
problems like the trajectory of a sphere fired by a slingshot can be solved using ODEs 
This is a very basic example. In fact, we can give the analytical solution (parabula). The paperball is accelerated by the catapult and flies some distance until it hits the groud.
1. The variables of the problem are
 the position of the paperball in space r(t) = ( x(t), y(t) )
 the velocity v(t) = ( vx(t), vy(t) )
 the acceleration a(t) = (ax(t), ay(t) )
2. Now we have to consider how our variables change:
 the change in position is proportional to the velocity at a certain moment: dr/dt = (vx, vy)
 the velocity is changed by the acceleration: dv/dt = (ax, ay)
 for our case, the acceleration is constant. We consider the gravitational force: F=m a where a = g = 9.81 m/s^2 is the gravitational acceleration. It only acts on the ycomponent of the velocity (perpendicular to the ground)
The result should be of the form result(t) = [rx, ry, vx, vy]
What we now have to do is to give octave a vector containing the derivative of our result vector. We do this by defining a function. [drx, dry, dvx, dvy] = deriv_result( time,position)
3. if we did not already do so, we have to formulate the problem as a system of first order ordinary differential equations.
What does this mean? Instead of writing the second derivative of some value, we define some additional variable that is the first derivative of the initial value and consider its changes.
example:
 the acceleration is the second derivative of the position d^2r / dt^2 = a
 we introduce the velocity. it is the first derivative of the position dr/dt = v
 the change of the velocity is given by the acceleration dv/dt = a
 now we have to solve the problem both for the position and the velocity and such the second order ODE has been reduced to two first order ODE
4. solving the problem in octave
We define derivative function
function dr = dr_gravi(t,r)
g = 9.81; %m/s^2
vx = r(3);
vy = r(4);
ax = 0;
ay = g;
dr = [vx,vy,ax,ay];
endfunction
In the program we have to set some initial value. For our problem this is the initial position r0=(x0,y0) and the initial speed v0 = (vx0,vy0)
X0 = 0
Y0 = 0
VX0 = 2
VY0 = 3
initialVector = [ X0, Y0, VX0, VY0]
Also, some integration window is necessary: for our problem, this is the starting time and the ending time for time integration:
StartT= 0 %s
StopT = 0.7 %s
then we let octave do the job:
[T,Result] = ode45( @dr_gravi,[StartT, StopT], initialVector)
Here we used the ode45 solver, but there are many more. See the odepkg documentation for details.
Our result looks like the following
T =
0.00000
0.07000
0.14000
0.21000
0.28000
0.35000
0.42000
0.49000
0.56000
0.63000
0.70000
0.70000
Result =
0.00000 0.00000 2.00000 3.00000
0.14000 0.18597 2.00000 2.31330
0.28000 0.32386 2.00000 1.62660
0.42000 0.41369 2.00000 0.93990
0.56000 0.45545 2.00000 0.25320
0.70000 0.44914 2.00000 0.43350
0.84000 0.39476 2.00000 1.12020
0.98000 0.29231 2.00000 1.80690
1.12000 0.14179 2.00000 2.49360
1.26000 0.05679 2.00000 3.18030
1.40000 0.30345 2.00000 3.86700
1.40000 0.30345 2.00000 3.86700
T is the vector that gives us the time in seconds for some row. Results contains the position and velocity components for the times stored in T. Using this data, we can plot the trajectory and other dynamic information.
The following code plots the trajectory of the ball and the time dependence of the velocity components:
The following code plots the trajectory of the ball and the time dependence of the velocity components:
figure(1)
plot(Result(:,1),Result(:,2),'o');
xlabel 'x position / m'
ylabel 'y position / m'
figure(2)
plot(T, Result(:,3),'.',
T, Result(:,4),'x');
xlabel 'time / s'
ylabel 'velocity / (m /s)'
legend('vx','vy')
legend('vx','vy')
Is there a way to influence the calculation? Of course there is. We can do that using odeoptions
options = odeset( 'RelTol',1e4,
'AbsTol',1e4,
'InitialStep',StopT/1e3,
'MaxStep',StopT/1e3)
[T,Result] = ode45( @dr_gravi,[StartT, StopT], initialVector , options )
There are many options to play around with, of which most are selfexplanatory . You have to find a tradeoff between accuracy and calculation time. Again, see the odepkg documentation.
A more complex example: slingshot with friction.
Additional to the gravity, let's consider the drag force influencing the flight of our paperball. It is a force proportional to the velocity v of an object.
F = 1/2 rho Cd A v^2
Additional parameters appear, which we set here:
 rho, the density of air ( rho = 1.2 kg / m^3)
 Cd  the drag coefficient . We consider a sphere here, which has Cd = 0.45
 A is the crosssection of the sphere
Again, we can split the force into two components: Fx and Fy. To calculate the acceleration components, we also need the mass of our object, as F=ma which leads to a = F/m.
We can give the derivative function as follows with additional parameters (area and mass of the sphere):
function dr = dr_gravi_friction(t,r,area,mass)
g = 9.81; %m / s^2
cdSphere = 0.45;
rhoAir = 1.20; %kg / m^3
frictioncoefficient = 1/2 * rhoAir * cdSphere * area / mass;
vx = r(3);
vy = r(4);
ax =  frictioncoefficient * vx^2; % only friction
ay = ( g + frictioncoefficient * vy^2 ) ; % friction and
% gravitation
dr = [vx,vy,ax,ay];
endfunction
Now we can try out how friction affects the trajectory. We calculate this for a relative heavy sphere, a light sphere and without friction. The radius of the spheres shall be r=2 cm
mass1 = 0.1 %kg
mass2 = 0.001 %kg
area1 = pi * 0.02^2 %m^2
area2 = area1;
options = odeset( 'RelTol',1e4,
'AbsTol',1e4,
'InitialStep',StopT/1e3,
'MaxStep',StopT/1e3)
% trajectory of heavy sphere
[T1,Result1] = ode45( @dr_gravi_friction,[StartT, StopT], \
initialVector , options, area1, mass1 )
% trajectory of light sphere
[T2,Result2] = ode45( @dr_gravi_friction,[StartT, StopT], \
initialVector , options,area2, mass2 )
% trajectory without friction
[T3,Result3] = ode45( @dr_gravi,[StartT, StopT], \
initialVector , options)
When we look at the results, we see the following: The relative impact of the friction force on the light sphere is quite strong. It considerably gets slowed down and hits the ground earlier. The trajectory deviates more and more from the parabula shape while propagating. Contrarily, the heavy sphere almost acts as if friction was absent. This basically is what we expected.
trajectories of spheres 
Subscribe to:
Posts (Atom)