Experimental Study of Simple Harmonic Motion

In the classroom, we often use a stopwatch to conduct the pendulum experiment, or simple harmonic movement experiment. Here is a challenge, can we produce a real graph of its movement and see what is the instantaneous angular position and velocity, that's a lot more information and fun.

First question, we need to decide the pendulum body is a weightless cord or a rigid uniform rod. The cord approach seems to be easier. From the practice of building one, I have the following tradeoff considerations: The easiest way to hang a pendulum system might be hang it to the upper edge of your door. That gives your ~2m pendulum length without doing any structural building work. But it does need the swing won't touch the door surface, which simply ruins the whole experiment. So the plane it swing should be accurately parallel to your wall/door surface. A weightless cord tends to be thin, it can spin easily, and complicate the measurement of swing angle. We want to use one measurement to represent the swing state. Thin cord, such as fish line, can be elastic and stretchable, that affects one of our most important constants measured by us and used in the equation, which is the length of the pendulum. Some may also be affected by the temperature. The weight mass hanging at the end of the cord has to be heavy enough so that the weight of the cord becomes negligible. Please comment if you agree or disagree with them, or you have other design tradeoff ideas. To study this problem, we need a device which is so light that its weight can be ignored and we still treat the pendulum system as a rigid uniform rod. I am using a COTS wearable electronic controller, which delivers the gyro, accelerometer, and angle information to us via bluetooth connection. These measurements will be stored into a mobile phone app data file. After that, we will analyze the data for our simple harmonic motion experiment. The numerical analysis focuses on the following topics: 1) Predict the pendulum oscillation period 2) Programmably collect the pendulum simple harmonic movement experiment data 3) Use kmean to group data and remove outliers in the analysis process 4) Use short-time FFT to estimate the pendulum oscillation frequency

Bluetooth measurement apparatus http://ebay.us/5BLNwm?cmpnId=5338273189

Android phone App: Go to Google playstore, search M2ROBOTS and install the control App. In case it is difficult to access Google playstore, visit my personal homepage for alternative app download method https://xiapeiqing.github.io/doc/intro_androidapp...

wood rod

few 3D printed parts

saw blades or similar metal material

There are many articles and books introducing the pendulum equation derivation, including your curriculum physics book. Such content might be better not being repeated here again. Only the final conclusion is listed here regarding the topic of "simple harmonic motion".
In order to know the period of a pendulum, all we need to know is the length of the pendulum, denoted as "l", in meters.

If we are reasonably sure the weight is located almost completely at the end of a weightless cord hanging in a pivot, and the pendulum is swinging at small angles θ, say less than 15°, the period T1 of such pendulum is given by:

T1 = 2*pi*(l/g)^0.5

g = gravity acceleration, approximately 9.8 m/s^2

If the weightless cord is replaced by a rigid uniform rod, again of length l, its simple harmonic motion period T2 is given by T1 = 2*pi*(2l/3g)^0.5

Effectively it has the same period as a weightless cord pendulum being two thirds of the rigid uniform rod length.

This is the background, and we can start preparing our experiment.

To build the pendulum structure, we 3D print some parts, and recycle something we already have.
The overall pendulum structure is shown in Fig.1. It's a mixture of 3D print parts together with some hand crafted parts and a long piece of wood rod from Lowe‘s.

The 3D printed part in Fig.2 is hang on the top edge of a door, because our door is an easy flat surface for us to hang something. STL file download link:

https://xiapeiqing.github.io/doc/kits/pendulum/pen...

The green part in Fig.3 connects the wood rod to a blade, and the blade sits on top of two pieces of rail mounted on the earlier 3D printed door hanger. STL file download link: https://xiapeiqing.github.io/doc/kits/pendulum/ro...

The two pieces of rail is made by breaking an old saw blade in half, see Fig. 4. The part in Fig.2 has prepared the right slot size for them. Ideally we can make a "V" shaped notch in those two saw blades using a file. A reasonably sharp edged metal, such as a single edge razor blade, or any hand made metal piece, can sit inside the "V" shaped notches. The reason we need a smaller contact area is to reduce the kinetic energy lost while swinging.

The last 3D printed part in Fig.5 is a small tray to hold the electronic measurement apparatus.

The download link: https://xiapeiqing.github.io/doc/kits/pendulum/ro...

The bluetooth measurement apparatus generates angle estimate, gyro measurement and accelerometer measurement. All these data is available to us through bluetooth wireless link.

