Making a Bode plot using the sinusoids data.

The latest version of the sinusoids program is sinusoids9

It does a slightly better job picking the frequencies to sample, so they are always exactly in one of the FFT bins. This page was written using sinusoids8 data.

I used sinusoids8 to collect some data where on each trial the motor was driven with a single sinusoid of a particular frequency. The resulting data is in lab1/sin8_1.log. I used parse-data (lab1/parse-data.c) to break that into individual files. On linux I compile parse-data.c by typing "make" in the lab1 directory. I run the program with "./parse-data sin8_1", assuming the putty.log file has been renamed sin8_1.log. The program creates a directory with that name (in this case sin8_1) and then takes each section of the putty log corresponding to a trial running a sinusoid with a constant frequency and creates a Matlab-readable data file with a name like f4_859.dat (for frequency 4.859Hz). The columns are: timestamp (in microseconds), left wheel angle, left wheel angular velocity, left motor command, right wheel angle, and right wheel angular velocity (in these experiments I sent the same command to both motors). Wheel angles are in encoder units and wheel velocities are in encoder units per second. parse-data also creates a text file (such as f4_859.txt) with text information relevant to that trial. The program also makes a file named files.txt, in the data directory which stores a list of all data files. For each trial there is the name of the data file, the name of the text file, the number of points sampled, the frequency, and the maximum amplitude of the sine function:

../sin8_1/f4_859.dat ../sin8_1/f4_859.txt 2059  4.8591 255.00

Processing the data in Matlab

I like to do my data processing and modeling in its own directory. In the lab1 directory there is a directory called "bode-plots", which I will use for this example. As you try different approaches, or process different data, you may want to create separate data processing and/or modeling directories with whatever name you choose.

It is useful to create a Matlab command file that loads all the data files created above into Matlab, so you don't have to do it manually. In the bode-plots directory (to get there type "cd bode-plots" on linux) there is an example command generator: create-commands.c, and the corresponding command file sin8_1.m, created with the command "./create-commands ../sin8_1". This program creates a file of matlab commands, in this case sin8_1.m. There is nothing magic about this program, it simply manipulates text from the files.txt file to create a Matlab command script. You can modify create-commands.c to match your own purposes, or edit the script it produces. Use "make create-commands" to compile create-commands.c.

Now you can start matlab. Make sure Matlab's current directory is the bode-plots directory, and then type "sin8_1". This makes Matlab look for the file sin8_1.m and execute the commands in it.

What do those Matlab commands do?

sin8_1.m applies process_pos.m and process_vel.m to the directory of data files (sin8_1) where each data file consists of driving the system at a particular frequency. The results are stored in an array fff (I am not good on making up useful names).

load ../sin8_1/f0_099.dat
[pml ppl pmr ppr] = process_pos( f0_099, 0.0992, 0.001 );
fff(1,1) = 0.0992;
fff(1,2) = pml;
fff(1,3) = ppl;
fff(1,4) = pmr;
fff(1,5) = ppr;
[vml vpl vmr vpr] = process_vel( f0_099, 0.0992, 0.001 );
fff(1,6) = vml;
fff(1,7) = vpl;
fff(1,8) = vmr;
fff(1,9) = vpr;
...

process_pos()

The function process_pos.m calculates the magnitude ratio and phase lag between input command and position of each wheel. The overall method is to take the FFT of both input and outputs, and calculate the magnitude ratio and phase lag at the frequency the system was driven at. process_pos.m uses fft2cga.m to call Matlab's FFT routine. See the comments in fft2cga.m for an explanation of what it does.

The program process_pos.m converts the encoder readings into radians. A full revolution (2*pi radians) is 84480 counts using the PicoEncoder library. There are 11 counts/revolution in the actual encoder. Quadrature encoding multiplies this by 4. The gear ratio is 30, and the PicoEncoder multiples by 64. (11*4*30*64 = 84480). So the scaling factor to turn encoder readings into radians is 2*pi/84480 (and the same scaling factor converts encoder counts/second into radians/second).

process_vel()

There is an analogous progam process_vel which calculates the magnitude ratio and phase lag of how the command affects the velocity in encoder units.

Want to see the plots and printouts?

The programs process_xxx.m print and plot a lot of stuff so fast you can't see what it is. Uncomment the pause statement at the end of process_vel.m to see the printouts and plots for each frequency. After each pause, you may need to click your mouse on the matlab window and then type a character (a space or newline) to continue.

