Vincent van Beveren

ITT’s IDL is a tool for data analysis, it consists of a workbench and a programming language. At FOM Rijnhuizen, where I work, it is used a lot of analysis. Unfortunately like almost all proprietary scientific software (Matlab, Origin, Mathematica) it has some issues:

• It has a constricting license:
  • Maximum nr. of users
  • Limited sharing
  • Limited usage of CPU’s
• Its algorithms can not be verified
  • IDL is a black box, how do you know it is mathematically correct?
• Vendor lock-in
  • Once you write your routines in IDL, migration becomes a difficult task
  • Limited 3rd party libraries
  • The vendor can increase the licenses fees to a ridiculous amount without customers leaving
• Its ‘programming language’ does not scale very well
  • Large IDL programs become unreadable
  • Bad constructs force messy code

This basically means that as a scientist you get limited by the usage of such a product, while paying more and more money to get it to work. Also when working somewhere else they might not have the product, and all custom written routines must be rewritten.

As an alternative I am a huge fan of Python, and then specifically all the add-on modules which it provides, like SciPy and Mayavi. These are rich, powerful and often open alternatives, which can be freely used (sometimes only free for non-commecial). However, they are open. You can verify every bit and byte.

So, how do you motivate scientists to migrate to another environment? Once which does not impose limitation, but doesn’t come in a box? Once which requires effort to start with? Well, that is not easy. First of all, you’d never migrate if you would loose all your previous routines. So, for this I started a new project, pyDL, because other packages could not cooperate with SciPy. Currently this codebase is being merged with Micheal McKerns’s pyIDL, becomming pyIDL 1.0. So, this would allow scientists to still use IDL and there routines, while using Python as the glue.

Second you would need to find compelling reason, beyond ‘hey its free and open’. So, what about tools. What if Python can make you more productive. Well, as a language Python certainly scales better. It actually retains its readability after one page of code.  There is a huge library of packages, for example it has a free finite element modeler, sfepy, integrated scientific suites like PythonXY and SAGE. But also a lot of non-science related modules. For example, it is possible to create neat user interfaces using QT, or WX. So, my strategy would be to promote its use.

Free, is actually not such a bad reason. I don’t really know how much license fees are payed, but I know there are a substantial amount. Saving on these fees could create new jobs, allowing more research to be performed.

A long time ago there was this contest, the question was the following:

Suppose as a measure of birth control, the state dictates that a pair can have only one boy. So, if they have a girl, they can try again, but if they have a boy they must stop. Now the question was: Would there eventually be (given that the chance for either a boy or a girl is 50%):

A. more boys

B. more girls

C. about an equal amount of boys and girls

The answer was, according to the contest, that it would not matter, and that the ratio of boys and girls would be 50/50.

After thinking about this for a bit, I can see the point. However, if you think about it a bit more, there are some issues:

- Suppose we do not deal with an infinite number of pairs, but just say, a 100. This would mean there would be a 100 boys after all pairs have tried to get as many children as possible. The chance for a girl however, does not have such a maximum, and could be anywhere between 0 and .. euh, I dunno… a 1000).

- Suppose each family does not want as many children as they can get, but only between 1 and 4. This would mean that in many cases a pair would not produce a boy.This would put the possible nr. of boys between 0 and 100 and the possible nr. of girls between 0 and 400.

So, I wrote a little simulation. It tests a 100 times, a 100 generations with a 100 pairs. Surprisingly almost always, the boys won… So, its not even distributed, but to the contrary to intuition it seems this measure has a bigger chance of producing more boys.

So the real answer is A. more boys. Note though that the bias is very small.

I’m currently in the process of writing a CHIP-8 emulator. If you want to learn how to program an emulator, its a good start, because its not too much work and gives you a good insight into whats involved in writing say a console emulator. You can download the emulator from the page here.

