Fourier Analysis of Acceleration Data from Accelerometer

Lets first read in the signal from a excel sheet in the working directory. This file includes 1000 samples at intervals of 30 milliseconds.

library(xlsx)
## Loading required package: rJava
## Loading required package: xlsxjars
walkingData <- read.xlsx("./data/SUB1_test_SUB2_training.xlsx",sheetIndex=1,header=TRUE)
head(walkingData)
##   subject1_test subject2_training
## 1        0.3142           -0.2067
## 2        0.2900           -0.0161
## 3        0.1477            0.2470
## 4        0.0671           -0.1074
## 5        0.0375           -0.3706
## 6        0.0295           -0.0939

Time signal exploration

let’s see their scatter plot

# plotting just the first 100 values in both time signals
plot(walkingData$subject1_test[1:100], walkingData$subject2_training[1:100], main = "Scatter plot", xlab="Subject 1 Test Acceleration in Y direction (RED)", ylab="Subject 2 Training Acceleration in Y direction (BLUE)", col = c("red", "blue"))

## Let just pass a regression line through to 
model <- lm(subject1_test ~ subject2_training, walkingData)
abline(model, lwd = 1, col = "green")

Looking at the spread of the scatter plot and the slope of the regression line it does not seem that Subject 1 test data can be used as a good predictor of Subject 2 training data.

Analysis of the Signals in the Frequency Domain

Now let’s look at the frequency spectrum of the time signals using Fast Fourier Transform (FFT).

Subject 1 Test Data FFT

# load the library for FFT related functions
library(stats)

# of samples
N <- 1000

# perform fft and take absolute values ignoring the imaginary part of complex numbers
fft_sub1_coef = abs(fft(walkingData$subject1_test))

# normalize the coeficients
# We scale by N and multiple by 4 to get total energy  (see diagnostics)
norm_fft_sub1_coef = 2* fft_sub1_coef/(N/2)

# get the first 500 coefficients, single sided spectrum
mainCoeffs_sub1 = norm_fft_sub1_coef[1:(N/2)]

# plot the coefficients 
plot(mainCoeffs_sub1, type = "h", main = "Frequency Spectrum for Subject 1 test data", xlab="Frequency Index", ylab="Magnitude", col = "red")

It is easy to see that the signal has one fundamental frequency and harmonics at 2f, 3f, 4f and so on..

Let’s zoom in to the first 100 points where the fundamental frequency component with the highest value occurs

# plot the first 10th of the spectrum only for better visibility
plot(mainCoeffs_sub1[1:(N/10)], type = "h", main = "1/10th Frequency Spectrum for Subject 1", xlab="Frequency Index", ylab="Magnitude", col = "red")

Extract highest value

# Find the maximum acceleration amplitude of the frequency component
max_sub1 <- max(mainCoeffs_sub1)
max_sub1
## [1] 0.301805

Extract the frequency index where the highest value occurs

# Find the bininded at which the maximum acceleration is found
maxindex_sub1 <- which.max(mainCoeffs_sub1)
maxindex_sub1
## [1] 59

Let’s look at the Harmonic frequencies only

magnitude_Subject_1 <- c(mainCoeffs_sub1[maxindex_sub1], mainCoeffs_sub1[maxindex_sub1*2], mainCoeffs_sub1[maxindex_sub1*3], mainCoeffs_sub1[maxindex_sub1*4], mainCoeffs_sub1[maxindex_sub1*5])

Harmonics_Frequency_Index_Subject_1 <- c(maxindex_sub1, maxindex_sub1*2, maxindex_sub1*3, maxindex_sub1*4, maxindex_sub1*5)

plot(Harmonics_Frequency_Index_Subject_1, magnitude_Subject_1, type = "h", main="Harmonic Frequencies of Subject 1", xlab="Frequency Index", ylab="Magnitude", col = "red")

Let’s look at the energy concertrated around the fundamental frequency +/- 8

sum(mainCoeffs_sub1[(maxindex_sub1-8):(maxindex_sub1+8)])
## [1] 1.115177

The above value seems to cover the range of amplitudes in the time domain

Now let’s interpret the fundamental frequency index value in terms of actual frequency. Each bin in the FFT result corresponds to a frequency of (binindex * fs)/N where binindex is the index of the bin, fs is the frequency with which the input signal is sampled and N is the number of points. In this case the sampling frequency is 1/0.03 which is 33.33 Hz. N is 1000.

# FFT Bin Spacing

# So in our case (binindex * fs)/N  is (binindex * 33.33/1000)