Checking for transmission (putty) errors

process_pos first plots the difference of the adjacent samples times (microseconds) to check for transmission errors or missing samples. Ideally, the sampling interval should always be 1000 (microseconds, which is 1 millisecond). You can see there is usually a little variation between 996 and 1004:

For lower frequencies you get solid blue:

Here is an example where a sample was lost in transmission:

Sometimes a line gets garbled during transmission with extra characters, typically at the begining of a line. This causes a huge error in the sample interval plot.

16690003 -1725850 205456 165 8832320 219369
616691000 -1725645 205068 166 8832541 222414
16692003 -1725437 206684 166 8832761 216734

83533004 1379750 -377778 -255 4442919 -409225
5330083534004 1379375 -376096 -255 4442516 -406160
83535001 1378994 -381851 -255 4442112 -409091

You can edit the log file xxx.log to remove the extra characters. Then re-run parse-data to regenerate the individual files for each trial.

Sometimes something is missing, and matlab flags an error:

43025003 901582 -249512 -139 2990091 -251721
43026002 901336 -247104 -138 2989833 -
43027001 901089 -247104 -138 2989581 -253214

You can edit the log file xxx.log to fix the glitch.

43025003 901582 -249512 -139 2990091 -251721
43026002 901336 -247104 -138 2989833 -252500
43027001 901089 -247104 -138 2989581 -253214

Then re-run parse-data to regenerate the individual files for each trial.

Sometimes there is one missing sample, or a small gap, and the sampling intervals are reasonable after the gap. You can interpolate to fill the gap. But sometimes the gap is large and the timing is screwed up after the gap. Note the gap below between 64417000 and 64543960 microseconds of about 127 milliseconds, and that the timing is off afterwards (the last 3 digits of the microsecond clock are not 00x as they usually are, and the intervals are not 1000+/-4). I would throw away a trial that had this big a glitch by deleting the corresponding line in the files.txt file.

64416000 593289 -389450 -254 5833838 -409249
64417000 592912 -379964 -254 5833434 -405704
64543960 543274 -390949 -248 5781864 -406172
64544810 542938 -394444 -248 5781522 -403509

Processing command to angles data

process_pos then plots the command (blue) and angles (left is reddish-orange and right is yellowish-orange). At 0.1Hz, the command and position are about the same size on the graph (I scaled the command to make this true). Note that the wheel positions each often have a different drift in addition to the sinusoidal movement.

At a medium frequency (1Hz), the angle excursions are reduced and more lagged relative to the command.

At a high frequency (10Hz), the angle excursions are too small to see, and more lagged.

Checking for bugs and glitches

The other two position plots are for debugging. They show the FFT power plot, and the command and angles with the angles scaled by the magnitude so they should have the same amplitude as the command (if not there is a bug).

FFT power spectrum plots

Here are the FFT power spectrum plots. To make sure the plots for the angles are clearly visible at high frequencies, I plot the log of the power spectrums. The spectrums are centered on the frequency of the sinusoidal command.
0.1Hz:

1Hz:

10Hz:

Checking the magnitude ratios

process_pos also plots the command and angles scaled by the estimated gain ratio, so the angles and command should have the same amplitude. If not, there is a bug.
0.1Hz:

At a medium frequency (1Hz), the angle excursions are more lagged relative to the command. The drift of each wheel is more visible.

At a high frequency (10Hz), the angle excursions are even more lagged to about 90 degrees (the sinusoids are almost perfectly out of phase).

Saving the data in the array fff

The frequency, magnitude ratio, and phase lag for commands to positions (and later commands to velocities) are stashed in the fff array.

Making the command to positions Bode plots

The Bode plot figures for the motor command to wheel angle are generated by the following commands (which are in plot_bode.m):

figure(1)
hold off
loglog( fff(:,1), fff(:,2), '.' )
hold on
loglog( fff(:,1), fff(:,4), '.' )
legend( 'left', 'right', 'Location', 'SouthWest' ) 
title( "Magnitude ratio for command to position (sin8_1)" )
xlabel( "Log frequency (Hz)" )
ylabel( "Log ratio" )           

figure(2)
hold off
semilogx( fff(:,1), fff(:,3), '.' )
hold on
semilogx( fff(:,1), fff(:,5), '.' )
legend( 'left', 'right', 'Location', 'SouthWest' )
title( "Phase for command to position (sin8_1)" )
xlabel( "Log frequency (Hz)" )
ylabel( "Phase" )

