page 1 of 8 pages

MTRX1701 Introduction to Mechatronic Engineering

Assignment 3: Working with sensor data in Python (Rev A) 2021

Notes on the Assignment

1. This assignment is worth 12% of your final mark in MTRX1701.

2. You will have the opportunity to ask the tutors for assistance during the tutorial

sessions in weeks 8–10. You should spend approximately one hour preparing for each

of these tutorials so that you can benefit fully from them.

3. The tutors will not assist you further unless you can provide real evidence you have

attempted the questions prior to the tutorials. Beyond the tutorial sessions, it is

estimated that you will need sixteen hours to complete the assignment (about 21

hours total).

4. You may discuss the assignment with your peers, however all written work submitted

for assessment must be strictly your own unless cited.

5. All work must be completed electronically (i.e. typed) including all diagrams.

6. You may complete this assignment using Python (this document) or MATLAB (see

the document “Assignment 3: Working with sensor data in MATLAB.pdf”.

Assignment Submission

1. Submit your assignment electronically as a PDF file via the “Assignments” page on

the MTRX1701 Canvas site.

2. Your assignment submission should be in the form of a report which includes answers

to the various questions and uses plots and code fragments when necessary, to back-up

your answers.

3. All your code must be submitted as appendices to your report, with one appendix per

question. Start each appendix on a new page, and designate appendices as A, B, etc.

4. Note that:

a. All plots in your report must have axis labels including units, a plot title and a

legend if there are multiple lines on one plot.

b. All plots and code fragments in the body of your report must be presented as

numbered figures with figure captions.

c. All code must be set in a fixed width font such as Courier.

5. The assignment is due at 11:59 pm, Saturday 15 May (end of week 10).

6. Late submission will be penalised by deducting 5% of the maximum mark for each

calendar day after the due date. Zero will be awarded after ten calendar days late.

7. Submit early enough to ensure that your submission is processed by 11:59 pm.

8. Keep the email that will be generated automatically by the system as proof of your

submission.

page 2 of 8 pages

Preparation

To complete this assignment, you will need to download and install some version of Python 3

along with the packages numpy, scipy and matplotlib. You can either use the official

distribution of Python or an alternative distribution such as Anaconda.

To install the official distribution of Python (download < 1GB):

1. Visit: https://www.python.org/downloads/

2. Select the installer according to your operating system and architecture.

3. Run the installer and follow the instructions.

4. Use pip (bundled with the official distribution of Python) to install and numpy,

scipy and matplotlib by opening a command window (terminal) and running:

pip install numpy scipy matplotlib

To install the Anaconda distribution of Python (2.5–3 GB):

1. Visit https://docs.anaconda.com/anaconda/install/

2. Download the installer according to your operating system and architecture.

3. Run the installer and follow the instructions accordingly. The installation should be

straightforward.

4. Anaconda comes pre-bundled with numpy, scipy and matplotlib. Once the

installation is complete, you can begin programming by launching the Spyder IDE or

your favourite code editor.

The current image on PCs in the School of IT and School of AMME is Anaconda3 5.0.1

(Python 3.6.3) x86.

The University also provides the option to use Anaconda Navigator on the Citrix Workspace.

To use this option,

1. Go to https://uniconnect.cloud.com/

2. Login using your unikey and password. The username must be entered in the

following format: shared\unikey; for example. shared\abcd1234

3. Select Anaconda Navigator from the App list.

4. Launch the Spyder IDE to start programming.

You can also use your Ed workspace on the ENGG1810 Ed web site to do this assignment if

you really want to, but it is recommended that you to the assignment on your own computer.

Context of the Assignment

Mechatronic engineering involves the management of data. A sensor is a source of data that

must be read, converted into appropriate engineering units, and then used in some way. If the

sensor is part of a control loop the data processing must be completed well within the loop

time. Understanding and working with data is critical to ensuring the reliability and

correctness of most mechatronic systems.

This assignment explores some of the important aspects of managing and interpreting data by

examining data from a variety of sensors associated with a mobile robot, Groundhog.

Groundhog was designed and built at the ACFR (Australian Centre for Field Robotics). Each

student is expected to work individually to understand the concepts and produce answers to

the questions. Students may be asked by their tutor to demonstrate their solutions during the

tutorials. You will use Python as a tool to load, modify and visualise data.

page 3 of 8 pages

When a signal is sampled by a computer-controlled data acquisition or control system the

sampling frequency is only approximately constant; deviations from a constant sampling rate

can be expected to decrease with faster, higher-quality computer hardware, just as the cost of

such hardware will inevitably be higher. A further complication can arise if data are being

transmitted on a network since in that case packets of data may occasionally be lost. When

working with data it is always wise to understand what the data represent, and to inspect the

data for anomalies such as irregular sampling or data losses. Your code should always allow

for such imperfections in data.

Question 1: Signal sampling [22 marks]

A sensor ‘sample’ is a single reading from a sensor that measures some signal – one datum.

When a sequence of samples is acquired over time the sampling rate is very important. This

question investigates the effect of sampling frequency on the accuracy with which the original

signal can be reconstructed from the sampled data.

For a constant sampling rate, the sampling interval ts is defined as the time in seconds

between successive samples. The sampling frequency fs is the number of samples acquired per

second (1/sec, or Hertz), or the reciprocal of the sampling interval:

fs = 1 / ts.

For example, if a sample is taken every 10 ms, the sampling frequency is = 1/0.010 = 100 Hz.

How should the sampling frequency be chosen? One may choose to sample a signal at a faster

rate (frequency) so that:

The sampled data contain1 more detail, allowing ‘zooming in’ on regions of interest, or

Any noise that corrupts the signal can be averaged out more precisely.

Conversely, it may be necessary to sample less frequently because:

A sensor may have a limited sampling frequency, or

The amount of data at a higher rate may overwhelm processing capacity or overflow the

available storage capacity.

The minimum rate at which a signal needs to be sampled depends on the highest frequency

of interest in the signal that we are measuring. The following questions explore why this is

important. Follow the steps outlined and address the important points in your report.

Creating a dataset using Python

You will use Python’s NumPy library to help generate an idealised, noise-free set of data, and

Python’s Matplotlib to help create plots.

0. The libraries first need to be imported as follows.

import matplotlib.pyplot as plt

import numpy as np

1. Generate a vector of data t representing a set of time values in seconds. Create five

seconds of time data at an interval of 0.001 seconds (1 ms) as follows.

# np.arange(t_start, t_end+t_step, t_step)

np.arange(0, 5+0.001, 0.001)

1 Yes, “data” is plural; the singular form is “datum”.

page 4 of 8 pages

2. Generate a vector y containing a sinusoidal signal with a frequency of 8 Hz.

y = np.sin(np.multiply(8*(2*np.pi), t))

3. Because the vector t only contains data at intervals of 1 ms, the sinusoidal signal y is

effectively sampled at that interval. What sampling frequency corresponds to this

sampling interval? How many samples are taken of each complete period of the sine

curve?

4. Plot the sampled signal amplitude y against time t using

plt.plot(t, y)

Add labels to the plot axes and add a title to the plot. Search the Matplotlib

documentation at https://matplotlib.org/ > [Documentation] to find out how to do this.

5. Zoom in to the graph by creating a new graph with a limited ranges on the axes as

follows, to display the first second of data. How many times does the signal repeat? Is

this what you expected?

# plt.xlim(left_lim, right_lim)

# plt.ylim(bottom_lim, top_lim)

Sampling from this Dataset using Python

What will happen when the 8 Hz signal is sampled at different sampling frequencies?

Investigate this by creating more time vectors, say t1, t2, etc., and generating responses from

the functions as in step 2 to create the corresponding outputs y1, y2, etc.

6. What happens if the 8 Hz signal is sampled at 20 Hz, 9 Hz, 8 Hz, 7 Hz and 5 Hz?

To see additional plots overlaid on the original plot axes, keep plotting after the

command plt.figure(fig_num, figsize=FIG_SIZE). This command retains the

previous plot and draws any subsequent plots on the same axes. The command

plt.subplot creates multiple plots in one figure environment.

7. Explain what is happening. Hint: try plotting with markers: plt.plot(t, y,

marker='o') to see where the samples are falling. Use plt.legend() to show

which markers are which.

8. Research the Nyquist rate and find what the minimum sampling frequency should be for

an 8 Hz signal. Sample the function at the Nyquist rate plus a small amount and show

your results.

Question 2: Noisy data and filtering using Python [26 Marks]

In this section, we will look at some data collected from a millimetre-wave radar. The data is

available in the comma-separated variable (CSV) file radar-2021.csv. This file contains

512 columns of numbers, with each column corresponding to the power received by the radar

(in arbitrary “radar units”) at a particular range. The radar rotates in azimuth at a constant

speed and the samples are taken at a constant sampling rate that corresponds to an increment

of 0.144 degrees of azimuth.

Important note: The outline solution provided below is not as complete as the one provided

in question 1; more, and larger, ‘gaps’ have been left for you to fill in...

1. Find out what millimetre-wave radar is and what “time of flight” means in that context.

2. Research how to convert a time-of-flight measurement in nanoseconds into a range in

page 5 of 8 pages

metres.

3. Load the data from the CSV file.

import numpy as np

raw_data = np.genfromtxt("radar_2021.csv", delimiter=",")

4. Find out how many samples are in the dataset using the command

(num_angles, num_ranges) = np.shape(raw_data)

5. Generate a vector angles containing the azimuth angle of each sample set (each row)

del_angle = np.radians( 0.144 ); # angles must be in radians!

angles = np.arange(num_angles)

angles = angles * del_angle

6. Create a vector ranges containing the ranges of each radar measurement using the

equation you found in part 2 of this question, with the time difference between

successive radar measurements set at 0.633 ns.

7. We want to be able to plot these data for visualisation. Start by associating each range

measurement with an angle using meshgrid():

rr, tt = np.meshgrid( ranges, angles )

zz = np.zeros(np.shape(raw_data))

8. Convert these data into Cartesian coordinates and visualise them as a scatter plot. Label

each axis (including the colour axis) and give the graph a title. Axis labels must include

units of measurement.

import matplotlib.pyplot as plt

# Convert to Cartesian, etc.

plt.scatter(...)

# Add labels, etc.

plt.show()

Hint: If your scatter points appear too large, consider passing the additional parameters

marker=’.’and s=1 to the plotting function.

9. The radar measurements have distinct returns, however there is noise introduced due to

thermal noise, imperfect power-matching, and manufacturing variations in component

fabrication. It is important to extract return signals that we are certain come from real

objects rather than being a result of noise. To do this we will use a CFAR detector.

Research CFAR detectors to find out what they are, and what the different variants are.

10. NumPy does not contain CFAR functions so we have supplied these (thanks Jacob!)

CFAR Detector Implementation Notes

a. We have supplied two functions CFAR1D() and CFAR2D()in the file CFAR.py

which implement 1-D and 2-D CFAR detectors.

b. The supplied detectors have variable window lengths at the edges of the data.

c. An example using CFAR2D() is shown below

x = raw_data;

hits2D, thresholds2D = ...

CFAR2D( x , pFA, nGuard , nTrain )

page 6 of 8 pages

11. Try applying a 2-D filter to the data using CFAR2D() and compare the outputs to those

from the 1-D filter CFAR1D(). How are they similar or different? Why is this so?

Question 3: Correction of sensor bias [26 marks]

This question relates to the ACFR Groundhog Robot, shown in Figure 1. Groundhog

navigates autonomously using data from an inertial navigation system (INS). The INS

contains a triaxial accelerometer and three gyroscopes oriented in the same directions as the

accelerometer axes. An accelerometer is a sensor that measures acceleration relative to an

inertial frame of reference; a triaxial accelerometer measures acceleration components in

three orthogonal directions.

Figure 1: The ACFR Groundhog robot.

The data for this question are from the triaxial accelerometer which is part of the INS

installed on Groundhog. The INS is fixed to the Groundhog frame with the accelerometer

oriented so the +x-axis points forwards, the +y-axis to the right and the +z-axis downwards.

INS data were acquired as the robot drove around a field. The data are in a file called acc-

2021.csv with the columns ordered as time (seconds), then the x, y and z components of

acceleration (m/s2).

1. In Python, load the data into a vector acc. Generate a time vector t from the data and

scale this appropriately.

2. Plot acc against t for each of the axes. What do each of the peaks represent? What is

the greatest acceleration experienced by the robot in each axis direction? Overall?

3. Acceleration data represent the time rate of change of velocity. If acceleration values are

integrated with respect to time, what do we end up with? What if we do it again?

page 7 of 8 pages

Integration of a signal involves summing the area under the curve. In sampling this

signal, we are sampling a continuous (analogue) signal—acceleration—at a set of

discrete times. There is no information available about the signal except at the discrete

times when the signal was sampled. In the time between samples, it could be assumed

either that the signal is constant or that it changed linearly. Assuming a linear change is

arguably the better assumption. Under this assumption the sampled signal can be

integrated by summing the trapezoidal areas defined by connecting adjacent sample

points with straight lines. This is the trapezoidal rule for numerical integration.

4. Calculate the cumulative sum of the individual accelerations using the trapezoidal rule

to generate a vector of velocity measured in metres per second.

Note: Python has a cumulative trapezoidal sum function called cumtrapz(). To use

cumtrapz(), add from scipy.integrate import cumtrapz at the top of the

file, and use the syntax cumtrapz(x,initial=0), where x is the variable being

passed to the function and initial=0 ensures an output matrix of the same size.

5. Do the same for the velocities to generate a vector of position in metres.

6. Plot the cumulative sums against time. When does the robot start moving? When does it

stop? What is its speed whilst moving?

7. Turn on the grid for the plot using the Python command plt.grid(). Can you see any

long-term trends that should not be there? When the robot is not moving, its position

should be constant, and its velocity should be zero; neither the position nor the velocity

cumulative sums should change. Are the position and velocity lines exactly horizontal

during periods when the robot was stationary? Can you see a trend where (supposedly)

horizontal lines have a non-zero slope?

A trend where the data ‘drifts’ in an unexpected way often comes from integrating a

signal bias. The bias is noticeable here when we expect a horizontal line (constant

position), but the bias in fact influences the entire dataset. An INS has a bias that

changes for several reasons, mostly related to temperature. Because of this, we need to

calculate and remove this bias for each dataset. A good way to identify the bias is to

leave the INS stationary for a short time, usually at the start or at both the start and

finish of data collection. If accelerometer bias is present, the average of the samples

collected will be non-zero2.

8. The samples between 4400 and 4497 inclusive were acquired while the robot was

stationary. Plot a histogram of these values. Python can plot a histogram using the

plt.hist() function; read the Matplotlib online documentation on Histograms to find

out how to use it.

Note: A subset of a vector can be selected in Python using subset =

vector[start:end]. For the current example, since the index starts from zero in

Python, use stationary_data = scaled_data[4399:4496].

9. Does the average value of the histogram look like it is zero? Find the actual average

value using Python’s np.mean() function and store in the variable bias.

10. This calculated mean value indicates that on average, each sample is moved away from

zero by this amount. Make a new variable acc_unbiased where you subtract bias

2 Provided that all noise present is zero-mean.

page 8 of 8 pages

from the biased data acc.

11. Repeat steps 3-8 using the new vector acc_unbiased. Generate a set of plots

containing the cumulative sums of acc and acc_unbiased and comment on the

difference.

Question 4: GPS data [26 marks]

Groundhog was also receiving satellite navigation data during the same mission as in part 3.

We wish to compare the inertial measurements to the GPS tracks in this question.

1. The GPS data were collected from a differential GPS system. Research differential GPS

and compare it to regular GPS.

2. Load the GPS data from the file GPS_2021.csv. From the latitudes and longitudes in

the file, where in the world did the measurements take place?

3. The GPS unit provides navigation solutions as coordinates in the latitude, longitude, and

altitude frame of reference. These coordinates are difficult to work with for most

ground-based robotic navigation tasks, so an alternative coordinate system—Universal

Transverse Mercator, or UTM—is used. Research this coordinate system and find the

transform from LatLonAlt to UTM. There are several different UTM spatial reference

systems in use; you should preferentially use Map Grid of Australia 2020 (MGA2020)

or WGS 84.

4. Apply this transform to the GPS data provided.

5. Plot the data in the UTM frame of reference and compare it to the results for the

position found in Question 3. How are they similar/different? What are some ways that

these two sets of data could be used together?

Revision history

A DCR 04/05/2021 Clarified “Assignment Submission” section

— DCR 25/04/2021 First release

Revision Made by Date Detail

欢迎咨询51作业君