MLS Generator for ATMEGA8

Maximal length sequences (MLS or m-sequences) are useful for system identification and digital communication. For example if a MLS is used to excite the input to a linear time-invariant (LTI) system, the system's impulse response can be determined from the cross correlation of the input and output. The impulse response can be useful in its own right or can be windowed to obtain a quasi-anechoic frequency response.

This code is written for an AVR ATMEGA8, which is an easy-to-use, low-cost microcontroller, but should work with others. It outputs the MLS on pin 18 (PB4). A clock signal is output on pin 19 (PB5), transitioning high on each value. Pin 17 (PB3) is a synchronizing signal, going high for one sample at the start of each loop of the sequence.

The MLS is generated by a 16-bit linear feedback shift register. The "feedback" variable can be changed to other values to generate other sequences. See http://www.ece.cmu.edu/~koopman/lfsr/index.html for other values.

The code is written for gcc with -O3 optimization and the loop runs in a minimum 19 clock cycles if the delay function is removed. The "NOP" (No OPeration) lines must be tuned if you use a different compiler, to ensure that each branch of each "if" statement takes the same number of cycles.

"N_delay" can be used to change the loop's speed. It is currently set for a loop of 73 cycles or approx 219 kHz for a 16 MHz clock.

Download from here. You have to sign up (free), but if you're interested in AVR microcontrollers you should be a member anyway.

Acoustic Tomography

One of the projects that I'm currently working on is a method of measuring temperature and velocity fields using sound, called acoustic tomography. This method relies on the two properties of sound: it's speed is temperature dependent, and the sonic speed is relative to the motion of the fluid (in our case, air). The principle is similar to that used for sonic anemometers, the time taken for a sound to cross a distance is measured in two directions and the mean wind speed and temperature can be determined. If we set up a large number of speakers and microphones around a measurement area, it is possible to reconstruct the temperature and/or velocity fields within the area. This can be used to generate an image or video of the temperature and wind distribution.

This has been done before, however we are looking at much smaller scales in order to measure the flow in street canyons.

For more information, see my paper at the New Zealand Acoustic Society conference. Example code is also available from my section of the Matlab File Exchange.

The image above shows the temperature (colours) and wind velocity (arrows) field for votex shedding behind a cylinder. The top image shows the simulated flow, and the bottom one shows the reconstruction. The video below shows multiple frames of similar data. Unfortunately the resolution is low, if you would like to see the full video, please contact me.

New Slide Background

In preparation for next week's presentation, I got Amy Templeman to design some new Powerpoint slide backgrounds. I really like them, so I thought I'd share. The top one is for Nutaksas; the other is for the University of Auckland.

Upcoming Conferences

Next week (Nov 24-28), Travis will be speaking at two conferences. He will be presenting some recent developments in acoustic tomography at ENZCON 2008, "New Zealand’s leading national conference in electronics", as well as the New Zealand Acoustic Society Conference 2008. Acoustic tomography is a non-invasive method of measuring temperature and flow fields by measuring the time taken for sound to cross an area or volume.

Colleagues from the University of Auckland will also be presenting at both conferences.

Watch this space for details.

PhD Thesis Available

My PhD Thesis, "Online Learning of a Neural Fuel Control System for Gaseous Fueled SI Engines," is now available at lulu.com. This is a beautifully bound 153-page hardcover version of the dissertation including all the papers and a listing of the computer code used. If the price is too expensive for you, please feel free to email me for a free electronic copy. I'm not posting direct links for the first couple years so I can get some feedback as to who is interested in the work.

The cover was designed by Amy Templeman.

Internoise 2008

Travis will be presenting two papers at the Internoise 2008 conference next week. The papers are entitled "Experimental Characterization of Sound Propagation in a Dense New Zealand Forest" and "Long Range Identification of Wildlife Using Phased Arrays of Microphones: A Feasibility Study".

The first is a study which measures the sound attenuation of the very dense forests on New Zealand; that is, how quickly a sound decreases in sound pressure level over distance. We found that sound is attenuated more in the native forests of the Hunua Ranges than in any forest in the literature, at approximately 0.5 dB/m.