We are going to conduct multiple experiments by deploying this apparatus at different position of the pendulum arm, and see the differences.

There are two doable methods for the experimental data collection before we analysis the acquired dataset:

1) Use the Android phone App specified in the requirements section to log all the measurement produced by the apparatus into a data file stored in your phone's SD card. We can copy the file and post process the information.

2) Use a bluetooth enabled computer, a PC, a laptop or a RaspberryPi mini-computer to establish a bluetooth connection to the apparatus and read the data for either realtime or offline analysis.

There exist both pros and cons for each method, we are going to try both and tell the difference in this instructable.

For method (1) using the android app, once we are in the android App control interface, the telemetry data sent from the bluetooth measurement apparatus to the android phone will be recorded into a datalog file named m2flightDatayyyymmdd_hhmmss.txt. It can be found in your android phone's Download/m2LogFiles folder. The folder "Download" is a pre-existing folder in your phone's android OS and "m2LogFiles" is a folder the App created. The filename content yyyymmdd_hhmmss is the way to encode the experiment starting time (year, month, day, hour, minute and sec) in the file name.

Each line in the log file is one record. It starts with the event timestamp, preamble string "eam:", followed by 4 triplet data, which are:

Accelerometer XYZ axis reading in raw sensor hardware register readback values

Gyroscope XYZ axis reading in raw sensor hardware register readback values

Magnetometer XYZ axis reading in raw sensor hardware register readback values

onboard estimated Roll/Pitch/Raw in degree

The datafile created using computer python program will use identical data file format, so that the program we use in the data analysis step won't be bothered with the data source being produced by our python program or android app.

Let's start coding using method (2).

To interact with the bluetooth measurement apparatus, two flavors of SDK are provided:

1) Python SDK, which can be installed by "pip3 install m2controller", python3 is the language used. The examples of user application code is stored in https://github.com/xiapeiqing/m2robots/tree/maste... For this experiment, we will use the python script pendulum1.py

2) Java SDK, which is not used in this instructable because we want later visualization and analysis of the acquired pendulum data, which might take a little bit more effort for us to program in Java.

The python3 data collection program source code contains many comments for the code functionality details. A snapshot of the source code is provided here.

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
from m2controller import m2controller
from m2controller import m2Const
import signal
import time
import datetime
import usrCfg
import pendulum2
requestExit = False
################################################################
# we want to use the same log file naming convention so that the data analysis module, pendulum2.py, can be agnostic to how we get the log data file 
################################################################
logfilename = "m2flightData%s.txt"%(datetime.datetime.fromtimestamp(time.time()).strftime('%Y%m%d_%H%M%S'))
dataLogfile = open(logfilename,"w") 
def signal_handler(sig, frame):
    global requestExit
    print('user Ctrl-C to exit program execution')
    requestExit = True
    
signal.signal(signal.SIGINT, signal_handler)
################################################################
# upon every measurement data becomes available at 20Hz rate, this "callback" function will be summoned 
################################################################
def callbackfunc(telemetry):
    strTimeStamp = datetime.datetime.fromtimestamp(time.time()).strftime('%H:%M:%S.%f')[:-3]
    dataStr = "%s,eam:%d,%d,%d, %d,%d,%d, %d,%d,%d, %2.1f,%2.1f,%2.1f\n"%(strTimeStamp,
        telemetry['m_fAccelHwUnit'][0],telemetry['m_fAccelHwUnit'][1],telemetry['m_fAccelHwUnit'][2],
        telemetry['m_fGyroHwUnit'][0], telemetry['m_fGyroHwUnit'][1], telemetry['m_fGyroHwUnit'][2], 
        telemetry['m_fMagHwUnit'][0],  telemetry['m_fMagHwUnit'][1],  telemetry['m_fMagHwUnit'][2], 
        telemetry['m_fRPYdeg'][0], telemetry['m_fRPYdeg'][1], telemetry['m_fRPYdeg'][2])
    ################################################################
    # we print the data string to the screen and save them into the log file
    ################################################################
    print(dataStr)
    dataLogfile.writelines(dataStr) 
################################################################
# initialize the controller, remember to set the field BleMACaddress to be your device's MAC address
################################################################
# TODO: let's initialize the BleMACaddress if not being set by user.
controller = m2controller.BleCtrller(m2Const.etDebian,callbackfunc,usrCfg.BleMACaddress)
controller.connect()
    