There is good agreement between the left and right sides.

Processing command to angular velocity data

process_vel.m does the same thing for angular velocities that process_pos.m does for angles. The angular velocities in encoder counts/second are converted to radians/second. A full revolution (2*pi radians) is 84480 counts when using the PicoEncoder library, so the scale factor is 2*pi/84480.

process_vel plots the command (blue) and angular velocities (left is reddish-orange and right is yellowish-orange). At 0.1Hz, the command and velocity are almost the same (I scaled the command to make this true).

At a medium frequency (1Hz), the angular velocities are reduced and more lagged relative to the command, but much less than the positions.

At a high frequency (10Hz), the angular velocities are even smaller and more lagged.

Checking for bugs and glitches

The other two velocity plots are for debugging. They show the FFT power plot, and the command and angular velocities with the angular velocities scaled by the magnitude so they should have the same amplitude as the command (if not there is a bug).

FFT power spectrum plots

Here are the FFT power plots. To make sure the plots for the angular velocities are clearly visible at high frequencies, I plot the log of the power. The spectrums are centered on the frequency of the command sinusoid. 0.1Hz:

1Hz:

10Hz:

Checking the magnitude ratios

process_vel also plots the command and angular velocities scaled by the estimated gain ratio, so the angular velocities and command should have the same amplitude. 0.1Hz:

At a medium frequency (1Hz), the angular velocities are a little bit more lagged relative to the command.

At a high frequency (10Hz), the angle excursions are more lagged.

Saving the data in the array fff

The command to velocity magnitude ratios and phase lags are also stashed in the array fff.

Making the command to velocities Bode plots

The Bode plot figures for the motor command to wheel angular velocity are generated by the following commands (which are in plot_bode.m):

figure(3)
hold off
loglog( fff(:,1), fff(:,6), '.' )
hold on
loglog( fff(:,1), fff(:,8), '.' )
legend( 'left', 'right', 'Location', 'SouthWest' ) 
title( "Magnitude ratio for command to velocity (sin8_1)" )
xlabel( "Log frequency (Hz)" )
ylabel( "Log ratio" )           

figure(4)
hold off
semilogx( fff(:,1), fff(:,7), '.' )
hold on
semilogx( fff(:,1), fff(:,9), '.' )
legend( 'left', 'right', 'Location', 'SouthWest' )
title( "Phase for command to velocity (sin8_1)" )
xlabel( "Log frequency (Hz)" )
ylabel( "Phase" )

There is good agreement between the left and right sides.

Saving the fff array

To save the frequency response data, so I can compare it to that of possible models, I type:

fff8_1 = fff;
save('fff8_1.txt','fff8_1','-ascii')
save('fff8_1.mat','fff8_1','-mat')

Here is the ASCII (text) file fff8_1.txt and the more compact binary file which represents numbers at their full resolution: fff8_1.mat. These files can be loaded into a new Matlab session using

load fff8_1.txt
    or
load fff8_1.mat

The variable fff8_1 will have the frequency response data for sin8_1.

Is the PicoEncoder velocity an accurate estimate of the true velocity?

If the velocity sinusoid is cos(freq*t), the corresponding position (the integral of velocity) is sin(freq*t)/freq. The gain of velocity to position is 1/freq, and the phase lag is -90degrees (-pi/2=-1.57). Let's check if this is true for the PicoEncoder velocity estimate.

Plotting the magnitudes of the command to position and command to velocity on the same plot:

The difference of these two plots (on the log scale) should be a straight line on the log scale, which is almost true except at high frequencies when the magnitude and phase estimates aren't very good.

Plotting the phases of the command to position and command to velocity on the same plot:

The difference of these two plots should be a constant (-pi/2) which is true up to about 1Hz, and then the command to velocity has too much lag, so its lead over the position decreases as the frequency increases.

The PicoEncoder velocity estimate does a good job estimating the magnitude of the velocity, but it is delayed relative to the true velocity. What should we do about that?

MATERIAL BELOW NOT UPDATED, LAST YEAR'S INFO FOLLOWS

What changes as the battery runs down?

Run32 ran the battery down from 6.9 to 4.13. DO NOT DO THIS YOURSELF. RECHARGE THE BATTERY ONCE IT GOES BELOW 7 VOLTS!!!! Running the battery this low damages the battery. Fully recharge your battery early and often. The battery voltage is the 3rd column below (first is file name, second is frequency). Once the battery voltage goes below 6.2volts, it collapses to roughly 4V. From the run32 files.txt:

