Summary

엑셀을 이용한 FFT 분석 예제가 다음과 같이 구현되었다.

https://blog.naver.com/skkong89/90102712165

R을 이용한 신호처리 라이브러리등이 있지만, 엑셀을 이용한 방식과 동일한 순서로,이번에는 R 을 이용해서 결과를 확인한다. 또한 10Hz, 25Hz 의 합성함수일 경우와 10Hz 에 noise 가 추가되었을 경우에도 DFT 변환 (Discrete Fourier Transform) 을 통한 주파수 정보가 정상적으로 구해지는지 테스트한다.

library(ggplot2) # 그래프
library(mrbsizeR) # fft 
## Loading required package: maps
N <- 256 # 샘플 수, 2^8
T <- pi # 전체 샘플링 시간 T = 0 ~ pi(3.14) sec
fs <- T/N # 샘플링 주파수(초당 샘플수)
fr <- 1/T # 주파수 해상도
tr <- T / N # 샘플 사이의 시간 증가분

f <- 10 # 주파수 Hz

# 0~pi sec 사이에 256개 샘플링 데이터 생성함
x_time <- seq(from = 0, to = T, by=tr) # 257개가 생성된다. 256개로 맞춰준다.
x_time <- x_time[1:N]
head(x_time) # 엑셀에서 1. Time 컬럼과 동일한 값 출력
## [1] 0.00000000 0.01227185 0.02454369 0.03681554 0.04908739 0.06135923
# 샘플 데이터 생성: 10Hz로 먼저 테스트해보고, (10Hz, 25Hz) 합성함수에 대해서도 테스트 해본다.
sampling_df <- sin(2*pi*f*x_time) # 샘플 데이터를 생성한다.
#sampling_df <- sin(2*pi*f*x_time) + rnorm(N, mean=0, sd=0.05)  # 샘플 데이터 + noise
#sampling_df <- sin(2*pi*f*x_time) + sin(2*pi*25*x_time) # 10Hz, 25Hz
head(sampling_df) # 엑셀에서 1. Data 컬럼과 동일한 값 출력(10Hz 일 경우)
## [1]  0.00000000  0.69689787  0.99958903  0.73685341  0.05730986 -0.65465154
# 위상의 변화를 주었을 경우도 정상 동작하는지 확인한다.
#test1 <- c(sampling_df[201:256], sampling_df[1:200])
#sampling_df <- test1


# FFT freq: 샘플링 시간에 따른 주파수 해상도 정보를 구한다. 
# 0부터 시작해서 주파수 해상도만큼 더해간다.
fft_freq <- seq(from = 0, by = fr, length.out = 256)
head(fft_freq) # 엑셀에서 1. FFT freq 컬럼과 동일한 값 출력
## [1] 0.0000000 0.3183099 0.6366198 0.9549297 1.2732395 1.5915494
# 0~0.5 sec 사이의 그래프를 그려본다.
g <- ggplot(data = data.frame(x = x_time, y = sampling_df), aes(x = x, y = y)) + geom_line() 
g <- g + xlim(0, 0.5) + xlab('Time(sec)') + ylab('Amplitude') + ggtitle('Data in Time Domain')
g

# DFT 행렬을 생성한다.
D <- dftMatrix(N) # 샘플링 갯수 만큼의 DFT matrix 생성

# DFT (Discrete Fourier Transform) 을 수행한다.
coeff <- D %*% sampling_df # coeff: Fourier coefficients
head(coeff) # 엑셀에서 2. FFT complex 컬럼과 동일한 값 출력
##                    [,1]
## [1,] 2.044075+0.000000i
## [2,] 2.046523-0.021890i
## [3,] 2.053896-0.043908i
## [4,] 2.066284-0.066183i
## [5,] 2.083844-0.088849i
## [6,] 2.106796-0.112052i
# FFT mag(진폭 크기): 위에서 구한 FFT 계수 정보에 대한 크기를 구한다.
fft_mag <- Mod(coeff) / (N/2) 
head(fft_mag) # 엑셀에서 3. FFT mag 컬럼과 동일한 값 출력
##            [,1]
## [1,] 0.01596934
## [2,] 0.01598937
## [3,] 0.01604973
## [4,] 0.01615113
## [5,] 0.01629482
## [6,] 0.01648261
# frequency 도메인에서의 주파수를 그린다.
g <- ggplot(data = data.frame(x = fft_freq, y = fft_mag), aes(x = x, y = y)) + geom_line() + xlim(0, 40)
g <- g + xlab('frequency(Hz)') + ylab('amplitude') + ggtitle('Data in Frequency Domain')
g