Ofcourse its also fun if you just wish to relive old black-and-white, huge pixel, times and want to play PONG with your friends! The emulator comes loaded with games, like U-BOAT, PONG and BRIX! Real more about CHIP-8 on wikipedia.

Details are on the next page. Have fun!

I have added the game Othello (or Reversi) to this site. Its implemented in JavaScript and you can play against another player or the computer. Its AI is based on a backtracking algorithm, and quite hard to beat (for a rookie like me) :D. It can kick your butt too! Try it on experienced or master

You can find the game here, and a document I wrote a while ago about how to program boardgames opponents here.

In my quest to beat PNG compression I’ve side-tracked a bit. I was working on my old Haartransform-like algorithm and tinkered a bit with it when I came up with a algorithm which could potentially also handles soft-waveforms (haar transform is usually only good with crisp data , like a square waveform). This can be well used for image compression but also for audio. So as an experiment I decided to create an lossy audio codec for it. At the end of this post you can download the codec, and listen to some samples of it. Note that it does NOT beat MP3…

The algorithm is a divide-and-conquer recursive procedure which stores the average and the difference. The idea behind the algorithm is that there are more lower frequencies than higher in a wave-form, and that you should find a efficient way of coding these. The algorith operates somewhat like this:

average = (sample1 + sample2) / 2
difference = sample1 - average

You can now reconstruct the actual samples using average and difference like this:

sample1 = average + difference
sample2 = average - difference

Now this is very similar to a haar transform. When using real values no information is lost in this process.

Now image in a set of 8 samples:

0, 2, 3, 2, 0, -2, -3, -2


Which is a rather crude waveform. When we take for each the average and the difference, we have 4 averages, and 4 differences. We then store the averages in the left-side and the differences in the right side. Effectivly we then seperated the lower frequenies (the averages) from the higher frequencies (the differences).

Applying the function one time the the sequence above will result in something like:

low: 1, 2.5, -1, -2.5, high: -1, +0.5, +1, -0.5

Now in the first 4 samples you might recognize the old wave-form, only its frequencies twice that of the original. The second part are the ‘high’ frquencies. Now we can apply the same transformation to the first part of the of the image, but also to the second part. Which will cause the first part, 4 samples to be split into two halves: 2 low frequences, and 2 high frequences. Applying then the transformation again to the blocks of 2 samples will cause it also to be split into a high-frequency part and a low frequency part. I will show this with only the first block:

1, 2.5, -1, -2.5 => low: 1.75, -1.75 and high: 0,75, 0.75

Then again on the low frequencies:

1.75, -1.75 => low: 0, high: 1.75

The lowest frequency is always the average. Now normally I would apply the algorithm on all blocks in the transformation, not only the lowest. This will cause an interesting decomposition. Reconstructing the wave-form is the exact opposite. Lets get back to the two 4 sample blocks with high and low frequencies:

low: 1, 2.5, -1, -2.5, high: -1, +0.5, +1, -0.5

Imagine I throw away all the high frequencies, make them 0:

low: 1, 2.5, -1, -2.5, high: 0, 0, 0, 0

And reverse the decomposition process, I would get the following wave-form:

1, 1, 2.5, 2.5, -1, -1, -2.5, -2.5

What happend! Your smooth waveform became a square wave! Shown together with the original samples it looks like this (where the bars are the new data)

Now this is bad news for any audio-wave-form. Not only did we remove the old high-frequency component, we also introduced a new high-frequency, that of a square wave. It will make your audio sound like a cheap musical greetings-card!

So the secret to the solutions lies in using interpolation when cacluating the differences. Not only do you take the average of s1 and s2, but also the average of your left neigbour to and your right neigbour. Using the interpolation you use a smooth-upscaling algorithm, removing the annoying square wave-forms.

To upscale smoothly we assume some relation to the average of its neigbour. This means that the left sample is one time its left neigbour average + 3 times its own average divided by 4. This is because the left sample we try to find is 0.5 removed from the mid-point of the nearest average and 1.5 removed from the mind-point of the left average. This means that the nearest average is 3 times as important as the left-average. Since we upscale only by 2 anything beyond 1 sample away is not important.

