Understanding data science: feature extraction with R

Keep reading for a walkthrough of how to:

- Read data into R

- Generate simple stats and plots for initial visualisation

- Perform a Fast Fourier Transform(FFT) for frequency analysis

- Caculate key features of the data

- Visualise and analyse the feature space

 사용하는 데이터는 University of Cincinnati에서 제공하는 bearing dataset이다. 이 데이터는 2,000rpm으로 회전하는 4개의 베어링을 연구하며 산출되었다. 또한 이 실험과정에서 일부 베어링이 고장이 났으며 데이터를 보는 것으로 고장을 감지해낼 수 있는지 살펴볼 것이다.

 베어링이란 무엇인가? 데이터를 분석하기 위해서는 사용할 데이터셋을 이해하는 것이 중요하다. 데이터가 나타내는 시스템이나 프로세스에 대해 지식을 가지고 있지 않으면 분석에서 유의미한 결과를 찾아내기 힘들다. 베어링이란 Rotating machine에서 핵심 구성 요소이다. 베어링의 상태가 Shaft가 얼마나 부드럽고 효율적으로 회전할 수 있을지를 결정한다. 

 1st_test.rar 파일은 2,156개의 파일을 가지고 있다. 파일은 timestamp로 네이밍되어 있으며 8 X 20,480의 매트릭스 형식으로 되어있다. 8개의 columns는 4개의 베어링이 2개씩 가지고 있는 가속도계(accelerometers)이며 rows는 1초(20kHz) 샘플이다. 한 번 데이터를 체크해보자.

setwd("C:/R/project")
getwd()

setwd를 통해 프로젝트 위치를 지정해준다. 

basedir <- "D:/Victoria/bearing_IMS/1st_test/"
data <- read.table(paste0(basedir, "2003.10.22.12.06.24"),header=FALSE,sep="\t")
head(data)

basedir을 설정해주고 이 곳에서 data를 가지고 온다. data의 일부분만 살펴보자.

> head(data)
      V1     V2     V3     V4     V5     V6     V7     V8
1 -0.022 -0.039 -0.183 -0.054 -0.105 -0.134 -0.129 -0.142
2 -0.105 -0.017 -0.164 -0.183 -0.049  0.029 -0.115 -0.122
3 -0.183 -0.098 -0.195 -0.125 -0.005 -0.007 -0.171 -0.071
4 -0.178 -0.161 -0.159 -0.178 -0.100 -0.115 -0.112 -0.078
5 -0.208 -0.129 -0.261 -0.098 -0.151 -0.205 -0.063 -0.066
6 -0.232 -0.061 -0.281 -0.125  0.046 -0.088 -0.078 -0.078

실제로 어떤 데이터를 가지고 있는지 궁금하다면 열어보는 것이 좋다.

 이름은 날짜로 되어있고 5분 주기로 1초 값이 저장된 것을 확인할 수 있다. 총 8열인데 4개의 베어링을 x,y축 기준으로 수집하였기 때문이다. 샘플링속도가 20kHz(1kHz = 1,000Hz) 라면 20,000개가 있어야 하는데 실제로는 20,480개가 있다. 이는 20kHz라고 해놓고 20.48kHz에서 하는것은 문제가 있지만 20kHz속도로 데이터가 오버해서 들어오는 것이므로 괜찮은 것으로 판별하였다. 

컬럼명이 명확하지 않으므로 적당히 수정해준다.

> colnames(data) <- c("b1.x","b1.y","b2.x","b2.y","b3.x","b3.y","b4.x","b4.y")
> head(data)
    b1.x   b1.y   b2.x   b2.y   b3.x   b3.y   b4.x   b4.y
1 -0.022 -0.039 -0.183 -0.054 -0.105 -0.134 -0.129 -0.142
2 -0.105 -0.017 -0.164 -0.183 -0.049  0.029 -0.115 -0.122
3 -0.183 -0.098 -0.195 -0.125 -0.005 -0.007 -0.171 -0.071
4 -0.178 -0.161 -0.159 -0.178 -0.100 -0.115 -0.112 -0.078
5 -0.208 -0.129 -0.261 -0.098 -0.151 -0.205 -0.063 -0.066
6 -0.232 -0.061 -0.281 -0.125  0.046 -0.088 -0.078 -0.078

