Making a Bode plot using the sinusoids data (2025 version)

I used sinusoids to collect some data where on each trial the motor was driven with a single sinusoid of a particular frequency.

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 sin_14.m, created with the command "./create-commands ../sin_14". This program creates a file of matlab commands, in this case sin_14.m. There is nothing magic about this program, it simply manipulates text 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 "sin_14". This makes Matlab look for the file sin_14.m and execute the commands in it.

What do those Matlab commands do?

sin_14.m applies process_pos.m and process_vel.m to the directory of data files (sin_14) 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 ../sin_14/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:

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:

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

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

load fff14.txt
    or
load fff14.mat

The variable fff14 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.