The following diagram tries to clarify the upscaling algorithm:

Probably not doing a good job, but study the source code okay!

So what do you get when we upscale like this using the previous decomposition:

low: 1, 2.5, -1, -2.5, high: 0, 0, 0, 0

We get the following results

0: 1 (we don't have a left neighbour - normal algorithm uses extrapolation here)
1: (1 * 3 + 2.5) / 4 = 1.38
2: (1 + 2.5 * 3) / 4 = 2.12
3: (2.5 * 3 - 1) / 4 = 1.625
4: (2.5 - 1 * 3) / 4 = -0.125
5: (-1 * 3 - 2.5) / 4 = -1.38
6: (-1 - 2.5 * 3) / 4 = -2.125
7: -2.5 (we don't have the right neigbour - normal algorithm uses extrapolation here)

This resuls in the following diagram:

Its not perfect yet, but its an improvement to the square resize. The actual algorithm also uses some average compensation and extrapolation around the averages. Here is a diagram showing the output after the correction and extrapolation have been added

Note that the recreated wave only uses the average data, meaning its constructed of only half of the information of the original wave.

This is because the average of the two interpolated samples should when averages still be the same as the original nearest neighbour average. Why is that? Well,we are trying to store the sample data in the same number or data elements as the original samples. When the interpolated average is the same as the nearest neigbour the difference between the guessed samples are the real samples is the same, e.g:
|guessedSample1FromAverageA - realSample1| = |guesedSample2FromAverageA - realSample2|
|guessedSample3FromAverageB - realSample3| = |guesedSample4FromAverageB - realSample4|

This is important because when we only need to store |guessedSample1FromAverageA - realSample1|. The table below shows this. Its the data shown in the graph above + the difference from the original

Now it clear to see that the wave can be reconstructed using the average and the difference, and that it would only take up a total of 8 samples, 4 for the lower (averages) and 4 for the higher (differences):

1, 2.5, -1, -2.5 and 0.38, -0.25, -0.63, -0.38.

Now applying the same transformation to both the high and the low components again, until each block to transform is only 2 samples.

The inverse is simply the oposite, reconstructing the weighted averges and adding or substracting the highfrequency components recusivly.

Now we have a reversable transform which also handles smooth data well because of weighted averages instead of nearest neighbour. Now for my image format I am creating a 2D version of this algorithm. For audio data you just need a 1D version.

So now we have some ‘Transform’ thingy. Whats up with that?

Well, as said before this algorithm basically decomposes a set of samples into frequency-amplitudes of wave-forms. The key lies in descarding low-amplitude frequency components, which are unhearable anyway. The important frequencies stay intact, while the unimportant are zeroed. Often its possible to zero the larger part of the frequencies. After that you can push the output of the transform through a statistical compression algorithm like the DEFLATE, and you get pretty good compression.

Here is an example of the audio compressors output (using ogg as media format ):
The original tune:
- original-bigyellow.ogg
The tune compressed at normal setting and decompressed:
- normalcompression-bigyellow.ogg
The tune compressed at high setting and decompressed (sounds aweful):
- highcompression_bigyellow.ogg

Now as you’ll notice the highly compressed audio has more plings and ploings. This is because the subtleties of the sound have gone, only the louder frequencies remain. You can download the original Java source codec here , named ‘Olaf’ (Open Lossly Audio Format).I guess the compression ratio lies around 1/4 for CD quality (note that MP3 is 1/10 for near CD quality).

Though its not very good it could be improved quite a bit. For example overlaying could remove the ~40Hz buzz which comes to play when compressing at higher ratios. Smarter encoding of the remaining frequencies could increase the compression ratio quite a bit.

As far as I know it does not make use of any patented algorithms, and no freaking descrete cosine transform (which everyone seems to love). Olaf is extremly quick to decode and requires no floating point operations

Maybe I’ll turn Olaf into a usable format. For now its just a proof of concept, and contains many bugs and inefficiencies.

I’ve been thinking for a while about lossless photo compression. I believe M$ just released its ‘open standard’ (yeah right…) lossless photo format. There is, in my opinion, not yet such an alternative in the open-source community.

Sure you have PNG, but reading through the PNG documentation I noticed it has a few weaknesses, which make it ideal for illustrations and web-images, but less ideal for photos. PNG is compressed in two steps. First for each scanline it uses a delta function to reduce the amount of different samples. In practice it means it uses the previous pixel to determine the next pixel, and only store the difference (or the pixel above). Which uses the correlation between one pixel and its neighbours to reduce the number of possible samples. Second it uses the deflate algorithm to compress the output. Deflate is used in GZIP for example, and uses Huffman and LZ77 compression.

This works well for drawing and web-illustrations. However ,for images with a ‘high’ amount of noise it might be less ideal, because the noise is somewhat unrelated. Photos often have this noise. PNG has some pre-filters which you can use to actually remove noise, increasing the compression factor of PNG quite a bit. However, this effectly turns PNG into a lossy compression (e.g. the original image can not be completly restored and some data is lost).

I have found that PNG compresses a normal photo to about 1/3 of its original size, give or take some depending on the picture. Compared to JPG, which compresses without much visual disturbance between 1/10 to 1/25, there is a huge gap.

I think it would be possible to create a lossless photo format with better performance (read compression ratio) than PNG by taking into account the special properties of real-life photos. Now all I have to do is figure out is how to I have come up with some points:

• The image processing should take place in subbands (higher frequencies should be split from lower frequencies). This way it might be possible to process the high-frequency noise separate from the lower-frequency image-data.
• It may chop the picture up in smaller subblocks for easier processing. For example 64 x 64 blocks.
• The format focusses primairly on Photos. It would perform poorly on drawings.
• The format should use correlation between pixels, but also between different color-channels.
• The compression happens in two steps. The first prepairs the image-data for compression, reducing the number of possible samples. The second is your plain and average compression like deflate. Somewhat similar to PNG.

I have some ideas about how to approach this. I already build a prototype which failed misserably (compressing only 2/3), however I have some new ideas. I’ll publish the results of the new format here as soon as I have some :).

If you have some ideas for this format, I’d be very interested!

A few posts before this one I posted an item on approximating sines. Unfortunatly the images are gone, but the idea is still the same :-). Now there is an ‘interesting’ detail I wish to share about the sine approximation system. You can use it to extrapolate a single sine without knowing its phase, frequency or amplitude.

