Development

The Team has been hard at work conceptualising, testing and refining the product. Highlights of the process will be posted here.

Part 1. 15th June ’19:

Coding started on mid-June 2019 and we focused on isolating dominant frequencies from recorded audio samples from computer mics. The Fast Fourier Transform (FFT) function was built in to the software and was extensively used in this particular process.

Several tests were done on determining the dominant frequencies of certain sounds.

Dominant frequency for banging table with hand is approximately 78.56Hz.
Dominant frequency of hitting a steel cup is approximately 3401.9 Hz.

As observed, noise was a major issue in these readings. Another analysis was done with the same steel cup.

Dominant frequency of 4162.8Hz for a different sample of the same steel cup.

Part 2. 20th June ’19:

It was hypothesized that the noise in the environment could have played a role in the inconsistent results. It should also be noted that infrasounds (the sounds that we intend to analyse in the prototype) are below 20Hz and require specialised microphones. The decision was made to purchase and experiment with SenseHAT. SenseHAT is an add on board that connects to all the GPIO pins on the Raspberry Pi and contains several sensors, such as a gyroscope and an accelerometer. It even has a humidity sensor but we will not be using that.

Top View of Raspberry Pi with SenseHat Module installed.

Since the Sense HAT module already used up all of the available GPIO pins, it was imperative that we made use of the USB I/O interface. We then planned to explore on USB connectors that adapt to compatible pins with the infrasound mics. A keyboard and HDMI cable was obtained for the purposes of programming within the Raspberry Pi itself.

Part 3. 25th June ’19:

Leveraging on the Audio Toolbox by MATLAB, we created a script which could decouple a mixture of sounds into their individual components. This was done through binary masking, where we trained the code to pick out noise signatures of the Short Time Fourier Transformation (STFT) of our desired component, and applied the mask to the mixed audio. 

We were successful in retrieving an audible, albeit muffled, sentence spoken by one person out of an audio containing 2 speakers and noisy background. 

Part 4. 28th June’ 19:

However, on further advice by our mentor, Prof Benoit, we decided not to use audio and infrasound triggers in our sensor, as they are not reliable predictors of an upcoming eruption. They also require vast computational and electrical power, and will not be a good fit for our sensor – our sensor aims to be autonomous and affordable.

We will be harnessing the SenseHat module on top of the Raspberry Pi to detect changes in motion. The code was written in Python to allow the system to distinguish between types of motion, and output different patterns on an LED array. The different movements include: 1. Sliding but no tilting; 2. Tilting in different directions but no sliding; 3. Sliding and tilting in different directions.

Before an eruption, magma will accumulate in the chambers around the volcano. The area under ground will swell and elevate slightly. The sensor can be orientated to face the volcano. As a result, we can reduce the likelihood of false positives by triggering the alarm only when the tilt is detected.

Different volcanoes have large variations in seismic activity, and it is rather difficult to collect data and test the algorithm. Professor Benoit hence advised us to come up with a fixed trigger which mimics the activity of the volcano: comprising of tilt and relatively high frequency vibration. We decided to use a manually tilted Chaldni plate and mechanical vibrator as our test bed.

Part 5. 1st July ’19:

We used the SenseHat module and Raspberry Pi modules to detect vibrations on the Chaldni plate. The frequency generator allowed us to precisely control the frequency of vibration and we introduced tilt by placing a items on one side of the stand and using trigonometry to measure angle.

For the code, we introduced a Python list which stored data of acceleration in X, Y and Z directions; and yaw, roll, and pitch, then plotted the points. Additionally, we noticed that it was important to optimise the sampling rate. Too low a rate gives incomplete data, while overly-high rates increases the difficulty of data handling.

Pythagorean Sum of accelerometer data from Washing Machine

The raw data was initially recorded separately in the X, Y and Z directions. To be able to have a more representative visualization of the data, we used the Pythagorean sum. It is essentially the square-root of the sum of squares of the data in all 3 directions, similar to how the length of a diagonal of a right-angled triangle is found with the length of two other sides.

Part 6. 10th July ’19: Data collection, Data Cleaning, Machine Learning

Being more familiar with the parameters and functions of Sense HAT module and Raspberry Pi, we wanted to use it to detect seismic motion from a slightly more chaotic source – bomb bags. With the initial features of a bomb bag explosion unknown, we aim to train a Machine Learning (ML) algorithm that is capable of detecting if a disturbance is caused by a bomb bag explosion or other sources (such as us manually dropping the sensor). 

The bomb bags that we used for our project

However, to train our ML algorithm, we first had to come up with enough data sets. And hence began the long and painful (to the ears) process of data collection. Before we dive into the actual data collection, however, let us discuss about the setup.

As shown in the photo above, the setup involved a Raspberry Pi (with SenseHat mounted) taped to a plate. The Raspberry Pi was then connected wirelessly to our computer via WiFi. This allowed us to control Raspberry Pi on our computer through a VNC viewer. Under the plate, we had placed a bomb bag which would provide the disturbance needed once it exploded. The setup itself looked simple but the coding that was needed turned out to be quite challenging indeed (more on that at the end of this chapter).

With the Python code and physical setup prepared, we were ready for data collection: Bob presses the button to run the Python script while Zheng Hao bursts the bomb bag  placed under the plate. The bomb bag would then expand rapidly, making the SenseHat experience high amounts of tilt in a short span of time. Then BOOM! The bomb bag explodes and the plate, along with the sensor, drops back to the table.