# Frequency bin indexes
F <- (0:N)*33.33/1000

# frequency in Hz at the index value
F[maxindex_sub1]
## [1] 1.93314

This value in Hz seems reasonable for a person to have a periodic walking signal at roughly 2 cycles per second

Subject 2 Training Data FFT

# of samples
N <- 1000

# perform fft and take absolute values ignoring the imaginary part of complex numbers
fft_sub2_coef = abs(fft(walkingData$subject2_training))

# normalize the coeficients 
# We scale by N and multiple by 2 to account for total energy in the negative and positive peaks of the time signal 
norm_fft_sub2_coef = 2* fft_sub2_coef/(N/2)

# get the first 500 coefficients, single sided spectrum
mainCoeffs_sub2 = norm_fft_sub2_coef[1:(N/2)]

# plot the coefficients 
plot(mainCoeffs_sub2, type = "h", main = "Frequency Spectrum for Subject 2", xlab="Frequency Index", ylab="Magnitude", col = "blue")

It is easy to see that the signal also has one fundamental frequency and harmonics at 2f, 3f, 4f and so on..

Let’s zoom in to the first 100 points where the fundamental frequency component with the highest value occurs

# plot the first 10th of the spectrum only for better visibility
plot(mainCoeffs_sub2[1:(N/10)], type = "h", main = "1/10th Frequency Spectrum for Subject 2", xlab="Frequency Index", ylab="Magnitude", col = "blue")

Extract highest value

# Find the maximum acceleration amplitude of the frequency component
max_sub2 <- max(mainCoeffs_sub2)
max_sub2
## [1] 0.1250704

Extract the frequency index where the highest value occurs

# Find the bininded at which the maximum acceleration is found
maxindex_sub2 <- which.max(mainCoeffs_sub2)
maxindex_sub2
## [1] 60

Interestingly, this is the same value we got for Subject 1 !!!

Let’s look at the Harmonic frequencies only and their magnitude

magnitude_Subject_2 <- c(mainCoeffs_sub2[maxindex_sub2], mainCoeffs_sub2[maxindex_sub2*2], mainCoeffs_sub2[maxindex_sub2*3], mainCoeffs_sub2[maxindex_sub2*4], mainCoeffs_sub2[maxindex_sub2*5])

Harmonics_Frequency_Index_Subject_2 <- c(maxindex_sub2, maxindex_sub2*2, maxindex_sub2*3, maxindex_sub2*4, maxindex_sub2*5)

plot(Harmonics_Frequency_Index_Subject_2, magnitude_Subject_2, type = "h", main="Harmonic Frequencies of Subject 2", xlab="Frequency Index", ylab="Magnitude", col = "blue")

Let’s look at the energy concertrated around the fundamental frequency +/- 8

sum(mainCoeffs_sub2[(maxindex_sub2-8):(maxindex_sub2+8)])
## [1] 0.4301378

The above value seems to cover the range of amplitudes in the time domain

Now let’s interpret at the fundamental frequency in terms of actual frequency.

# frequency in Hz at the index value
F[maxindex_sub2]
## [1] 1.96647

Summary of results

Subject 1 Test Data

Fundamental Frequency Binindex: 59 equivalent to 1.93314 Hz with magnitude 0.301805 \(m/s^{2}\) that eventually decays with odd and even harmoics

Subject 2 Training Data

Fundamental Frequency Binindex: 60 equivalent to 1.96647 Hz with magnitude 0.1250704 \(m/s^{2}\) eventually decays with odd and even harmoics

This is an interesting finding. While the fundamental frequency is the same, the amplitude is considerably different for the fundamental frequency. This finding suggests that we may more or less walk at the same fundamental frequency, but the amplitude of this frequency will wary considerably.

Last updated on Monday Feb 4nd, 10:00 AM, 2015.

References:

  1. “Documentation.” Fast Fourier Transform. N.p., n.d. Web. 01 Feb. 2015. http://www.mathworks.com/help/matlab/ref/fft.html.

  2. “The Scientist and Engineer’s Guide ToDigital Signal Processing By Steven W. Smith, Ph.D.” Harmonics. N.p., n.d. Web. 02 Feb. 2015. http://www.dspguide.com/ch11/5.htm.

  3. Binindex explanation on stack overflow: http://stackoverflow.com/questions/4364823/how-do-i-obtain-the-frequencies-of-each-value-in-a-fft/4371627#4371627

  4. Neto, João. “Fourier Transform: A R Tutorial.” N.p., n.d. Web. http://www.di.fc.ul.pt/~jpn/r/fourier/fourier.html.