while True:
    ################################################################
    # wait for measurement data created and sent from the pendulum measurement apparatus
    ################################################################
    controller.m_CommsTunnel.waitForNotifications(1.0)
    if requestExit:
        ################################################################
        # house keeping works here when we finish data logging
        ################################################################
        controller.stop()
        dataLogfile.close()
        break
################################################################
# data collection completed, now let's analyze the log data
################################################################
pendulum2.parseDataLogFile(logfilename)

For long term update, please checkout https://github.com/xiapeiqing/m2robots/blob/maste...

Now let's explain its operation method. This python program is written on top of a pip installable package, named m2controller. The lower level package offers callback mechanism, so that every received measurement update will trigger the callback function we wrote, and save the data into a local log file. The format of log file data content is identical to what is produced by android companion app, so that the data log file created by either python program or andriod companion app is exchangeable.

The user ctrl-C signal, captured by the operating system, is passed to the program, and stop the infinite loop waiting for the new arrival of measurement data.

Till now, log file has been successfully created, and this program will call the analysis program to study our experiment results.

Here are two experiments, and the comparison shows the very noticable difference by attaching a 7gram device at different locations.

In Fig.2, we use a scale to determine the actual weight of this bluetooth measurement apparatus.

Fig.3 depicts the pendulum setup where the 7gram device is attached to the lower end of the pendulum. Setup configuration in Fig.4 has the 7gram mass located much closer to the swinging pivot.

Fig.5 is a closeup view of the pendulum structure.

The bluetooth measurement apparatus weighs ~7gram, which weighs much less than a ~1.6meter long wood stick. Use the assumption of "rigid uniform rod", and we have this pendulum period equation, T1 = 2*pi*(2l/3g)^0.5

To get the gravity constant, we can use 9.8m/s^2. But a more accurate gravity constant at any given geolocation can be retrieved from this web service:

https://www.wolframalpha.com/widgets/view.jsp?id=e...

For san francisco, it is 9.81278m/s^2

The pendulum length is measured to be 64.5''

2*pi*sqrt(2*64.5*0.0254/(3*9.81278)) gives the expected pendulum period of 2.0962(sec).

Let's see whether it agrees with our experiments.

In the 1st experiment, the pendulum setup has the 7gram device attached to the lower end of the pendulum. My log file could be downloaded in:

https://xiapeiqing.github.io/doc/kits/pendulum/pen...

Rename it to "PendulumTestData.txt" and put it in the same folder of the python analysis program. A snapshot of the source code is provided here.

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
import csv
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
from datetime import datetime, timedelta
import seaborn as sns
from sklearn.cluster import KMeans
from collections import Counter
################################################################
# this function runs the data file analysis work
################################################################
def parseDataLogFile(datafilename):
    ################################################################
    # extract data in the comma seperated data log file (CSV) and save the content in each column into one float-type variable 
    ################################################################
    with open(datafilename) as csvfile:
        readCSV = csv.reader(csvfile, delimiter=',')
        timestampS = []
        fAccelHwUnit_x = []
        fAccelHwUnit_y = []
        fAccelHwUnit_z = []
        fGyroHwUnit_x = []
        fGyroHwUnit_y = []
        fGyroHwUnit_z = []
        fMagHwUnit_x = []
        fMagHwUnit_y = []
        fMagHwUnit_z = []
        fRPYdeg_r = []
        fRPYdeg_p = []
        fRPYdeg_y = []
        for row in readCSV:
            try:
                x = datetime.strptime(row[0].split(',')[0],'%H:%M:%S.%f')
                timestampS.append(timedelta(hours=x.hour,minutes=x.minute,seconds=x.second,microseconds=x.microsecond).total_seconds())
                fAccelHwUnit_x.append(float(row[1][4:]))
                fAccelHwUnit_y.append(float(row[2]))
                fAccelHwUnit_z.append(float(row[3]))
                fGyroHwUnit_x.append(float(row[4]))
                fGyroHwUnit_y.append(float(row[5]))
                fGyroHwUnit_z.append(float(row[6]))
                fMagHwUnit_x.append(float(row[7]))
                fMagHwUnit_y.append(float(row[8]))
                fMagHwUnit_z.append(float(row[9]))
                fRPYdeg_r.append(float(row[10]))
                fRPYdeg_p.append(float(row[11]))
                fRPYdeg_y.append(float(row[12]))
            except:
                pass
        timestampS = np.asarray(timestampS)
        timestampS = timestampS - timestampS[0]
        fAccelHwUnit_x = np.asarray(fAccelHwUnit_x)
        fAccelHwUnit_y = np.asarray(fAccelHwUnit_y)
        fAccelHwUnit_z = np.asarray(fAccelHwUnit_z)
        fGyroHwUnit_x = np.asarray(fGyroHwUnit_x)
        fGyroHwUnit_y = np.asarray(fGyroHwUnit_y)
        fGyroHwUnit_z = np.asarray(fGyroHwUnit_z)
        fMagHwUnit_x = np.asarray(fMagHwUnit_x)
        fMagHwUnit_y = np.asarray(fMagHwUnit_y)
        fMagHwUnit_z = np.asarray(fMagHwUnit_z)
        fRPYdeg_r = np.asarray(fRPYdeg_r)
        fRPYdeg_p = np.asarray(fRPYdeg_p)
        fRPYdeg_p = fRPYdeg_p - np.mean(fRPYdeg_p)
        fRPYdeg_y = np.asarray(fRPYdeg_y)  
    ################################################################
    # we need accurate estimate of the sampling frequency for precise oscillation period estimate
    ################################################################
    FsHz = getSamplingIntervalS(timestampS)
    ################################################################
    # use pitch component in the attitude heading reference system output for pendulum period analysis
    ################################################################
    analyze_timeSequence(timestampS,fRPYdeg_p,FsHz,'pitch')
    ################################################################
    # use acceleromter raw measurement output for pendulum period analysis
    ################################################################
    analyze_timeSequence(timestampS,fAccelHwUnit_x,FsHz,'accel')
    ################################################################
    # use gyro raw measurement output for pendulum period analysis
    ################################################################
    analyze_timeSequence(timestampS,fGyroHwUnit_y,FsHz,'gyro')
    print('done, congratulations :-)')
    plt.show()
   