The video above demonstrates the collection of one positive data set. We repeated this procedure 81 times to get 81 sets of positive data. Following that, we collected 81 sets of null data by manually raising the plate and dropping it.

Being the curious reader, I’m sure the question you have right now is: what type of data did we collect?

Well, the data that we collected was the linear/angular acceleration values of the sensor for/about the:

  1. X-axis
  2. Y-axis
  3. Z-axis
  4. Pitch-axis
  5. Roll-axis
  6. Yaw-axis

Take note that we were only concerned with the values of these characteristics during the period when the sensor drops (ie. right after the bomb bag explodes).

To collect each set of data, the Python script runs for 3 seconds with 100 data points taken. This means that each data set consists of 100 values each of X-acceleration,Y-acceleration, Z-acceleration, pitch-acceleration, roll-acceleration and yaw-acceleration. We chose such a high sampling frequency because we needed a large number of data points for our Machine Learning model to be accurate.

However, we soon discovered that a run time of 3 seconds presented a fresh challenge: the bomb bag took longer than 3 seconds to explode. This was obviously a problem because we were interested in the acceleration values when the sensor dropped back to the table (00:57 in the above video). Two solutions came to mind but they each had their flaw.

  1. Run the python script just before the bomb bag exploded. Unfortunately, this did not work because the bomb bag exploded in a very unpredictable fashion (even after it was fully inflated, it took a varying amount of time before exploding).
  2. Increase the run time from 3 seconds to 2 minutes. This approach guaranteed that we would capture the falling of the sensor (since almost every bomb bag exploded within 2 minutes). However, this meant that we would have thousands of data points for each data set which would take up lots of memory space in Rasp Pi. Furthermore, given that we had 81 data sets, the amount of data collected would have been massive.

And hence, we had a huge challenge facing us…

Main challenges:

Getting our sensor to only record data when the sensor is dropping given that: 1. explosion of bomb bag is unexpected and 2. we could only store 3 seconds’ worth of data points.

Solutions:

After brainstorming and analyzing our preliminary data, we discovered that when the sensor dropped, there was a sharp decrease in acceleration-value for the Z-axis since the contact force between the plate and SenseHat dropped (note that SenseHat measures acceleration by using the amount of G-force acting along that particular axis). Knowing this, we wrote a code that would continuously run 3-seconds loop. At the end of the 3-seconds, we set the condition that the data collected would only be recorded into a CSV file if the difference between the maximum and minimum Z-acceleration values was greater than 0.6. Otherwise, the data would be deleted and the loop would run again.

