hello guys,
I’ve been running the script below in R to plot the frequency vs RMS acceleration diagrams. How can I combine the dataX, Y & Z curves into one diagram, and plot each line without colour into various form of dotted line since my publication can only be in black?
So far I can only ran the script to analyze the X, Y, and Z separately. Modify the line "data<-dataZ" as needed for each axis.
#Gulf Coast Data Concepts
#June 22, 2016
#file name: "car_vibration.r"
#R script for processing X2 accelerometer data into
#a frequency spectrum plot
#----------------------------------------------------
#This script takes a segment of data (a block) and converts
#the time series into a frequency series (Fast Fourier Transform)
#The FFT output represents the energy level for each frequency bin
#The analysis repeats by pulling successive blocks of data
#throughout the data set. Each block is plotted during the analysis.
#The results of each FFT are summed and averaged. The output data is
#saved to a text file
#-----------------------------------------------------
######################################################################
#Functions must be declared first before the script can utilize them #
######################################################################
#fileList is a simple function to build an array of text strings containing
#the files to analyze. This particular script will only use one data file
#as outlined in the Helmholtz test procedure. The fileList is kept here to
#maintain continuity with other scripts. The list will contain only one file
######################################################################
#these lines will build a list using the path, directory, and file numbers defined
#at the beginning of the script
fileList<-function(start, finish, dr, dir){
dataFiles<-array()
#prefix name of the data files
filename<-"DATA-"
g<-as.numeric(start)
n<-1
while(g<(as.numeric(finish)+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(dr,":", "\\", dir, "\\", name, sep="")
dataFiles[n]<-path
g<-g+1
n<-n+1
}
return(dataFiles)
}
#the next few lines of code build a Hann Window filter
#for processing successive FFTs. It's good to filter the
#input data such that there are no discontinuities in the
#transitions between blocks of data. Hann Window is a common
#way to do this.
hannWindow<-function(block){
#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
}
return(hann.window)
}
#################################################################
#Function definitions are complete. Before starting the main #
#script, we need to collect information from the user regarding #
#files and analysis parameters. #
#################################################################
#First, print a note to the user that this script is intended
#for Windows PCs and the file directory path is not compatible
#with Mac or Linux.
print("##############################################")
print("Note, this script is intended for Windows based PCs.")
print("The file path locations are not compatible with Mac/Linux")
print("##############################################")
#ask for the drive location of the data files
#this is the drive or path that contains the GCDC directory
drive<-(-1)
while(drive<0){
drive<-readline("Enter drive location letter: ")
drive<-ifelse(grepl("[a-z,A-Z]", drive), as.character(drive), -1)
}
#the directory containing the data files
#ask for the drive location of the data files
#this is the drive or path that contains the GCDC directory
directory<-(-1)
while(directory<0){
directory<-readline("Enter drive location letter(GCDC): ")
if(directory==""){
directory="GCDC"
break
}
directory<-as.character(directory)
}
#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
#ask for the sequential files to analyze
fileNumStart<-(-1)
while(fileNumStart<0){
fileNumStart<-readline("Enter the file to analyze: ")
fileNumStart<-ifelse(grepl("[0-9]", fileNumStart), as.numeric(fileNumStart),-1)
}
#the script uses only one file so the fileList will include only
#one entry
fileNumFinish<-fileNumStart
#build a vector list of files to analyze
#the first file is used to determine the start_time
dataFiles<-fileList(fileNumStart, fileNumFinish, drive, directory)
#open the first file and read the first five lines
#the third line contains the start_time
lines<-readLines(dataFiles[1],6)
#parse the lines
date<-unlist(strsplit(lines[3],","))
#combine the 2 and 3 entries into a date/time
d<-paste(date[2],date[3])
#convert the date/time text into a R type date
startDate<-as.POSIXlt(d)
#read the sample rate setting from the first data file
#the X2 logger records sample rate on line 6 of the file header
temp<-unlist(strsplit(lines[6],","))
SR<-as.numeric(temp[2])
#is user ready to start analysis
test<-(-1)
while(test<0){
test<-readline("Press Enter to start analysis! ")
if(test==""){
break
}else{
test<-(-1)
}
}
############################################################
#the following parameters affect the analysis algorithm #
############################################################
#block is a segment of data processed through the FFT algorithm
#The script will process 4 seconds of data thru the FFT
block<-4*SR
#the FFT function can't handle blocks of data larger than 4096
if(block<0 | block>4096){
print("Invalid block size")
block<-4096
}
#the script progresses through data set pulling
#blocks of data. Each block of data can overlap the
#previous block of data. Overlapping the data improves
#the signal-to-noise of the results but increases
#computations and time
#define the FFT overlap, 1=0%, 2=50%, 4=75%, etc
overlap<-4 #75% of data will overlap previous block of data
#and 25% will be new data from the file
######################################################################
# Analysis parameters finished. Now, set up variables before #
# starting the main part of script #
######################################################################
#check system clock for start of analysis
startTime<-Sys.time()
#cfactor is used to convert the raw data into "g". The factor
#is different depending on the logger type and sensor range.
#X16 series -> 2048
#X2, low gain -> 6554
#X2, high gain -> 13108
#X200 -> 163.8
cfactor<-6554
#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)
#build a history of FFTs
fftHistory<-array(0)
#build a Hann Window filter
#when processing successive FFTs it's good to filter the
#input data such that there are no discontinuities in the
#transitions. Hann Window is a common way to do this.
hann.window<-hannWindow(block)
#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))
#####################################################################
# START OF MAIN SECTION #
#####################################################################
#this loop will read the input file into three arrays
#there is only one data file so this loop will complete once
f<-1
while(f<(length(dataFiles)+1)){
print(paste("Opening data file: ", dataFiles[f]))
colNames<-c("time", "Ax", "Ay","Az")
dataIn<-read.table(dataFiles[f], sep=",", comment=";", col.names=colNames, fill=TRUE)
#initialize three arrays to store the xyz data
dataX<-array(0)
dataY<-array(0)
dataZ<-array(0)
#use spline to resample the data into evenly spaced samples. The accelerometer
#sensor streams data at the set sample rate but this may vary depending on
#the accuracy of the sensor's clock. For example, configuring to sample at 400Hz
#may actually result in data recorded at 410Hz. We need the resulting data to be
#very accurate and consistent. Therefore, as data arrives to the logger CPU, it is time
#stamped using an independent and accurate real time clock. These time stamps
#are assumed to be accurate. We will resample the data based on
#the time stamps to ensure the time series data is consistent and accurate.
################################################################################
#determine the total time recorded in the file by finding the
#largest time stamp and subtracting the smallest time stamp
timeFile=max(dataIn[,1])-min(dataIn[,1])
#total number of samples that should be in the file assuming
#the sample rate used. In reality, the sensor may produce data
#slower or faster than the set rate, that's why we are resampling
numberSamples=timeFile*SR
#The time stamps are based on the real time clock so these values are
#accurate. However, the arrival of the data could be different than the
#set sample rate due to the sensor's clock error. Therefore, we will
#use a cubic spline algorithm to resample each axis such that the
#sample timing is consistent. This will improve the accuracy of the FFT
#by correcting the sample rate error from the sensor.
resample<-spline(dataIn[,1],dataIn[,2],n=numberSamples,xmin=min(dataIn[,1]),xmax=max(dataIn[,1]))
dataTime<-unlist(resample[1])
dataX<-unlist(resample[2])
resample<-spline(dataIn[,1],dataIn[,3],n=numberSamples,xmin=min(dataIn[,1]),xmax=max(dataIn[,1]))
dataTime<-unlist(resample[1])
dataY<-unlist(resample[2])
resample<-spline(dataIn[,1],dataIn[,4],n=numberSamples,xmin=min(dataIn[,1]),xmax=max(dataIn[,1]))
dataTime<-unlist(resample[1])
dataZ<-unlist(resample[2])
#Convert the raw data into g's
dataX<-(dataX)/cfactor
dataY<-(dataY)/cfactor
dataZ<-(dataZ)/cfactor
#Evaluate the z-axis data which is vertical to the road surface
data<-dataZ
#Another way to look at the data is to use the RMS magnitude all three axis
#data<-sqrt(dataX^2+dataY^2+dataZ^2)
#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 extracts the real portion
#of the FFT. The results of the FFT are summed and then averaged.
i<-1
while((i+block)<length(data)) {
datablock<-data[i:((block-1)+i)] #get the data
data.hann<-hann.window*datablock #apply the filter
data.fft<-fft(data.hann) #perform FFT
#data.fft<-fft(datablock) #no filter
data.mod<-Mod(data.fft) #Mod separates the real and imaginary values
data.mod.half<-data.mod[1:((length(data.mod))/2)] #use the real half of FFT
#correct output so values are in acceleration units (g)
fftOutput<-data.mod.half/(block/2)
#calculate the time stamp for this particular block of data
time<-timeF+(i+block)/SR/60
#combine the FFT and x.axis into one array for plotting purposes
output<-array(c(x.axis, fftOutput), dim=c(length(x.axis),2)) #combine arrays
#plot the output in log scale on the y-axis
#this plot shows the spectral response of the resonating bottle
plot(output, type="l", xlim=c(0,SR/2), ylim=c(1e-6,1e6), log="y", tck=1, xlab="Frequency", ylab="Total RMS Acceleration (g)")
text(round(time, 1), x=70, y=50)
text("elapsed minutes", x=70, y=10)
#Sys.sleep(0.1) #pause so the plots don't flash by too quickly
#add the fft result to the history
fftHistory<-fftHistory+fftOutput
#save the time stamp 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)
#back to the top of loop and repeat everything for
#the next block of data
}
#update timeF to track time into the next data file
timeF<-time
#move to the next file
f<-f+1 #increment f
#back to top of main loop to repeat everything
#for next data file
}
#end of main loop
#convert time into date
timeDate<-startDate+timearray*60
#average out the FFT history
fftHistory<-fftHistory/(length(timearray))
#plot averaged spectral response for the entire data set
output<-array(c(x.axis, fftHistory), dim=c(length(x.axis),2))
plot(output, type="l", tck=1, xlab="Frequency (Hz)", ylab="RMS Acceleration", xlim=c(0,SR/2), ylim=c(0,max(fftHistory[5:length(fftHistory)])), main=dataFiles[1])
text(paste("Data File:",dataFiles[1]), x=SR/4, y=0.005, col="red")
#finished analysis
finishTime<-Sys.time()
#print(paste("total time for analysis: ",(finishTime-startTime)/60," minutes"))
#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_File_", as.character(fileNumStart),".txt", sep="")
columnOutput=c("freq","accel (g)")
write.table(output, outputDataFile, sep="\t",row.names=FALSE, col.names=columnOutput)