summary를 통해 데이터의 간략한 정보를 보고 그래프를 그려본다.

> summary(data$b1.x)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.72000 -0.14600 -0.09500 -0.09459 -0.04200  0.38800 

> plot(data$b1.x,t="l")

이 그래프는 b1.x이므로 베어링1의 x axis이다. 평균이 0보다 작음을 알 수 있고 대부분이 -0.2 ~ 0.0에 걸쳐있다. min, max의 값이 튀긴 하지만 나머지 데이터와 다른 크기의 순서가 아니므로 나쁜 데이터라고 생각하기는 어렵다. 플룻마다 베어링의 1초단위 스냅샷을 보여주고 있으므로 데이터 분석을 수행하려면 모든 스냅샷을 함께 고려해야 한다. 이 것은 어려운 일이기 때문에 데이터 셋으로 작업할 방도를 찾아야 한다. 참고로 실제 데이터에서 다음 값을 그려낸 것이다.

 첫번째로는 전통적인 엔지니어링 방식이다. feature 추출을 통해 데이터를 압축하고 또 압축하는 것이다. 앞서 한 파일에 summary를 계산해보았다. 이처럼 평균, 최대 및 최소값, 주요 정보를 대표적으로 계산하여 압축된 데이터 모델에서 feature를 추출해낼 수 있을것이다. feature추출은 잡음이 없는 신호를 추출하여 처리해야하는 데이터의 양을 줄이는 것을 목표로 하기 때문에 원본 데이터에서 중요한 정보를 잃을 수도 있다. 한 번 extract some feautres를 진행해본다.

 분석에 유용한 형식으로 데이터를 수집하려면 2,156개의 시간순으로 정렬된 데이터를 4개의 파일로 처리한다. 베어링 분석에 대한 기존의 연구에 따르면 회전 및 스핀과 관련된 주파수에서의 진동 수준이 가장 중요하다. 하지만 기본적인 소양을 가지고 있지 않으면 주파수를 계산하기 어려우므로 데이터 패턴에 중점을 두고 분석한다. 

 진동 신호는 고속 푸리에 변환(FFT)를 사용하여 주파수 성분으로 분해할 수 있다. R은 FFT를 호출해서 사용할 수 있으며 직관적인 버전으로도 변환이 가능하다.

b1.x.fft <- fft(data$b1.x)
head(b1.x.fft)

> head(b1.x.fft)
[1] -1937.262000+ 0.000000i   -11.119305+10.116574i
[3]    -5.012298+ 7.024364i     6.090647+10.141552i
[5]    -2.122373- 2.104667i    -0.143628- 6.089364i

length(b1.x.fft)

> length(b1.x.fft)
[1] 20480

b1.x.fft[1:10240]
amplitude <- Mod(b1.x.fft[1:(length(b1.x.fft)/2)])
# mod(n, m) = n - m * floor(n/m)
frequency <- seq(0, 10000, length.out=length(b1.x.fft)/2)
plot(amplitude ~ frequency, t="l")

중간 과정이 잘 이해가 안되니 b1.x.fft[1]을 한 번 살펴보자. 

> b1.x.fft[1]
[1] -1937.262+0i
> Mod(b1.x.fft[1])
[1] 1937.262
> amplitude[1]
[1] 1937.262

이 때에 1에서 1/2까지의 데이터만 사용하는 이유는 좌우 대칭인 그래프가 그려지기 때문이다. 조건을 달리하여 그래프를 그려보자.

temp_amplitude <- Mod(b1.x.fft)
temp_frequency <- seq(length(b1.x.fft))
plot(temp_amplitude ~ temp_frequency, t="l")

이제 frequency를 볼 수 있으며 낮은 주파수에 초점을 맞춰 다시 그려본다.