################################################################
# in bluetooth communication process, there is a rare chance that the data comm packet could be lost
# we use K-mean to isolate the 20Hz meaurement data from outliers, which are caused by dropped packet 
# dive into "signal and system for more details"
################################################################
def getSamplingIntervalS(timestampS):
    plt.figure()
    sampleIntervalS = np.diff(timestampS)
    sns.distplot(sampleIntervalS)
    plt.ylabel('histogram')
    plt.xlabel('measurement interval(s)')
    clusterCnt = 5
    km = KMeans(n_clusters = clusterCnt)
    km.fit(sampleIntervalS.reshape(-1,1))
    centroids = km.cluster_centers_
    elemCnt = Counter(km.labels_)
    occurrenceCnt = []
    for ii in range(clusterCnt):
        occurrenceCnt.append(elemCnt[ii])
    FsHz = 1/centroids[occurrenceCnt.index(max(occurrenceCnt))]
    return FsHz
################################################################
# use spectrometer, i.e., short time FFT to get the frequency component, peak bin is our best estimate of pendulum oscillation
################################################################
def analyze_timeSequence(timestampS,timeSeqData,FsHz,strComment):
    fig, (ax1, ax2) = plt.subplots(nrows=2)
    ax1.plot(timestampS, timeSeqData, marker='o', markerfacecolor='blue', markersize=2, color='skyblue', linewidth=1)
    ax1.set_title("pendulum time domain measurement -- %s"%strComment)
    ax1.set_xlabel("sampling time(second)")
    ax1.set_ylabel(strComment);
    NFFT = 2048  # the length of the windowing segments
    Pxx, freqs, bins, im = ax2.specgram(timeSeqData, NFFT=NFFT, Fs=FsHz, noverlap=NFFT/2)
    ax2.set_title("Spectrogram")
    ax2.set_xlabel("samples")
    ax2.set_ylabel("frequency(Hz)");
    # The `specgram` method returns 4 objects. They are:
    # - Pxx: the periodogram
    # - freqs: the frequency vector
    # - bins: the centers of the time bins
    # - im: the matplotlib.image.AxesImage instance representing the data in the plot
    pkresult = np.where(Pxx == np.amax(Pxx))
    oscFreqHz = freqs[pkresult[0][0]]
    print('pendulum oscillation Freq(Hz)=%f, Period(Sec)=%f, estimation data source: %s'%(oscFreqHz,1/oscFreqHz,strComment))
    return 1/oscFreqHz
################################################################
# should we run this program independently, i.e., not being called by pendulum1.py,
# we define a default log data file name to be analyzed
################################################################
if __name__ == "__main__":
    defaultFilename = './PendulumTestData.txt'
    import os.path
    if os.path.isfile(defaultFilename):
        parseDataLogFile(defaultFilename)
    else:
        print ("default log file %s not existing"%defaultFilename)