../run32/f02x0 2.0 6.32 767 3000
../run32/f02x1 2.1 6.2 753 2951
../run32/f02x2 2.2 6.28 762 2901
../run32/f02x3 2.3 6.28 763 2851
../run32/f02x4 2.4 6.19 752 2801
../run32/f02x5 2.5 4.01 487 2751
../run32/f02x6 2.6 6.2 753 2701
../run32/f02x7 2.7 4.03 489 2651
../run32/f02x8 2.8 6.3 765 2601
../run32/f02x9 2.9 6.2 753 2550
../run32/f03x0 3.0 4.04 490 2501
../run32/f03x1 3.1 4.04 491 2452
../run32/f03x2 3.2 4.03 489 2400
../run32/f03x3 3.3 4.05 492 2351
../run32/f03x4 3.4 4.05 492 2301
../run32/f03x5 3.5 4.06 493 2251
../run32/f03x6 3.6 4.04 490 2201
../run32/f03x7 3.7 4.05 492 2151
../run32/f03x8 3.8 6.24 758 2101
../run32/f03x9 3.9 4.08 495 2051
../run32/f04x0 4.0 4.09 497 2001

Here are the Bode plots for left side for run32 (fff32.txt or fff32.mat) and run33 plotted together. I stopped run32 early because of the low voltage, and fully recharged the battery before run33. The low battery voltage lowered the magnitude plot by 15%. It should have been lowered by the ratio of the starting voltages: 6.9/8.3 (17%), which is close. So it appears the voltage affects the magnitude plot linearly. Above 2.5Hz we see the collapse to 4V reduces the magnitude more. The voltage seems to have little effect on the phase plot until the collapse to 4V, which we don't care about since we won't ever go this low again!

These plots were generated by:

figure(1)
loglog( fff33(:,1), fff33(:,2), '.' )
hold
loglog( fff32(:,1), fff32(:,2), 'r.' )
title( "Magnitude ratio angle/torque" )         
xlabel( "Log frequency (Hz)" )          
ylabel( "Log ratio" )                           
legend('run33','run32')
figure(2)
semilogx( fff33(:,1), -fff33(:,3), '.' )        
hold
semilogx( fff32(:,1), -fff32(:,3), '.' )
title( "Phase lag of angle wrt torque" )
xlabel( "Log frequency (Hz)" )          
ylabel( "Phase lag" )                     
legend('run33','run32')

However, here is some contradictory evidence. Run 33 did the trial from low to high frequencies, starting with a fully charged battery. Run34 did the trials in the reverse order, from high to low frequencies, also starting with a fully charged battery.

Here are the magnitude plots for both the left and right sides.

Here are the phase plots for both the left and right sides.

The phase plots do show a dependence on voltage here, as do the magnitude plots. Note that the center frequencies were done at about the same voltages. A better experiment would be to do several sweeps low to high frequencies as the voltage went down, to get a voltage difference at the center frequencies. These plots were generated using:

figure(1)
semilogx(fff33(:,1), fff33(:,6), '.' )
hold
semilogx(fff34(:,1), fff34(:,6), '.' )
title( "Voltage when the frequency sampled" )
xlabel( "Log frequency (Hz)" )          
ylabel( "voltage" )
legend('run33','run34','Location','west')
figure(2)
loglog( fff33(:,1), fff33(:,2), '.' )
hold
loglog( fff34(:,1), fff34(:,2), '.' )
loglog( fff33(:,1), fff33(:,4), 'k.' )
loglog( fff34(:,1), fff34(:,4), 'm.' )
legend('run33 left','run34 left','run33 right','run34 right')
title( "Magnitude ratio angle/torque" )
xlabel( "Log frequency (Hz)" )
ylabel( "Log ratio" )
figure(3)
semilogx( fff33(:,1), -fff33(:,3), '.' )
hold
semilogx( fff34(:,1), -fff34(:,3), '.' )
semilogx( fff33(:,1), -fff33(:,5), 'k.' )
semilogx( fff34(:,1), -fff34(:,5), 'm.' )
legend('run33 left','run34 left','run33 right','run34 right')
title( "Phase lag of angle wrt torque" )
xlabel( "Log frequency (Hz)" )
ylabel( "Phase lag" )

Unless we can model the dependence on voltage, the spread of these frequency responses gives us an idea of how much uncertainty about the model we should include.