The second paper investigates the feasibility of locating birds by their calls, with the particular application of finding the location of kokako in the Hunua Ranges. Due to the high attenuation of sound through the forest, we investigated a scheme which would place arrays of microphones on towers above the canopy. This was compared to a more traditional method of installing a grid of individual microphones near the forest floor. This work used the forest reverb model seen earlier on this blog.

Details to follow.

Thesis Defense Update

I am happy to announce that, as of Wednesday, I am now officially Dr. Travis Wiens. My thesis defense went well and had a lot of good questions. My dissertation, entitled "ONLINE LEARNING OF A NEURAL FUEL CONTROL SYSTEM FOR GASEOUS FUELED SI ENGINES" will be available shortly. If you'd like a copy, please feel free to contact me.

I'd also like to announce that Nutaksas will now be adding minor surgery to our list of services available.

University of Reading Seminar

Travis recently presented a seminar at the Department of Meteorology at the University of Reading. The topic was the introduction of the acoustic tomography system currently in development at the University of Auckland. The slide show is available here.

Bath/ASME Symposium on Fluid Power and Motion Control

Travis recently presented the Recurrent Generalized Neural Network with training via the Complex Method at the Bath Conference. Details are available in a previous post. The slides may be found here.

One of the many uses of artificial neural networks is for nonlinear black-box dynamic modeling. This is the process of predicting the response of a dynamic system to its inputs. A number of researchers have attempted this using gradient methods of selecting the neural network parameters, to varying levels of success. However, most will tell you that the most difficult (and accuracy effecting) process is approximating the derivatives. For example, Backpropagation Through Time uses an algebraic "unrolling" of the IIR filter's derivative into its approximate FIR form, although the algebra quickly gets out of hand for systems with long impulse responses.

An easier (and arguably better) method is to use a non-derivative optimization method. My favorite is the Complex Method (matlab code), which is a simpler method of the Nelder-Mead Simplex Method. This method requires the ability to calculate a "fitness" for each set of parameters (network weights), which will be maximized. In the case of dynamic system modelling, we want to minimize the error between the neural network estimate of the system response to a given input and the real system's measured response, so we use the negative error. The fitness of a large number of paremeter sets are calculated. The worst is removed from the set, and a new set of parameters is generated, based on the fitness of the previous points. This process is repeated for a number of generations until a desired accuracy is reached.

We (myself and University of Saskatchewan colleagues Rich Burton, Doug Bitner and Greg Schoenau) have put together a paper showing how this works, to be presented at the Bath Symposium on Power Transmission and Motion Control. This paper shows the performance of the complex method when applied to the modelling of real data from a load-sensing hydraulic pump. We found the complex method to be more accurate and much faster than a previous method used on the same data.

You can get a preprint of the paper here, and the Matlab code here. As always, comments are appreciated.

Lowest Common Mulitple of Periods within an Octave

In digital signal processing, it is sometimes required that one calculate the lowest common multiple of a number of periods. For example, imagine we have two sine waves with periods of 6 and 8 samples, and wish to calculate the discrete Fourier transform to determine the phase and amplitude of these two frequencies. Since the DFT assumes that the wave is periodic, we should use a DFT length of 24, the lowest common multiple of 6 and 8.

I came upon the inverse problem when designing a system to generate square waves, which would be analyzed with a DFT. However this problem had a further constraint, I needed to avoid interfering harmonics. A square wave has odd harmonics of the fundamental frequency. That is, for a 1 Hz square wave, there will be components at 3, 5, 7, ... Hz. In order to avoid interference between frequencies, I needed to find the shortest DFT length k that was a multiple of N periods, where the ratio of the largest and smallest periods was less than 3.

