#Gulf Coast Data Concepts
#August 25, 2013
#file name: "train_vibration.r"
#R script for processing X2 accelerometer data
#to find vibration level of a freight train
#define the drive location of the data files
#this is the drive or path that contains the GCDC directory
drive<-"h:"
#the directory containing the data files
directory<-"GCDC"
#prefix name of the data files
filename<-"DATA-"
#the data logger creates files sequentially numbered
#the following parameters define the file numbers to analyze
#for example, files 377 through 380
fileNumStart<-377
fileNumFinish<-380
#to create an animated GIF of the spectral response, first download
#and install the "animation" package into R
#also install graphicsmagick, which is used to combine the png images into GIF
#within R enter the following command:
#saveGIF(source("train_vibration.r"), movie.name="c:\\train.gif", img.name="train", convert = "gm convert", cmd.fun=system, clean=TRUE)
#this next line sets the pause interval between frames
#ani.options(interval=0.05)
#the following parameters affect the analysis algorithm
#block is a segment of data processed through the FFT algorithm
block<-1024
#overlap sets how much of the next block of data
#overlaps the previous block of data
overlap<-4 #defines the FFT overlap, 2=50%, 4=25%, etc
#sample rate of the data set as collected by the X2-1 data logger
SR<-256 #this is the sample rate
#defines the number of trailing FFT results to average when
#comparing to the noise level. Large number increases signal-to-noise
#but introduces a longer time lag
fftTrailSize<-4
#number of initial FFT results to collect at the beginning of
#analysis (first file) to use in creating a noise floor as a
#basis for comparison
fftNoiseSize<-64
#count how many FFT results have been accumulated for noise calc
fftNoiseIncr<-0
#cutoffL and cutoffH define the frequency band to monitor
#everything outside the frequency band is removed
#in the case of a train, we area looking at the infrasonics between 0-20Hz
cutoffL<-1
cutoffH<-20
#define array to hold list of data files
dataFiles<-array(0)
#these lines will build a list using the path, directory, and file numbers defined
#at the beginning of the script
g<-as.numeric(fileNumStart)
n<-1
while(g<(as.numeric(fileNumFinish)+1)){
#add leading zeros if need to make valid file name
if(g<10)
temp<-paste("00",as.character(g),sep="")
else if(g<100)
temp<-paste("0",as.character(g),sep="")
else
temp<-as.character(g)
name<-paste(filename, temp, ".csv", sep="")
path<-paste(drive, "\\", directory, "\\", name, sep="")
dataFiles[n]<-path
g<-g+1
n<-n+1
}
#initialize three arrays to store the xyz data
dataX<-array(0)
dataY<-array(0)
dataZ<-array(0)
#a time value is calculated for each FFT, time increases with
#each new block of data as the script proceeds through a file.
#The time values are stored to timearray. When a new file is opened,
#timeF is used to carry the last time value into the new file data
time<-0
timeF<-time
#timearray stores all the FFT time values and is used later for
#plotting purposes
timearray<-array(0)
#store the total RMS Acceleration for the current FFT
vibRMS<-0
#build a history of total RMS acceleration for each FFT
vibrationHistory<-array(0)
#create filter factor table
#this loop creates a list of numbers used
#to filter the FFT output to just the frequencies
#within the cutoffL/cutoffH band
n<-1
filterFactors<-array(0,dim=c(block/2))
while(n<(block/2+1)){
if((SR*n/block)>cutoffL){
filterFactors[n]<-1
#anything higher than cutoffH is multiplied
#by zero and removed from output
if((SR*n/block)>cutoffH){
filterFactors[n]<-0
}
#lower than cutoffH is removed
}else{
filterFactors[n]<-0
}
n<-n+1
}
#the next few lines of code build a Hann Window filter
#first, initialize the array with a zero
hann.window<-0
#the first element of the hann.window array has a 0 so the
#index will be initialized to '2' to address the next array element
z<-2
#start a loop that will build an element array of factors
#representing a Hann filter
while(z<(block+1)){
#this is the formula for a Hann filter
temp<-0.5*(1-cos((2*pi*z)/(block-1)))
#take the new factor and append it to the hann.window array
hann.window<-append(hann.window, temp)
#increment the index by 1
z<-z+1
}
#initialize a matrix to hold the trailing results of the
#spectral energies
fftTrail<-matrix(0, nrow=block/2, ncol=fftTrailSize)
#this will hold the sum of the trailing data
fftTrailSum<-array(0, dim=c(block/2))
#fftAvg holds the avg spectral energy for the trailing ffts
fftAvg<-array(0, dim=c(block/2))
#array to hold the first FFT results used to
#determine the noise floor
fftNoiseTrail<-matrix(0, nrow=block/2, ncol=fftNoiseSize)
fftNoiseSum<-array(0, dim=c(block/2))
fftNoiseAvg<-array(0, dim=c(block/2))
#anything above the average is considered a noteworthy signal
vibrationSignal<-array(0, dim=c(block/2))
#freqarray saves the dominant frequency (ignoring DC to 0.1 hertz)
freqarray<-array(0)
#build an array of frequency values based on fft size
#this is used for plotting the spectral output
x.axis<-0:(block/2-1)
x.axis<-x.axis*((SR/2)/(block/2))
#this loop will read the input file into three arrays
f<-1 #start at the second line, first line has output file name
while(f<(length(dataFiles)+1)){
message(paste("Opening data file: ", dataFiles[f]))
data<-read.table(dataFiles[f], sep=",", comment=";")
dataX<-c(data[[2]])
dataY<-c(data[[3]])
dataZ<-c(data[[4]])
#X2-1 conversion
#convert raw data into g's
#in low gain, divide raw data by 6554
#in high gain, divide raw data by 13108
dataX<-(dataX)/6554
dataY<-(dataY)/6554
dataZ<-(dataZ)/6554
#calc the rms
data<-sqrt(dataX^2+dataY^2+dataZ^2)
#or use the vertical z-axis only
#data<-dataZ
#this is the main loop to calculate each FFT. First, a block data points
#is pulled from the input data and the Hann filter applied to it. An FFT is
#performed on the filtered data. The "Mod" function calculates the magnitude
#using the real and imaginary output of the FFT. A trailing average is calculated.
i<-1
while((i+block)0.1){
text("TRAIN DETECTED!!!", x=70, y=200)
}
text(round(time, 1), x=70, y=50)
text("elapsed minutes", x=70, y=10)
#calculate the RMS energy
vibRMS<-sqrt(sum(vibrationSignal^2))
#add the new value to the history array
vibrationHistory<-append(vibrationHistory, vibRMS)
#save the time associated with this block of data
timearray<-append(timearray, time)
#increment the index by a percentage of the block
#defined by the overlap variable.
i<-i+(block/overlap)
#pause for 0.1 seconds so the plots don't flash by too quick
#Sys.sleep(0.1)
#uncomment the next line if you are building a animated GIF
#ani.pause()
}
#update timeF to track time into the next data file
timeF<-time
#move to the next file
f<-f+1 #increment f
}
#end of main loop
#combine the vibration history data and time array into one data array
vibrationVStime<-array(c(timearray, vibrationHistory), dim=c(length(timearray),2))
#plot history versus time
plot(vibrationVStime, type="l", tck=1, xlab="Time (minutes)", ylab="Total RMS Acceleration (g)", ylim=c(0,1))
#write the results to the output file
#uncomment the next two lines to have the results stored to a file
#outputDataFile<-paste(drive,"\\","analysis_output_of_Files_", as.character(fileNumStart),"-", #as.character(fileNumFinish),".txt", sep="")
#write.table(vibrationVStime, outputDataFile, sep="\t")