So, if you have a sine-wave-form in samples, you can use this algorithm to see how it continues, without any advanced phd level math.The only requirement is that the sine is centered around 0 (e.g. the average of all samples is about 0), and its a single sine. How does this work. The algorithm which was finally found to produce a sine was:

y = 1;
i = 0;
while (thereismoretime) {
i += -y * f;
y += i;
outputSample(y);
}


Where y is the current position of the wave, i is the current inertia, and f the flexibility (or inflexibility actually - because the higher f is the quicker it will try to return to nutral - but i was already used). Now the trick in extrapolation is it get the initial values of i and y, and the constant value of f. To do this you will just need the 3 last samples of the waveform.

Let say that s(t) is the last sample and s(t - 1) is the previous sample.

The y, is simple. It will just be y = s(t). Since the starting position of y is the last sample (or maybe the next-to-last sample).

The initial value of i, is not much more difficult. Since every turn y += i is invoked, it is basically the derative (the difference between this sample, and the previous). This means that i = s(t) - s(t - 1), where t is the last sample index. Lets turn this into a function (handy for later):
i(t) = s(t) - s(t - 1).

Now we only need to deduct f. F is used in the following function:
i += -y * f;.
There are no strangers here, except for f. We can calculate the value of i for almost any place in the array of samples, we can calculate the difference of i too (double derivative): i’(t) = i(t) - i(t - 1). So now we have i’(t) we know that i’(t) = -y * f. Where -y can be replaced with -s(t - 1). Its t - 1, becuase the change i’(t) is over a change which has already happend. Draw the steps out and you will see. So i’(t) = -s(t - 1) * f.