It took a bit of time to brute force the problem, so I'll include the results here:
N=3 k=12 p=[ 2 3 4 ]
N=4 k=60 p=[ 2 3 4 5 ]
N=5 k=120 p=[ 4 5 6 8 10 ]
N=6 k=240 p=[ 8 10 12 15 16 20 ]
N=7 k=360 p=[ 8 9 10 12 15 18 20 ]
N=8 k=720 p=[ 8 9 10 12 15 16 18 20 ]
N=9 k=840 p=[ 14 15 20 21 24 28 30 35 40 ]
N=10 k=1680 p=[ 14 15 16 20 21 24 28 30 35 40 ]
N=11 k=2520 p=[ 24 28 30 35 36 40 42 45 56 60 63 ]
N=12 k=2520 p=[ 24 28 30 35 36 40 42 45 56 60 63 70 ]
N=13 k=5040 p=[ 28 30 35 36 40 42 45 48 56 60 63 70 72 ]
N=14 k=5040 p=[ 28 30 35 36 40 42 45 48 56 60 63 70 72 80 ]
N=15 k=10080 p=[ 56 60 63 70 72 80 84 90 96 105 112 120 126 140 144 ]
N=16 k=10080 p=[ 56 60 63 70 72 80 84 90 96 105 112 120 126 140 144 160 ]
N=17 k=15120 p=[ 48 54 56 60 63 70 72 80 84 90 105 108 112 120 126 135 140 ]
N=18 k=25200 p=[ 60 63 70 72 75 80 84 90 100 105 112 120 126 140 144 150 168 175 ]
N=19 k=27720 p=[ 55 56 60 63 66 70 72 77 84 88 90 99 105 110 120 126 132 140 154 ]
N=20 k=30240 p=[ 96 105 108 112 120 126 135 140 144 160 168 180 189 210 216 224 240 252 270 280 ]
N=21 k=50400 p=[ 120 126 140 144 150 160 168 175 180 200 210 224 225 240 252 280 288 300 315 336 350 ]
N=22 k=55440 p=[ 60 63 66 70 72 77 80 84 88 90 99 105 110 112 120 126 132 140 144 154 165 168 ]
N=23 k=55440 p=[ 60 63 66 70 72 77 80 84 88 90 99 105 110 112 120 126 132 140 144 154 165 168 176 ]
N=24 k=83160 p=[ 105 108 110 120 126 132 135 140 154 165 168 180 189 198 210 216 220 231 252 264 270 280 297 308 ]
N=25 k=110880 p=[ 120 126 132 140 144 154 160 165 168 176 180 198 210 220 224 231 240 252 264 280 288 308 315 330 336 ]
N=26 k=110880 p=[ 120 126 132 140 144 154 160 165 168 176 180 198 210 220 224 231 240 252 264 280 288 308 315 330 336 352 ]
N=27 k=166320 p=[ 210 216 220 231 240 252 264 270 280 297 308 315 330 336 360 378 385 396 420 432 440 462 495 504 528 540 560 ]
N=28 k=166320 p=[ 210 216 220 231 240 252 264 270 280 297 308 315 330 336 360 378 385 396 420 432 440 462 495 504 528 540 560 594 ]
N=29 k=166320 p=[ 210 216 220 231 240 252 264 270 280 297 308 315 330 336 360 378 385 396 420 432 440 462 495 504 528 540 560 594 616 ]