By doing so, this allowed us to capture the data concerning the period where the Sensor was dropping without having excessive data points (:


Spyder would only “get” out of the loop and proceed to record the data in a CSV file if the max range of Z-acceleration exceeded 0.6

 

Inputting the measured data into a CSV file

Data Cleaning

After collecting the required data (81 positive sets and 81 negative sets), we moved on to data cleaning. Remember that each data set had 100 data points for each characteristic (X, Y, Z accelerations, Pitch, Roll, Yaw) so that’s 600 values per data set (and we had 81 positive data sets). The following screenshot of part of our CSV file shows just a fraction of the total raw data that we collected..The raw data we collected.(Note that the each column represents a different data point while each row represents a particular characteristic of a particular data set)

It was therefore crucial that we extracted the key information from these values to be used for our Machine Learning algorithms. After due research, we concluded that the relevant features would be range of accelerations (max acceleration – min acceleration), standard deviations of accelerations and mean acceleration. So, from the initial 600 values (6 characteristics x 100 data points) per data set, we condensed it into 18 values (6 characteristics x 3 features). Looking at the chunks of data above, you would probably be wondering how could we split them nicely into the 18 values per data set? The solution is simple.

We knew that there was a pattern behind how the data was tabulated in the CSV file: Each data set took up 6 rows. The first row would be the pitch values, second to fourth rows the Z, X and Y acceleration values and the last two rows the roll and yaw values. This pattern repeats itself every 6 rows. Hence, we wrote a code that did the following: Take the row number and find its remainder after dividing by 6. Using this remainder, Raspberry Pi could then know which characteristic that row of code corresponds to. Then, we used methods contained within the NumPy library to find the range, standard deviation and mean. Finally, we added all these data into a new CSV file.

The code we wrote to neatly arrange the massive amount of data collected before putting them into a new CSV file. This CSV file is then used in our ML algorithm.  

This CSV file- with the data neatly sorted out- was then used to train and test our Machine Learning algorithm.

Cleaned data sets with the relevant features extracted out to be used in our Machine Learning algorithm

Machine Learning

With the cleaned data sets, we then had enough “ammunition” to proceed with our task of training and testing our ML algorithms. After some research, we concluded that the most appropriate algorithm would be Random Forest which utilizes a number of decision trees to come up with a decision (in this case, whether the drop was caused by the bomb bag explosion or manual drop). Each decision tree in the forest would use different features to make an educated guess.

An example of a decision tree (of course, the one used by our algorithms are a lot more complex)
Source: https://chunml.github.io/ChunML.github.io/tutorial/Decision-Tree/

Each decision tree has different nodes with different criteria. Ultimately, the algorithm would pick the decision chosen by the majority of the trees.

Even after selecting an algorithm, we had to spend time fine-tuning the relevant hyper-parameters (number of trees, depth of each tree) before we were satisfied that the model was good. Take note that this fine-tuning of hyper-parameters was different from training our model.

Then we used the acquired cleaned data to train our model. Since we already decided on using a particular model, there was no need for us to compare the accuracy of different models and hence no need to use data for validation purposes. Therefore, we had the luxury of using 80% of the data for training and the remaining 20% to test the accuracy of our model (as compared to using 60%, 20% and 20% for training, evaluation and testing respectively in conventional methods). Using the test data sets, we evaluated our model and we are proud to announce that the model was able to correctly predict the disturbance as being caused by bomb bags explosion 96% of the time i.e. an accuracy of 96%.

Our Machine Learning Model

The picture above shows the code to train our ML model. The brief outline of the code is as follows:

  1. Importing the relevant methods (RandomForestClassifier, cross_val_score etc.) from the sklearn library and the CSV file
  2. Splitting the data contained in the CSV file randomly into training sets and test sets.
  3. Having our model predict the y-values (0 if manual and 1 if bomb bags)
  4. Comparing the predicted values against the actual values to determine the accuracy of our model.

With an accuracy of 96%, we were satisfied that our model was optimal. Lastly, what remained to be done was adding a path that allowed Rasp Pi to use back this model to predict the source of subsequent disturbances.

In the code above, Rasp Pi can quickly predict the source of the disturbances by using the data encapsulated in df2 and the model we just created, gs. That brings us to the end of our adventures with bomb bags and ML. (We will meet again in part 9)

Part 7. 17th July ’19: Prototype 1 of SensorX

After we successfully trained our machine learning algorithm to predict bomb bags explosion with 96% accuracy, we started work on making a complete package that utilizes both the tilt and the vibrations detected by the SenseHat to predict whether there is a volcanic eruption.

Schematic diagram of full package

As can be seen from the diagram above, the Arduino board acts as the micro-controller that controls the separate electronic components (buzzer,led and tilt sensor).

Brief explanation on the working principle of our full sensing package

The video above offers a brief explanation on what triggers the LED: tilt (measured via tilt sensor) and vibrations(measured using Raspberry Pi)

  1. Tilt. To make the tilt sensor (00:45), we pasted a sheet of aluminium foil onto the inner surface of a plastic wine glass. Inside the dome of the wine glass, there is a metal washer that is again wrapped in aluminium foil and suspended by a metal wire. The other end of this wire is connected to a breadboard. We then used another wire, with one end connected to the breadboard and the other end to the aluminium foil that is pasted on the wine glass. When the tilt sensor is on a level surface, the metal washer hangs freely and it does not touch the aluminum surface. In this case, the circuit is opened. However, when the sensor is on an inclined surface, the washer swings and touches the aluminium surface.This causes the circuit to be closed.
  2. Vibrations. Using the SenseHat mounted onto the Raspberry Pi, we are able to measure vibrations.

Condition 1 + Condition  2 = LED lit up!

Using a code that we wrote, the Arduino board is able to simultaneously recognize that the circuit is closed (condition 1) and that the Raspberry Pi is experiencing vibrations (condition 2) and hence command the LED to light up.

Demonstration of our sensing package

Main Challenges:

  1. Getting Arduino to “receive information on whether there is vibrations” from Raspberry Pi (recall that the Raspberry Pi has to recognize that there is vibrations before conveying this information to Arduino). This is especially difficult because Arduino and Raspberry Pi operate on entirely different platforms (C++ and Python respectively).
  2. Getting Arduino to recognize that the circuit is closed.

Solutions:

Highlighted is the code that allows Arduino to “know from Rasp Pi that there is vibrations


1. Looking at the highlighted lines in the first set of code (which is written for Raspberry Pi), using the “Serial’ function, Raspberry Pi will encode a value of “s” into the serial port “ttyACM0” if it detects a vibration (defined as the range of Z-acceleration being greater than 0.6). Then, in the next set of code (written for Arduino), Arduino will execute a simple if-else statement: If the value encapsulated in serial port “ttyACM0” is “s”, then it will run the subsequent digitalRead() and digitalWrite() functions. (See next paragraph for more details). Otherwise it will just run digitalWrite(led,LOW) i.e. LED is off

Highlighted is the code to recognize that there is tilt

2. As can be seen from our code above, we utilized the digitalRead and digitalRead functions to trigger the lighting of the LED. Using the digitalRead function, it will return a value of “HIGH” if there is a voltage detected at pin 2 and “LOW” if there is none. Naming this value as “state”, we then pass it as a parameter into digitalWrite, which then causes the LED to light up (if state=HIGH) or remain off(if state=LOW).

Part 8. 19th July ’19: Attaching a bluetooth module to SensorX

After completing the 1st prototype of SensorX which utilizes both tilt and vibrations to detect an oncoming volcanic eruption, we started brainstorming on how we could send out a warning to more people in the vicinity in the shortest amount of time. After some research, we concluded that using Bluetooth as a signal transmitter is the best solution primarily because a Bluetooth module consumes a low amount of power. Furthermore, by tapping onto the Bluetooth Low Energy technology, we were further able to reduce power consumption which equates to less human maintenance .

The video below shows the working principle of our Bluetooth module.

In the video, we see that when the circuit is closed (the white wire is inserted into the breadboard), the led lights up and the Bluetooth transmitter sends a message to the phone. For the phone to receive the message, we downloaded the app “Bluetooth Terminal HC05” from the Google Play Store.  The picture below is the code we used for Arduino to command the transmitter.

Arduino code for Bluetooth transmission

The flow of the code is as follow:

  1. Importing the BTserial library.
  2. Designating the Arduino pins to be used for receiving(10) and transmitting(11).
  3. (FOR LIGHTING UP OF LED) Having Arduino check the state of pin 2 (whether there is voltage across pin 2) using the method digitalRead(). If there is, Arduino will pass a voltage through pin 13 using the method digitalWrite(). This voltage causes the LED to light up.
  4. (FOR TRANSMISSION OF BLUETOOH SIGNAL) If Arduino detects a voltage across pin 2, it would send out a message via the transmitter.

Moving forward, we plan to incorporate the Bluetooth module into our SensorX package such that  when both tilt and the right vibrations pattern are detected, the transmitter would send out a warning message to all smartphones in the vicinity, warning them on an imminent volcanic eruption. However, since the range of Bluetooth transmission is short- we tested ours and found out that it had a max range of 9m- we plan to use Bluetooth routers to increase the transmission range and hence increase the warning distance.

Part 9. 18th July ’19: Prototype 2 of SensorX (combining Prototype 1 with Machine-Learning functionalities)

After discussion with Dr Ho, we concluded that we needed to build a model volcano as a test bed for SensorX. And hence, we began to brainstorm how to go about building one, keeping in mind that our model has to produce both tilt and vibrations, just like an actual volcano. Also, our volcano model must be able to produce both positive and negative data sets because our aim was to have the ML algorithm within SensorX distinguish between vibrations from an actual volcanic source and other sources (perhaps the movements of goats in the area). Eventually, we settled on a model that fulfilled all these requirements. We will now explain how our model works.

Design Setup :

We attached a piece of plywood onto a hinge which was then attached to a large piece of acrylic plastic(the base of our model). The hinge allowed the plywood to rotate freely. Then we inserted a hydration bag underneath the plywood and taped a mini-motor onto the plywood. Lastly, we taped SensorX onto the plywood.

Setup to simulate volcanic activity

Design principle

When we pump air into the hydration bag, the expansion of the bag causes SensorX to be tilted. For the positive vibration source, we use the taped mini DC motor to produce vibrations on the plywood. For the negative vibration source, we manually shook the plywood (causing SensorX to be shaken together with it).

Goal

We aim to have SensorX capable of distinguishing the positive vibration source (DC motor) from the negative sources (manual shaking of plywood) via Machine Learning. Hence, only when SensorX detects a vibration from the DC motor (and not just any vibrations!) and is tilted beyond a certain angle will it be triggered.

Important note: Prototype 2 is significantly different from Prototype 1 in that only when the vibrations signals detected by Rasp Pi matches that from a positive source will the vibrations condition be fulfilled. Recall that previously in Prototype 1 (mentioned in detail in part 7), as long as the amplitude of vibrations was above 0.6, the vibrations condition was deemed to have been fulfilled. You can think of Prototype 2 as an upgrade of Prototype 1 with Machine Learning capabilities.

With the setup prepared, the next step for us was to collect data to train/test our ML model, just like what we did for the bomb bags. The video below illustrates how we went about collecting one set of positive data (you can hear the noise generated by the vibrations of the DC motor in the background of the video).

With the data sets (196 positive and 196 negative sets) collected, we then proceeded with data cleaning. Just like the bomb bags data (part 6), we collected 100 data points for each of the 6 features (X, Y, Z acceleration, roll,pitch and yaw).  Since the process and code used were similar, we would not repeat them here. 

With the cleaned data, we were then ready to train our ML model. Again, we found Random Forest to be the most suitable algorithm. With a suitable model chosen, we then used 80% of the data to train the model and the remaining 20% to test it. This time round, our model achieved an accuracy of 94%, which is excellent. We were delighted with our accomplishment.

Since our mentor, Mr Kanesh, happened to be in the SPMS building, we decided to seek his feedback. He told us that our ML algorithm was great but a drawback of our model volcano was that by using a DC motor as our vibration source, we were only able to get vibrations of a very narrow range of frequencies. Given that an actual volcano vibrates over a wide range of frequencies (1 to 25 Hz), our model lacks realism.

Hence, after another round of brainstorming, we decided that it would be much better to use a mechanical wave driver as a vibration source since it is able to produce vibrations of varying frequencies and amplitudes.

Setup plus SensorX (inside box)

In the photo above, the mechanical wave driver(located on the right of the setup) provides the vibration source by oscillating back and forth. Through the sine wave generator (grey box next to the driver), we are able to control the frequency and amplitude of the vibrations and adjust these parameters to fit those of an actual volcano.

Being ever the observant reader, I am sure you would have noticed that our setup looks entirely different from the original one( using the DC motor). That’s because we changed our setup to accommodate the mechanical driver. To allow SensorX to vibrate with minimum impedance, we taped it onto a piece of plywood which was then placed on 2 PVC tubes which were allowed to roll. The plywood was attached to the mechanical driver through a thin metal sheet.

We are sure that a discerning and intellectually curious reader (the fact that you’re here proves that 😉 ) like you would know that changing the setup meant that the data sets that we collected previously are now obsolete. Which means….another round of data collection?

Ah yes, you’re right, we will have to do that. But that’s a story for another time! Cya in the next part!

Part 10. 19th July ’19: Data collection (with mechanical wave driver as vibration source) and feature extraction

In the last post, we mentioned that after consultation with Mr Kanesh, we decided that we would use the mechanical wave driver instead of the original DC motor to provide the vibrations. That meant that we had to recollect our data. Hence, the first part of this post is dedicated to detail how we went about doing so. Then we would move on to discuss how we extracted the key features from our raw data.

Data Collection

This time round, as we expected this volcano setup to be our final one, we decided to collect a large number of data sets in order for our ML model to achieve a high degree of accuracy. Hence, we went about collecting 800 sets of positive data and 800 sets of negative data. Here, you might ask how did we produce the conditions for our positive and negative data sets? To produce the positive results, we used the mechanical driver as the vibration source, varying the frequency (8-12 Hz) and amplitude while collecting the 800 data sets. We varied the frequency and amplitude of vibrations because we wanted to replicate actual volcanic conditions as much as possible. To produce the negative results, we manually shook the plywood that SensorX was resting on.

The video below shows how we went about collecting positive data.

Importantly, for this round of data collection we collected data only for the 3 axes of linear acceleration. We chose to neglect the 3 angular accelerations (pitch, roll, yaw) because we wanted to reduce the computational strain on Rasp Pi. Also, after analysing data from previous rounds, we discovered that there was little information contained in the angular acceleration values.

Since the code that we used for data collection was largely similar to previous rounds (only difference was that we removed the lines telling Rasp Pi to record the angular accelerations), we would not discuss about them here.

Feature extraction

With the collected raw data, the next step for us was to extract the key features. These were the steps we took:

  1. Apply filtering and windowing to raw data.
  2. Use Fast Fourier Transform (FFT) to produce 3 amplitude-frequency graphs of each data set. (1 each for X,Y, and Z acceleration)
  3. Extract dominant frequency for each graph.

“Chottomatte!” you say, why did we choose to extract the dominant frequency? Recall that in our previous 2 rounds of machine learning, the features that we had used were the standard deviation, range and mean values of the characteristics. But now we were going to use the dominant frequency?

The reason was simple: After reading up on volcanoes (and coincidentally we happened to have a volcano mini-expert in our group), we found out that most volcanoes erupted within a small vibration frequency range that is within the larger range of 1-25Hz. And most interestingly, we also discovered that vibrations signals are one of the most prominent tell-tale signs of an imminent eruption. In other words, if we could have SensorX analyse whether the vibrations it detected was from this small frequency range, then it could potentially predict whether there would be an eruption with a high degree of accuracy. And we planned for SensorX to be able to do this via a Machine Learning model.

Of course, to successfully train our ML model to predict on the basis of vibrations frequency, we would have to supply it with relevant data. Hence, we decided to provide it with the dominant frequencies of each data set.

“Ahhhh…..” you say, so that’s why we chose to focus on extracting the dominant frequency.

With a clearer understanding on the importance of dominant frequency, let us explain in detail the 3 steps we took for feature extraction.

Apply filtering and windowing to raw data

Low-pass filter:

In our data sets, we varied the frequency produced by the sine wave generator to be from 8-12Hz because the frequency generated by the seismic activity of a volcano averages at 10Hz. Since we knew that any high frequency signals (>20Hz) was noise, we applied Butterworth filter such that we only record signals that are less than 20Hz. We did this to prevent high-frequency noise from appearing as peaks and distorting our actual data.

Code for filtering our raw data

The picture above shows the code we used to filter our raw data. Firstly, we imported the signal toolbox from the scipy library (line is not shown in picture). Then we initialized the filter object, which we named as butter, in the first line of code. Next, we passed our raw data, X, along with butter into the sosfilt method. The method returned us the value filtered, which is an array containing our filtered data.

To visualize the effects of filtering, we printed two diagrams. One with FFT and windowing applied and the other with FFT, windowing and filtering applied. Check out the results below.

Effects of filtering

Beyond 20Hz, the amplitude of the red graph drops sharply. That was because we set the threshold frequency as 20Hz. Since, we knew that the vibrations experienced by SensorX was definitely lower than 20Hz, filtering provided us with cleaned data that was more accurate. 

Windowing:

If we could do a Fast Fourier Transformation(FFT) on a signal that is infinite, the results would tell us perfectly the frequencies (and their corresponding amplitudes) that are present in the signal. However, in reality, the data points collected are finite (100 data points in our case) and hence FFT resolves this problem by assuming that our 100 data points are infinitely repeated. This approach does not present any problem when the frequencies of the constituent waves (also known as input frequencies) are exact integer multiples of the frequency resolution. (Note that frequency resolution is sampling frequency divided by the number of sampling points. It is the smallest change in frequency that the FFT can detect. Hence, it is also the step interval between adjacent bins on the frequency domain graph produced after FFT.) The picture below complements our explanation. In this example, the input frequency is 125kHz, the sampling frequency is 1000kHz and frequency resolution is 62.5kHz.

Basics of FFT. Source:https://training.ti.com/ti-precision-labs-adcs-fast-fourier-transforms-ffts-and-windowing

As mentioned above, repeating our sample points infinitely is not problematic if our input frequencies are integer multiples of resolution frequency. But what if they are not? (In fact, it is far more likely for our input frequencies to not be integer multiples).

In that case, we will have the problem of waves (each with one single frequency) occupying multiple frequency bins. This phenomenon is called a spectral leakage. Recall that FFT repeats each set of data points (called a record) infinitely to produce an infinite time-domain signal. If the input frequencies are all integer multiples, then each record connects smoothly to the next. But if they are not, then there will be  discontinuities at the end of each record. This causes spectral leakage which is problematic because it makes it hard to differentiate the signal from noise. Hence, windowing reduces the magnitude of this problem by smoothing out these discontinuities. It does so by gradually reducing the amplitude of the signal to 0 at the start and end of each record. The picture below shows the effects of windowing.

Effects of windowing. Source:https://training.ti.com/ti-precision-labs-adcs-fast-fourier-transforms-ffts-and-windowing

With the importance of windowing explained, we will now discuss how we performed windowing on our data.

After some research, we narrowed down the available types of windows to 3: blackman, hanning and hamming. Each window has its own set of characteristics which we will not cover here. To determine which window is best suited for our data, we did the following:

  1. Of the 800 sets of positive data, we took out a test set.
  2. We performed the 3 types of windowing on the test set.
  3. We analysed the graphs produced by the different windows.

The following picture is the python code we used to perform windowing.

The flow of the code is as follow:

  1. Import blackman, hamming and hann methods from scipy.signal toolbox (not shown in picture)
  2. Using the above methods, initialize the windowing objects (bm,hn and hm)
  3. Pass the windowed data (bm* X, hn* X and hm* X ) into the fft method.
  4. Plot the 3 FFT graphs (each with a specific type of windowing performed) to compare the effectiveness of the windows

 

 

The results produced by the 3 windows are shown below. You can see that by comparing the FFT (blue) graph against the one with windowing performed (red), the amplitude of non-dominant frequencies is generally lower when windowing is performed. Again, by comparing the 3 red graphs, we could see that the hamming window was the most effective as it reduced the amplitude of non-dominant frequencies by the most. Hence, we decided to use the hamming window for the remaining data sets.

Effects of the different windows

Hence, we proceeded to use the hamming window for our remaining data sets.

Note: The observant reader in you would have realized that the effects of filtering and windowing were minimal. Indeed that is true. However, that is because the vibrations data that we collected were already “clean”. “Clean” in the sense that there were very little noise since we conducted the data collection in the SPMS physics lab, free from external disturbances. In the real world however, filtering and windowing are common industry practice due to the prevalence of noise. 

Also, you might ask: What is the method semilogy? What does it do?

Well, semilogy creates a plot using a base 10 logarithmic scale for the y-axis and a linear scale for the x-axis. This is crucial as we wanted to visualize the differences between the graphs – otherwise it will be difficult to notice the differences with a linear scale.

After filtering and windowing, the next step was to pick our dominant frequency.  This was a relatively easy step as we simply ran a Python code that went through the array of amplitude values and returned us the highest amplitude, along with its corresponding frequency.

With the 1600 (800 positive and 800 negative) data sets prepared and backed by our technical knowledge of filtering and windowing, we were ready for Machine Learning.

Part 11. 22nd July ’19: Prototype 3 (Fitting our extracted data into the Machine Learning algorithm)

After we successfully collected the 1600 data sets, the next step was to use them to train our ML model. When completed, this would be our Prototype 3. Before we move on, here’s a brief recap of our 3 prototypes.

Prototype 1: Uses both tilt and vibration to get triggered. But it’s simpler because long as the magnitude of acceleration is greater than 0.6, the vibrations’ condition is fulfilled.

Prototype 2: Also, it uses both tilt and vibration to get triggered. It’s more complicated than Prototype 1 because it uses a trained Machine Learning model to determine if the vibrations pattern resemble those on a volcano before an eruption.

Prototype 3: Again uses both tilt and vibration to get triggered. It is however the most complex because not only does it utilize a ML model for analysis of vibrations pattern, the feature that we fed into the model (and the feature it looks out for) is the dominant frequency of the vibrations. To understand why we chose to use the dominant frequency, please refer to previous post.

The following is the hardware of Prototype 3. (Note that the hardware is exactly the same as Prototype 2)

The training of our model was split into 2 parts. Firstly, we had to feed all the data into our model before hypertuning the parameters to find the most suitable model.

Feeding data into ML model

In this section, we will discuss how we used methods contained within the numpy and pandas library to input our data into the model.

Importing the relevant files and functions

The flow of the above code is as follows:

  1. We import the relevant libraries pandas and numpy to allow us to access their inbuilt functions. Unlike its animal namesake, the pandas (python data analysis library) library is a powerful tool for data manipulation and analysis while numpy is useful when we are dealing with large volumes of data.
  2. We retrieved the positive data sets (data_collection_18.csv) from our computer. Thanks to the read_csv model contained within numpy, our data was conveniently stored in the DataFrame format. Wait…before we move to the next point. Here’s a quick understanding check: When we refer to the positive data sets, what are they? Answer: They are the 100 sample values of each acceleration (Z,X or Y) for each data set. Since there are 800 positive data sets and 3 directions of acceleration, we have 2400 lines of values (with 100 values in each line).
  3. Then, we imported the signal, fft and hamming functions from the scipy library. Recall: We used the hamming window because we found that it is the most effective at reducing spectral leakage for our data sets. We performed butterworth low-pass filter to prevent high frequency noise from appearing as peaks on our FFT.
  4. Lastly, using the linspace and asarray functions, we created an array containing all the frequencies of the bins of our FFT. If you are observant, you would have realized that we limit the maximum frequency of our FFT to be half of the sampling frequency (100/1.5Hz). The maximum frequency of a complete FFT is the sampling frequency. So why did we essentially remove the second half of our FFT? Because the second half (aka the alias region) is essentially a duplicate of the first and it is redundant to include it. Hence, we simplified our code by removing it.

The next step was then to retrieve our data from the DataFrame and sort them out such that we can feed it to our model.Performing filtering, windowing and then FFT

The flow of the above code is as follows:

  1. We first created an empty array a with 1 row and 3 columns as well as an empty list called rows_list.  We will explain the purpose of doing so shortly.
  2. We created a loop that iterated the range from i=0 to i=2399. If you have no idea why we set the endpoint as 2399, please re-read point 2 from the above paragraph.
  3. Using the iloc[i] function, we retrieved the value contained within the ith row of our DataFrame, df_pos.
  4. We then performed filtering on it using sosfilt. Then, we performed hamming (windowing) on it and lastly FFT. With the asarray function, we converted it from a list data type into an array. This is important because an array is able to have many functions performed on it, but a list is far less versatile. At the end, what we got was an array containing all the amplitudes of the signal at different frequencies. 
  5. Using argmax, we were able to find out the index position of the maximum amplitude contained within the array. BUT since we were interested in the dominant frequency, we had to find it from the xf array (which contains all the frequencies of the bins in our FFT) using the index position.
  6. We added the value of the dominant frequency into rows_list. When the list reached 3 values (the 3 dominant frequencies for the X,Y and Z accelerations from 1 data set), we (and here’s the ingenious part) stacked them vertically onto the previous sets of 3 values using the vstack function. This meant that ultimately, we had an array with 3 columns and 800 rows. Each row contained the 3 dominant frequencies of each data set and each column contained all the X/Y/Z dominant frequencies. The picture below is a visual representation of our a array. With our data sorted out in such a neat way, feeding our model became a much simpler task.

Visual Representation of a array

7. Then the last line “processed_df_pos = pd.DataFrame({‘Z’:a[:,0], ‘X’:a[:,1], ‘Y’:a[:,2]}) ” created a DataFrame with 3 columns: Z, X and Y.  All the Z dominant freqs were in the first column, all the X dominant freqs were in the second column and the Y dominant freqs in the third column. This step was crucial as a ML model could only take in data when they are in a DataFrame and not in an array.

8. We added in another column called “event” with a value of 1 (the corresponding value for the negative data sets is 0) to let our model differentiate between the positive and negative data sets.

We then repeated these exact steps once more, but this time for our negative data sets. Once this was done, we combined both DataFrames into one.

After reading up on relevant materials online, we narrowed down our choice of Machine Learning algorithms to 3: Decision Trees, Logistic Regression and Support Vector Machines. However, we did not know which was the best among them. Also, each algorithm came with a set of hyperparameters that had to be defined by the user. Again, we were unsure of the optimal values.

Main Challenge

We did not know which ML algorithm would yield the most accurate results. Neither did we know the optimal values for our hyperparameters.

Solution

We performed a RandomizedSearch on all 3 ML algorithms. GridSearch essentially allowed us to test a user-defined range of hyperparameter values to see which set of values were the best. Then, using the best hyperparameter values for each algorithm, we could then determine which was the best algorithm out of the three.

The remainder of this post would be dedicated to detailing how we went about this endeavor. 

 

Fine-tuning Decision Trees (Random Forest) model:

Fine-tuning Random Forest model

The flow of our code is as follows:

  1. First, we imported the RandomForestClassifier ML algorithm and RandomizedSearchCV optimization algorithm from the sklearn library.
  2. Then, we set the range of values for our hyperparameters. The hyperparameters we decided to fine-tune were: The number of decision trees (n_estimators), the depth of each tree (max_depth) and whether the decision trees sampled with replacement (bootstrap).
  3. Using these hyperparameter values, we created rf_random, a collection of ML models. Each of these models contained a specific set of hyperparameter values.
  4. We then trained these models using our x_train and y_train data sets.
  5. Using these data sets, we then had a best set of hyperparameter values. We printed them out using the line “rf_random.best_params_“.
  6. With these best values, we created our best ML model, rf_max. We then trained it using the x_train and y_train data sets. With our model completed, what remained was to evaluate it

Retrieving the performance results for our RandomForest model, rf_max.

Results for our RandomForest model

 

Looking at the bottom at the bottom of the above picture, we could see that our model fared well in precision, recall and accuracy (There’s a difference between them!) The picture below explains the difference. The results (we were especially interested in recall) were good but we wanted to know if the other models could have done better…

Defintion of accuracy, recall and precision source:www.youtube.com

 

Fine-tuning logistic regression model:

Code for fine-tuning logistic regression parameters

The process of fine-tuning our LR model was largely similar to RF except for a few differences.

  1. There was an additional step of normalizing our features. This to to ensure that all the features, regardless of their absolute scale, have the same weightage assigned to them. The analogy is as such: Imagine you want to compare the cost of living between city A and city B. So you measured and summed up two values: 1. The difference between the price of a car in city A and city B and 2. The difference between the price of a burger in city A and city B.  Since the price difference of the cars is much greater than that of the burgers in absolute terms, in this approach the price difference of the burgers has negligible impact on the final results even though it is important. Likewise, this problem presents itself in logistic regression since some features contain much larger differences than others. Hence, by scaling all values to a range of 0 to 1, normalizing resolves this problem as all the values would then be on the same scale.
  2. The parameters that we tuned were “C”, which is the inverse of regularization strength, and maximum number of iterations.

And the results are shown below…

Results for logistic regression model

The results were good but subpar to the RandomForest model.

Lastly, we moved on to the Support Vector Machine Model (SVM).

Fine-tuning the Support Vector Machine model

Again the steps were largely similar except that the hyperparameters we tweaked were: kernel, C and gamma.

Results for SVM model

The results were good but again it was not good as RandomForest.

Hence, we decided that RandomForest would be our chosen algorithm for Prototype 3. With the algorithm chosen, what was left was to save  our trained RandomForest model such that we could make subsequent predictions easily. Saving ML algorithm with pickle

In the above code, we ‘dumped’ our model, rf_max into a txt file.

Loading our ML model using pickle

On the Rasp Pi, we then retrieved the text file before loading it using pickle. With that, GS was a fully trained RandomForest model, capable of predicting whether there would be an imminent eruption on the basis of past data.  This approach was brilliant because by having our computer do the tough work of training our ML model before passing the finished product to Rasp Pi , it drastically reduced the computational demands on the latter. This allows our Rasp Pi to remain fast and low-power consuming.

This concludes the making of Prototype 3. Adieus!

Part 12. 25th July ’19: Meeting with Prof Benoit and Mr Doray and the birth of our Final Prototype

At the initial stage of this M&T project, we were introduced to Prof Benoit by Dr. Ho. Prof Benoit asked us if we would be interested in creating a volcanic eruption sensor and if we were, he could use his expertise to guide us along. We took up his offer immediately. We were extremely fortunate to have Mr Arnold Doray, Chief Executive Officer of Terra Weather Pte Ltd, to guide us along. After countless hours of hard work, SensorX was born.

So we were glad that a few weeks ago, Prof Benoit and Mr. Doray both indicated their interest to view a live demonstration of SensorX. The live demo had just taken place this afternoon and this post would note down the feedback they gave.

We are glad to say that SensorX was well-received by our mentors. Moving forward, these were some pointers they gave:

  1. Remove our tilt sensor (an aluminum ring suspended in a aluminium cone) because of several reasons.
    1. Our tilt sensor was not precise enough. Yes, it could detect tilt but it could not sense the difference between a 15° and 20° tilt.
    2. Over time, the electrical conductivity of our aluminium cone would deteriorate due to formation of an oxide layer, and its sensitivity would be reduced.
  2. Remove the need for SensorX to be placed on a horizontal surface before starting it. Currently in the market, volcanic seismometers have to be placed level for them to detect seismic signals accurately. If we could eliminate that hassle, then SensorX would have an added advantage.
  3.  Bi-directional sensing. In our setup, the periodic motion of SensorX was only along the X-axis. On actual volcanoes, the movement tends to be multi-directional. Hence, if we moved SensorX along the Y-axis as well, it would be a much better sensor.

After brainstorming, we came up with solutions for the above pointers.

  1. To remove the need for a tilt sensor, we would utilize the value of z-acceleration which is measured by Rasp Pi. From basic physics, we know that when an object is resting on an inclined surface, the normal contact force is smaller as compared to when it is resting on a level surface. Likewise, when SensorX is resting on an inclined surface (i.e. there is tilt), it would record a smaller g force. Hence, we plan to use this value of z-acceleration to determine if there is tilt.
  2. Instead of using the absolute value of z-acceleration to determine tilt, we use the difference between the instantaneous and initial acceleration. In that case, SensorX would work regardless of initial tilt angle.
  3. At this juncture, we agreed that it is better to focus solely on uni-directional displacement. Once we perfect SensorX using uni-directional movement, we will then move on to bi-directional.

Finding initial value of tilt

Checking for tilt

In the above lines of code, we first measure the initial z-acceleration. Once the sensing starts, we would then periodically measure the instantaneous z-acceleration and take the difference between these 2 values of acceleration. If the difference exceeds 0.01 (which corresponds to roughly 20°), SensorX’s tilt condition will be triggered.

Part 13. 26th July ’19:  Improving our Final Prototype

As mentioned in the previous post, thanks to the invaluable feedback given by our mentors, we made significant improvements to SensorX. However, we did not stop there. We made some further changes which we will document in this post.

We removed the features dominant frequencies of Z and Y accelerations from our Machine Learning data sets.

Why?

We removed the Z dominant frequency because it was difficult to have SensorX simultaneously sense the Z-acceleration value for 2 purposes: to detect tilt and to detect vibrations via our ML model. Hence, we decided that we would only use the values of Z-acceleration for tilt detection

We removed the Y dominant frequency because we realized that it was distorting our ML model. We discovered this after first running our ML model based on the data of both X and Y dominant frequencies and running it based on the data of only X dominant frequencies. We realized that in the latter case, the ML model scored higher on accuracy, precision and recall.

At the bottom of the above picture, we can see that the new model achieved excellent results (:

Hence, we decided that for our final prototype, SensorX would only use the dominant frequency of X-acceleration for the vibrations condition.

Part 14. 29th July ’19:  Demonstration of Final Prototype.

Upon completion of our final prototype today, we decided to film SensorX in action. Here is the video.

From our video, you can see that when the vibrations condition was met (the frequency of vibrations was within the specific range), the output on the window was “1” and the line “checking for tilt” was printed. Then, when we tilted the setup beyond the threshold angle, the tilt condition was met. SensorX was triggered: SenseHat’s LED displayed “Triggered” and a warning message was sent via Bluetooth to the phone.

And that ladies and gentlemen, is our final prototype of SensorX. We hope you like it as much as we do. Now, to the last obstacle: FINAL PRESENTATION.