plot(amplitude ~ frequency, t="l", xlim=c(0,1000),ylim=c(0,500))
axis(1, at=seq(0,1000,100),labels=FALSE)

1kHz 아래에서 가장 큰 주파수가 보인다. 가장 큰 frequency를 보이는 상위 주파수 15개를 추출해본다.

sorted <- sort.int(amplitude, decreasing=TRUE, index.return=TRUE)
top15 <- sorted$ix[1:15] # indexes of the largest 15
top15f <- frequency[top15] # convert indexes to frequencies
top15
top15f

 

> top15
 [1]    1 1011 1018  506 1004 1019  994  996   60 1003  945   52
[13] 4527 3694 4432
> top15f
 [1]    0.00000  986.42446  993.26106  493.21223  979.58785
 [6]  994.23772  969.82127  971.77459   57.62281  978.61119
[11]  921.96504   49.80955 4420.35355 3606.79754 4327.57105

이를 통해 알 수 있는 이상치는 49.8Hz, 57.6Hz, 493Hz이다. 5개 구성 요소의 빈도를 중점적으로 살펴본다. x축과 y축이 존재하기 때문에 10개의 feature를 가질 것이며 이를 function으로 만들면 다음과 같다.

fft.profile <- function (dataset, n)
{
  fft.data <- fft(dataset)
  amplitude <- Mod(fft.data[1:(length(fft.data)/2)])
  frequencies <- seq(0, 10000, length.out=length(fft.data)/2)
  sorted <- sort.int(amplitude, decreasing=TRUE, index.return=TRUE)
  top <- sorted$ix[1:n] # indexes of the largest n components
  return (frequencies[top]) # convert indexes to frequencies
}

