Adding lines to my graph

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)

I can't run this code, because I don't have the data you're using. So, I'm providing an example of plotting multiple variables in same plot using the matplot function.

# simulating data
t <- sort(x = runif(n = 50))
x <- sort(x = runif(n = 50))
y <- sort(x = runif(n = 50))
z <- sort(x = runif(n = 50))

# combining simulated data in a matrix
d <- cbind(x,
           y,
           z)

# plotting multiple line diagram of x, y, z against t
matplot(x = t,
        y = d,
        type = "l",
        lty = 1:3,
        col = 1,
        xlab = "t",
        ylab = "x, y, z")

Created on 2019-02-02 by the reprex package (v0.2.1)

If that's what you're looking for, you'll just have to modify slightly. You may add a legend too.

Hope this helps.

Yes thank you very much, it is what I am looking for. I tried to modify the script by deleting the following, and adding the one you given. However, access is denied.

Here are the data I am using right now.

Besides, is it possible to add a box to the plot?

I'm confused about where am I supposed to put the data files. I renamed the CSV files to DATA-001.csv etc. and placed in C:\GCDC. Is that correct?

I believe your code creates plot separately for the three variables. I'm not sure, though.

If you want to plot the multiple line diagram in Line 344: output<-array(c(x.axis, fftHistory), dim=c(length(x.axis),2)), you can have three outputs corresponding to dataX, dataY and dataZ (changing accordingly in Line 279), and then use:

# outputX is obtained in Line 343 using data<-dataX in Line 279
plot(outputX, 
     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])

# outputY is obtained in Line 343 using data<-dataY in Line 279
plot(outputY, 
     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])

# outputZ is obtained in Line 343 using data<-dataZ in Line 279
plot(outputZ, 
     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])

############################################################

output_column1 <- cbind(outputX[,1],outputY[,1],outputZ[,1])
output_column2 <- cbind(outputX[,2],outputY[,2],outputZ[,2])

matplot(x = output_column1,
        y = output_column2,
        type = "l",
        lty = c(1, 3, 4),
        lwd = c(1, 1, 2),
        col = 1,
        tck=1,
        xlab = "Frequency (Hz)",
        ylab = "RMS Acceleration",
        xlim = c(0, SR/2),
        ylim = c(0, max(fftHistory[5:length(fftHistory)])),
        main = dataFiles[1])
legend("topright",
       legend = c("outputX", "outputY", "outputZ"),
       col = 1,
       lty = c(1, 3, 4),
       lwd = c(1, 1, 2))

This will yield the following:

Does this help?

This is unexpected. Which line creates this error?

What about rect?

Yes, correct!

I am still not there yet, tho I changed the Line 279 with the code you provided, there is an error...

#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
  # outputX is obtained in Line 394 using data<-dataX in Line 279
  plot(outputX, 
       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])
  
  # outputY is obtained in Line 394 using data<-dataY in Line 279
  plot(outputY, 
       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])
  
  # outputZ is obtained in Line 394 using data<-dataZ in Line 279
  plot(outputZ, 
       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])
  
  ############################################################
  
  output_column1 <- cbind(outputX[,1],outputY[,1],outputZ[,1])
  output_column2 <- cbind(outputX[,2],outputY[,2],outputZ[,2])
  
  matplot(x = output_column1,
          y = output_column2,
          type = "l",
          lty = c(1, 3, 4),
          lwd = c(1, 1, 2),
          col = 1,
          tck=1,
          xlab = "Frequency (Hz)",
          ylab = "RMS Acceleration",
          xlim = c(0, SR/2),
          ylim = c(0, max(fftHistory[5:length(fftHistory)])),
          main = dataFiles[1])
  legend("topright",
         legend = c("outputX", "outputY", "outputZ"),
         col = 1,
         lty = c(1, 3, 4),
         lwd = c(1, 1, 2))
  #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)

Please clarify me something.

In which line do you want to plot the multiple line diagram?

Inside the nested while loop in Line 306:

Or, at the end in Line 344:

I thought that you want in Line 344. So, basically what I did was:

  1. copy-pasted from Line 229 to Line 343 thrice (including the one you already had)
  2. changed original Line 279 to data <- dataX, data <- dataY, data <- dataZ
  3. modified original Line 343 to outputX, outputY and outputZ.

Next, I created output_column1 and output_column2, and plotted thereafter as this.

Does that make any sense? :confused:

Thank you very much, it makes sense tho it works differently in my hands :disappointed_relieved:

I want it at the end in Line 344, and I have followed all the steps you gave, and it comes out like this....

#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<-dataX
  data<-dataY
  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
# outputX is obtained in Line 343 using data<-dataX in Line 279
plot(outputX, 
     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])

# outputY is obtained in Line 343 using data<-dataY in Line 279
plot(outputY, 
     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])

# outputZ is obtained in Line 343 using data<-dataZ in Line 279
plot(outputZ, 
     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])

############################################################

output_column1 <- cbind(outputX[,1],outputY[,1],outputZ[,1])
output_column2 <- cbind(outputX[,2],outputY[,2],outputZ[,2])

matplot(x = output_column1,
        y = output_column2,
        type = "l",
        lty = c(1, 3, 4),
        lwd = c(1, 1, 2),
        col = 1,
        tck=1,
        xlab = "Frequency (Hz)",
        ylab = "RMS Acceleration",
        xlim = c(0, SR/2),
        ylim = c(0, max(fftHistory[5:length(fftHistory)])),
        main = dataFiles[1])
legend("topright",
       legend = c("outputX", "outputY", "outputZ"),
       col = 1,
       lty = c(1, 3, 4),
       lwd = c(1, 1, 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)

Ok. Let me present the code I suggested:

#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<-dataX
  #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
outputX<-array(c(x.axis, fftHistory), dim=c(length(x.axis),2)) 

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<-dataY
  #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
outputY<-array(c(x.axis, fftHistory), dim=c(length(x.axis),2)) 

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
outputZ<-array(c(x.axis, fftHistory), dim=c(length(x.axis),2)) 

output_column1 <- cbind(outputX[,1],outputY[,1],outputZ[,1])
output_column2 <- cbind(outputX[,2],outputY[,2],outputZ[,2])

matplot(x = output_column1,
        y = output_column2,
        type = "l",
        lty = c(1, 3, 4),
        lwd = c(1, 1, 2),
        col = 1,
        tck=1,
        xlab = "Frequency (Hz)",
        ylab = "RMS Acceleration",
        xlim = c(0, SR/2),
        ylim = c(0, max(fftHistory[5:length(fftHistory)])),
        main = dataFiles[1])
legend("topright",
       legend = c("outputX", "outputY", "outputZ"),
       col = 1,
       lty = c(1, 3, 4),
       lwd = c(1, 1, 2))

#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)

Using the inputs C, GCDC and 1, this code generates the following graph:

Rplot

I guess this is what you want.

After having a look on your script, i finally figure out what you were explaining. Sorry and thank you so much for your patience and time. :pray::pray::pray::pray:

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.