Now, while the square wave has no even harmonics, these harmonics can be generated by non-linear processes, such as the saturation of vacuum tubes. So the same analysis was performed, constraining the ratio of the maximum and minimum periods to be less than 2 (or within one octave). Here are the results:
N=3 k=60 p=[ 3 4 5 ]
N=4 k=180 p=[ 9 10 12 15 ]
N=5 k=720 p=[ 8 9 10 12 15 ]
N=6 k=840 p=[ 20 21 24 28 30 35 ]
N=7 k=2520 p=[ 35 36 40 42 45 56 60 ]
N=8 k=2520 p=[ 35 36 40 42 45 56 60 63 ]
N=9 k=5040 p=[ 35 36 40 42 45 48 56 60 63 ]
N=10 k=10080 p=[ 56 60 63 70 72 80 84 90 96 105 ]
N=11 k=15120 p=[ 70 72 80 84 90 105 108 112 120 126 135 ]
N=12 k=27720 p=[ 55 56 60 63 66 70 72 77 84 88 90 99 ]
N=13 k=27720 p=[ 55 56 60 63 66 70 72 77 84 88 90 99 105 ]
N=14 k=55440 p=[ 55 56 60 63 66 70 72 77 80 84 88 90 99 105 ]
N=15 k=83160 p=[ 165 168 180 189 198 210 216 220 231 252 264 270 280 297 308 ]
N=16 k=83160 p=[ 165 168 180 189 198 210 216 220 231 252 264 270 280 297 308 315 ]
N=17 k=110880 p=[ 198 210 220 224 231 240 252 264 280 288 308 315 330 336 352 360 385 ]
N=18 k=166320 p=[ 165 168 176 180 189 198 210 216 220 231 240 252 264 270 280 297 308 315 ]
N=19 k=221760 p=[ 252 264 280 288 308 315 320 330 336 352 360 385 396 420 440 448 462 480 495 ]
N=20 k=277200 p=[ 264 275 280 300 308 315 330 336 350 360 385 396 400 420 440 450 462 495 504 525 ]
N=21 k=332640 p=[ 198 210 216 220 224 231 240 252 264 270 280 288 297 308 315 330 336 352 360 378 385 ]
N=22 k=360360 p=[ 252 260 264 273 280 286 308 312 315 330 360 364 385 390 396 420 429 440 455 462 468 495 ]
N=23 k=554400 p=[ 264 275 280 288 300 308 315 330 336 350 352 360 385 396 400 420 440 450 462 480 495 504 525 ]
N=24 k=720720 p=[ 504 520 528 546 560 572 585 616 624 630 660 693 715 720 728 770 780 792 819 840 858 880 910 924 ]
N=25 k=720720 p=[ 504 520 528 546 560 572 585 616 624 630 660 693 715 720 728 770 780 792 819 840 858 880 910 924 936 ]
N=26 k=720720 p=[ 504 520 528 546 560 572 585 616 624 630 660 693 715 720 728 770 780 792 819 840 858 880 910 924 936 990 ]
N=27 k=720720 p=[ 504 520 528 546 560 572 585 616 624 630 660 693 715 720 728 770 780 792 819 840 858 880 910 924 936 990 1001 ]

So, for example, the period of the shortest repeating waveform with 12 frequencies within an octave is 15120. If we allow frequencies up to the 2nd harmonic of the lowest frequency, the shortest period is 2520.

Radiohead "House of Cards" Surface Fit

As you may know, Radiohead graciously released some of the 3D data from their "House of Cards" video on Google Code. Each frame of the data has approximately 6000 points of x,y,z position data, as well as a colour intensity value (like it's in black and white), spaced non-uniformly through the space.
I've been using some Radial Basis Function Networks to approximate some functions and wanted to see how they worked to interpolate the non-uniform data only a uniform grid. Using a Thin Plate Spline basis function, I created the surface in the image. A thin plate spline is a surface that fits exactly on the control points and the surface between the control points bends as if it was a thin plate of metal. There is also, luckily, a closed form solution for the parameters.
If you're interested, the Matlab code is available of the Matlab Central File Exchange.
Next: video.

Forest Reverb Model

For your auditory pleasure, Travis Wiens and Nutaksas Research are pleased to announce the release of a Forest Reverberation Modeller. This package of Matlab functions allows you to simulate the impulse response of an arbitrary forest of any number of trees of any size at any location, as well as source and listener lcoations.
Based on Morse's theory, each tree is treated as an acoustically hard cylinder and the effect of each path from source to receiver via any number of trees is calculated. The maximum number of scatterings may be limited to cut down on excessive calculation times.
The idea was based on Kyle Spratt's Treeverb proposal, although it was developed independently and gives different results.
Download the code from Matlab Central File Exchange, and feel free to check out the other files available. Watch this space for a more technical description.