Now we want f, so we divide the whole thing trough s(t - 1), which ends up being:
f = -i'(t) / s(t - 1)
En voila, we are done.

To put this in code (where the s[] array contains the samples):


y = s[t];
i = s[t] - s[t - 1];
f = -(i / s[t - 1]);
//.. loop here


Ofcourse if s[t-1]=0 you have a problem, but I’m sure that if you understood all the stuff above this far, you can find a solution for that either.

Now the real question is: what can you do with it. Simple anwser: No idea. As said before as soon as the wave-form is not a single sine, it won’t work. It would be handy if you have a partial wave-form and you wish to know how it looks. However, I’m sure much smarter people than me have found much better solutions for that using much more solid math :-).

Unfortunatly due to a server upgrade yesterday some images have been lost. I will attemp to recover them. Unfortunatly I do not have a backup of the resoures belonging to the ‘approximating sine’ post.

As posted before, I did som experimenting with 2 dimensional gravity fields. For the sake of speed and just out of curiousity I have ported the gravity field simulation to Chaos Pro, and with succes. Below are some picture of the first gravity field simlations in Chaos Pro.

You can clearly see the repeating patterns. However, I have seen instances of the actual black shape coming back in its subcomponents. Hopefully I’ll be able to find those again. For now enjoy the pictures

A sine is a very natural form. Because its so natural, you would expect it to be simple to generate. And it is! No magic cos or sin function is required if you understand the natural dynamics involved in creating a sine.

So, how is a sine naturally created? Think of a snare from a guitar. The surface (in this case a snare) has a tendency to get back to its nultral (balanced) state. So, the more the vibration devivates from the natural position of the surface the more it has the tendency to get back to the nutral state (or when tourmented enough, it will break). Lets call the position of the ’snare’ x. At position 0 its nutral. So when x = 0, it has no tendency to go anywhere.

However, if the position is 1 (whatever that may be, 1 centimeter, 1 inch, 1 something), it wants to get back with the exact opposite force, which is -1. So, how to code this:

x = 1;
will = -x;

Now we also need to realize that a wave is a ‘progressing’ form. The next state is dependend on the previous state (unlike the sine function which is uncorrelated). So, to generate the waveform we need to simulate time. This can be done by iteration:

x = 1;
while (thereismoretime) {
w = -x; // it wants to go back exactly the amount it deviates
x += w;
}

Now this will create a rather boring graph:

So, what are we missing? Yes! Inertia… the tendency for a body to continue the path it is taking, unless stopped. Lets rename w (will) to inertia (i):

x = 1;
i = 0;
while (thereismoretime) {
i += -x; // substract from the inertia the will to get back to the nutral
x += i;
}


This will cause the following graph:

So, we’re getting there. We now have inertia, and the tendency to always return to the center. However, it is rather crude. Suppose this was sampled at 44kHz then frequency of this wave would be 44kHz / 4, because one wave is 4 ’samples’ long. So what would be 11kHz. What can we do to lower the Hz? Decrease the amount the wave is corrected. So, suppose we change the following line:

i += -x;

into

i += -x * 0.1;

What do we get than?

That starts to look like it. Its not entirely accurate, but it is almost the same. However, lets compare this to a real cosine function.

Looks quite like it huh? I cheated a little bit, since the approximation runs behind, but its quite clear they are (almost) the same. And by decreating the factor it changes per iteration (the 0.1 to for example 0.01) it will probably get more accurate.

The sine and cosine are very natural functions, which occure everywhere in nature. This is why you can use a descrete fourier transform to analyse the frequencies in sampled data. Real life waves have this shape, its natural