# How many FFT components should I grab as features?
n <- 5
# Set up storage for bearing-grouped data
b1 <- matrix(nrow=0, ncol=(2*n+1))
b2 <- matrix(nrow=0, ncol=(2*n+1))
b3 <- matrix(nrow=0, ncol=(2*n+1))
b4 <- matrix(nrow=0, ncol=(2*n+1))
for (filename in list.files(basedir))
{
cat("Processing file ", filename, "\n")
timestamp <- as.character(strptime(filename, format="%Y.%m.%d.%H.%M.%S"))
data <- read.table(paste0(basedir, filename), header=FALSE, sep="\t")
colnames(data) <- c("b1.x", "b1.y", "b2.x", "b2.y", "b3.x", "b3.y", "b4.x", "b4.y")
# Bind the new rows to the bearing matrices
b1 <- rbind(b1, c(timestamp, fft.profile(data$b1.x, n), fft.profile(data$b1.y, n)))
b2 <- rbind(b2, c(timestamp, fft.profile(data$b2.x, n), fft.profile(data$b2.y, n)))
b3 <- rbind(b3, c(timestamp, fft.profile(data$b3.x, n), fft.profile(data$b3.y, n)))
b4 <- rbind(b4, c(timestamp, fft.profile(data$b4.x, n), fft.profile(data$b4.y, n)))
}
write.table(b1, file=paste0(basedir, "../b1.csv"), sep=",", row.names=FALSE, col.names=FALSE
write.table(b2, file=paste0(basedir, "../b2.csv"), sep=",", row.names=FALSE, col.names=FALSE
write.table(b3, file=paste0(basedir, "../b3.csv"), sep=",", row.names=FALSE, col.names=FALSE
write.table(b4, file=paste0(basedir, "../b4.csv"), sep=",", row.names=FALSE, col.names=FALSE

 4개의 베어링을 FFT 변환을 하여 상위 5개의 feature를 저장하는 코드이다.

b1.csv에는 다음과 같은 데이터들이 저장된다. 그러니까 b1의 상위 5개 feature가 시간일자로 정리되어 있는 것이다. 그래프로 그릴 때에는 b1.x의 그래프를 주르륵~ 그려내는 것이다.

 결과적으로 나타나는 데이터는 2156(측정 데이터) row , datetime + x.feature(5) + y.feature(5) = column 으로 이루어져있다. 값은 frequency(kHz)이다. 

# x axis components
plot(b1[, 2], pch=4, col="dodgerblue2", ylab="Frequency", main="1st strongest x", ylim=c(0, 2)) 
points(b2[, 2], pch=3, col="darkorchid2")
points(b4[, 2], pch=1, col="chartreuse")
points(b3[, 2], pch=20, col="coral2")

legend("topleft", c("bearing 1", "bearing 2", "bearing 3", "bearing 4"), col=c("dodgerblue2", "darkorchid2", "coral2", "chartreuse"), pch=c(4, 3, 20, 1))

plot(b1[, 7], pch=4, col="dodgerblue2", ylab="Frequency", main="1st strongest y", ylim=c(0, 2)) 
points(b2[, 7], pch=3, col="darkorchid2")
points(b4[, 7], pch=1, col="chartreuse")
points(b3[, 7], pch=20, col="coral2")

legend("topleft", c("bearing 1", "bearing 2", "bearing 3", "bearing 4"), col=c("dodgerblue2", "darkorchid2", "coral2", "chartreuse"), pch=c(4, 3, 20, 1))


plot(b1[, 3], pch=4, col="dodgerblue2", ylab="Frequency", main="2nd strongest x", ylim=c(0,10000))
points(b2[, 3], pch=3, col="darkorchid2")
points(b4[, 3], pch=1, col="chartreuse")
points(b3[, 3], pch=20, col="coral2")

plot(b1[, 8], pch=4, col="dodgerblue2", ylab="Frequency", main="2nd strongest y", ylim=c(0,10000))
points(b2[, 8], pch=3, col="darkorchid2")
points(b4[, 8], pch=1, col="chartreuse")
points(b3[, 8], pch=20, col="coral2")

plot(b1[, 4], pch=4, col="dodgerblue2", ylab="Frequency", main="3rd strongest x", ylim=c(0,10000))
points(b2[, 4], pch=3, col="darkorchid2")
points(b4[, 4], pch=1, col="chartreuse")
points(b3[, 4], pch=20, col="coral2")

plot(b1[, 9], pch=4, col="dodgerblue2", ylab="Frequency", main="3rd strongest y", ylim=c(0,10000))
points(b2[, 9], pch=3, col="darkorchid2")
points(b4[, 9], pch=1, col="chartreuse")
points(b3[, 9], pch=20, col="coral2")

plot(b1[, 5], pch=4, col="dodgerblue2", ylab="Frequency", main="4th strongest x", ylim=c(0,10000))
points(b2[, 5], pch=3, col="darkorchid2")
points(b4[, 5], pch=1, col="chartreuse")
points(b3[, 5], pch=20, col="coral2")

plot(b1[, 10], pch=4, col="dodgerblue2", ylab="Frequency", main="4th strongest y", ylim=c(0,10000))
points(b2[, 10], pch=3, col="darkorchid2")
points(b4[, 10], pch=1, col="chartreuse")
points(b3[, 10], pch=20, col="coral2")

plot(b1[, 6], pch=4, col="dodgerblue2", ylab="Frequency", main="5th strongest x", ylim=c(0,10000))
points(b2[, 6], pch=3, col="darkorchid2")
points(b4[, 6], pch=1, col="chartreuse")
points(b3[, 6], pch=20, col="coral2")

plot(b1[, 11], pch=4, col="dodgerblue2", ylab="Frequency", main="5th strongest y", ylim=c(0,10000))
points(b2[, 11], pch=3, col="darkorchid2")
points(b4[, 11], pch=1, col="chartreuse")
points(b3[, 11], pch=20, col="coral2")

헷갈려서 적어보자면 R은 index를 집어 넣어주는지 데이터[0]은 row 순서이다. 1이 date이고 2는.. 보시다시피... 토대로 추출한 feature의 그래프를 그려본다. 

 feature 그래프를 보면 어느 순간부터 문제점을 보이기 시작한다는 것을 알 수 있다. 안타깝게도 FFT는 '어떤 시점'에서 문제가 생겼는지는 알 수 없다. 그래도 베어링 3번 4번이 문제가 생겼다는 것을 어림짐작으로 알 수 있다.

 그래프를 그리고 패턴을 알아내는 것은 데이터 과학의 핵심 부분이다. 하지만 데이터의 양이 많은 경우 이 것이 어려울 수 있다. 그렇기 때문에 피쳐 추출 기술을 이용하여 데이터를 축소하였다. 이를 발전 시키기 위해서는 베어링 결함을 모델링하거나, 베어링 결함을 감지하게 하는 것이 있다. 

---

알아야 할 지식

- 푸리에 변환

- 차원축소와 특징선택(http://terryum.io/korean/2016/05/05/FeatureSelection_KOR/)

 머신러닝의 성능은 데이터의 양과 질에 굉장히 의존적이다. 때문에 충분한 데이터를 모으고, 어떤 feature가 유용한지 아닌지 확인하는 과정을 거친다. 이 것을 확인하는 과정을 특징 선택(특징 추출)이라고 부른다.

---

사용한 전체 코드

basedir <- "Victoria/bearing_IMS/1st_test/"
data <- read.table(paste0(basedir, "2003.10.22.12.06.24"), header=FALSE, sep="\t")
head(data)
colnames(data) <- c("b1.x", "b1.y", "b2.x", "b2.y", "b3.x", "b3.y", "b4.x", "b4.y")

summary(data$b1.x)
plot(data$b1.x, t="l") # t="l" means line plot

b1.x.fft <- fft(data$b1.x)
amplitude <- Mod(b1.x.fft[1:(length(b1.x.fft)/2)])
frequency <- seq(0, 10000, length.out=length(b1.x.fft)/2)
plot(amplitude ~ frequency, t="l")

plot(amplitude ~ frequency, t="l", xlim=c(0,1000), ylim=c(0,500))
axis(1, at=seq(0,1000,100), labels=FALSE)  # add more ticks

sorted <- sort.int(amplitude, decreasing=TRUE, index.return=TRUE)
top15 <- sorted$ix[1:15] # indexes of the largest 15
top15f <- frequency[top15] # convert indexes to frequencies

fft.profile <- function (dataset, n)
{
  fft.data <- fft(dataset)
  amplitude <- Mod(fft.data[1:(length(fft.data)/2)])
  frequencies <- seq(0, 10000, length.out=length(fft.data)/2)
  
  sorted <- sort.int(amplitude, decreasing=TRUE, index.return=TRUE)
  top <- sorted$ix[1:n] # indexes of the largest n components
  return (frequencies[top]) # convert indexes to frequencies
}

# How many FFT components should I grab as features?
n <- 5

# Set up storage for bearing-grouped data
b1 <- matrix(nrow=0, ncol=(2*n+1))
b2 <- matrix(nrow=0, ncol=(2*n+1))
b3 <- matrix(nrow=0, ncol=(2*n+1))
b4 <- matrix(nrow=0, ncol=(2*n+1))

for (filename in list.files(basedir))
{
  cat("Processing file ", filename, "\n")
  
  timestamp <- as.character(strptime(filename, format="%Y.%m.%d.%H.%M.%S"))
  
  data <- read.table(paste0(basedir, filename), header=FALSE, sep="\t")
  colnames(data) <- c("b1.x", "b1.y", "b2.x", "b2.y", "b3.x", "b3.y", "b4.x", "b4.y")
  
  # Bind the new rows to the bearing matrices
  b1 <- rbind(b1, c(timestamp, fft.profile(data$b1.x, n), fft.profile(data$b1.y, n)))
  b2 <- rbind(b2, c(timestamp, fft.profile(data$b2.x, n), fft.profile(data$b2.y, n)))
  b3 <- rbind(b3, c(timestamp, fft.profile(data$b3.x, n), fft.profile(data$b3.y, n)))
  b4 <- rbind(b4, c(timestamp, fft.profile(data$b4.x, n), fft.profile(data$b4.y, n)))
  
}

write.table(b1, file=paste0(basedir, "../b1.csv"), sep=",", row.names=FALSE, col.names=FALSE)
write.table(b2, file=paste0(basedir, "../b2.csv"), sep=",", row.names=FALSE, col.names=FALSE)
write.table(b3, file=paste0(basedir, "../b3.csv"), sep=",", row.names=FALSE, col.names=FALSE)
write.table(b4, file=paste0(basedir, "../b4.csv"), sep=",", row.names=FALSE, col.names=FALSE)

par(mfrow = c(5,2))
# x axis components
plot(b1[, 2], pch=4, col="dodgerblue2", ylab="Frequency", main="1st strongest x", ylim=c(0, 2)) 
points(b2[, 2], pch=3, col="darkorchid2")
points(b4[, 2], pch=1, col="chartreuse")
points(b3[, 2], pch=20, col="coral2")

legend("topleft", c("bearing 1", "bearing 2", "bearing 3", "bearing 4"), col=c("dodgerblue2", "darkorchid2", "coral2", "chartreuse"), pch=c(4, 3, 20, 1))

plot(b1[, 7], pch=4, col="dodgerblue2", ylab="Frequency", main="1st strongest y", ylim=c(0, 2)) 
points(b2[, 7], pch=3, col="darkorchid2")
points(b4[, 7], pch=1, col="chartreuse")
points(b3[, 7], pch=20, col="coral2")

legend("topleft", c("bearing 1", "bearing 2", "bearing 3", "bearing 4"), col=c("dodgerblue2", "darkorchid2", "coral2", "chartreuse"), pch=c(4, 3, 20, 1))

plot(b1[, 3], pch=4, col="dodgerblue2", ylab="Frequency", main="2nd strongest x", ylim=c(0,10000))
points(b2[, 3], pch=3, col="darkorchid2")
points(b4[, 3], pch=1, col="chartreuse")
points(b3[, 3], pch=20, col="coral2")

plot(b1[, 8], pch=4, col="dodgerblue2", ylab="Frequency", main="2nd strongest y", ylim=c(0,10000))
points(b2[, 8], pch=3, col="darkorchid2")
points(b4[, 8], pch=1, col="chartreuse")
points(b3[, 8], pch=20, col="coral2")

plot(b1[, 4], pch=4, col="dodgerblue2", ylab="Frequency", main="3rd strongest x", ylim=c(0,10000))
points(b2[, 4], pch=3, col="darkorchid2")
points(b4[, 4], pch=1, col="chartreuse")
points(b3[, 4], pch=20, col="coral2")

plot(b1[, 9], pch=4, col="dodgerblue2", ylab="Frequency", main="3rd strongest y", ylim=c(0,10000))
points(b2[, 9], pch=3, col="darkorchid2")
points(b4[, 9], pch=1, col="chartreuse")
points(b3[, 9], pch=20, col="coral2")

plot(b1[, 5], pch=4, col="dodgerblue2", ylab="Frequency", main="4th strongest x", ylim=c(0,10000))
points(b2[, 5], pch=3, col="darkorchid2")
points(b4[, 5], pch=1, col="chartreuse")
points(b3[, 5], pch=20, col="coral2")

plot(b1[, 10], pch=4, col="dodgerblue2", ylab="Frequency", main="4th strongest y", ylim=c(0,10000))
points(b2[, 10], pch=3, col="darkorchid2")
points(b4[, 10], pch=1, col="chartreuse")
points(b3[, 10], pch=20, col="coral2")

plot(b1[, 6], pch=4, col="dodgerblue2", ylab="Frequency", main="5th strongest x", ylim=c(0,10000))
points(b2[, 6], pch=3, col="darkorchid2")
points(b4[, 6], pch=1, col="chartreuse")
points(b3[, 6], pch=20, col="coral2")

plot(b1[, 11], pch=4, col="dodgerblue2", ylab="Frequency", main="5th strongest y", ylim=c(0,10000))
points(b2[, 11], pch=3, col="darkorchid2")
points(b4[, 11], pch=1, col="chartreuse")
points(b3[, 11], pch=20, col="coral2")