For long term update, please checkout https://github.com/xiapeiqing/m2robots/blob/maste...

The source code contains detailed comments, let's give a high level summary of the mathematical estimation here.

1) We first read the CSV file content into the computer, using a python package called "csv". We have periodical measurement.

21:34:26.362,eam:0,-128,14464, -8,144,-96, 2112,-1280,1664, -0.5,-5.5,40.5

21:34:26.373,eam:128,0,14272, -8,136,40, 2112,-1280,1664, -0.5,-6.5,40.0

21:34:26.412,eam:448,-64,14208, -8,136,24, 2176,-1280,1664, -0.5,-7.5,40.5

21:34:26.462,eam:448,-128,14272, -8,120,16, 2176,-1280,1664, -0.5,-8.0,40.5

2) Since the measurement rate is so critical and directly introduces pendulum period estimation error, we want to estimate them. Our nominal measurement interval is 50ms, i.e., 20Hz. Average over all measurements seems OK, but we occasionally lose data transmission packet, the update interval becomes 100ms or 150ms, ...

If we plot the occurrence of these data, see Fig.1, as a human, we can easily have an eyeballing value of 0.05sec. However, can we do better than that?

We need to use classification method to only select the good ones for averaging computation. Python has toolbox named KMeans to help us with clustering, or say classification. These concepts are used in many big data and AI areas.

3) Fig.2 contains two images. The top plot is a time-domain sequence of our swinging angle measurement in deg. By referencing to the x-axis timestamp in Second, we can read approximately 22.5 cycles in 50second, which translates to 2.22 Sec pendulum period. Is there a way to automate this process and have more accurate estimate? Yes, we can use mathematical tool called spectrogram, which uses a small chunk of measurement data and tell us its frequency, see the figure below. The y-axis reading for the darkest line is the pendulum oscillation frequency. Being a horizontal line confirms the pendulum oscillation didn't change at all throughout the experiment. The inverse value of oscillation frequency is the pendulum oscillation period.

The final report made by the program is a text summary:

pendulum oscillation Freq(Hz)=0.449224, Period(Sec)=2.226059, estimation data source: pitch

We can find our earlier eyeballing hand calculation result, 2.22sec, is fairly consistent with the program computed value.

Compared to 2.0962(sec) theoretically computed value, we have ~5% remaining error. How to get rid of them? Remember the assumption is "rigid uniform rod"? Even a 7 gram extra weight seems trivial, it is the biggest cause of the remaining error.

We now move the device, close to the pivot. See previous step for a close-up photo. The log file I created can be downloaded here:

https://xiapeiqing.github.io/doc/kits/pendulum/pen...

Run the same analysis steps, and we get Period of 2.089867(Sec), see Fig.3, which is almost identical to the theoretical prediction. Great!

Since we have not only swinging angle measurement, but also gyroscopic measurement and accelerometer measurement at the same rate. Run the same analysis for the other two measurement, we get results in Fig.4 and 5. Estimates from all three measurement sources agree, which make us more confident on the success of our experiment.

Here is the result as the final output of python program running:

pendulum oscillation Freq(Hz)=0.478499, Period(Sec)=2.089867, estimation data source: pitch

pendulum oscillation Freq(Hz)=0.478499, Period(Sec)=2.089867, estimation data source: accel

pendulum oscillation Freq(Hz)=0.478499, Period(Sec)=2.089867, estimation data source: gyro

Last thought in this step, how can the estimation results be exactly identical using different input data source? This is counter-intuition. I will leave this question to the readers. Here is a hint: remember we are using the short-time FFT to estimate the oscillation frequency? In digital domain, the frequency estimate is given in discrete frequency bins instead of a floating number estimate.

There are few categories of future work recommendations.

In earlier step, we manage to reduce our experiment error from ~5% to less than 1%, can we do better than that? Noticing the oscillation magnitude decreases exponentially, one contributing factor can be the air drag caused when swinging the pendulum. The cross section of the pendulum might need to be modified to be of streamline shape so as to reduce the aerodynamic drag.

Can we apply a time-varying gain learned using adaptive filter techniques to output a constant peak-magnitude signal. In the meantime, correlate the magnitude attenuation will external forces.

We can hardly find anything simpler than the "simple harmonic movement". Can we use the facilities we analyze the pendulum to analyze something more complicated, a sport activity, a water rocket launch sequence, etc?

Happy hacking