Signal and time series seen from eight miles high cloud - DFT & Simple digital filtering

Discrete Fourier Transform과 단순한 주파수 도메인 필터링 

DFT_ex

Hyun Bong Lee

A simplified view of Discrete Fourier Transform, and grossly naive take on filtering.

Key Words : Fourier Series, Signal, Discrete Fourier Transform, Spectrum, Digital Sampling
Nyquist Frequency, Frequency Domain Analysis

Some say Fourier Transform is one of the 10(17) equations that changed the world(from googling). I totally agree. When I first understood its meaning, it was an amazing experience.

T = 10    # It will be the length of the Signal/time series (10 sec)
dt <- 0.01 # sampling interval. So the sampling freq. is 100Hz, & Nyquist Freq. 50Hz 
t <- seq(0,T,by=dt) # t : time indices

# 3 sinusoidal frequencies to make composite signal 
f1=3   
f2=7
f3=2  
# So, fundamental freq is 1Hz

Create complex signal with three sinusoidals plus noise

# create complex signal 
s_clean = 5*sin(2*pi*f1*t) + 2*cos(2*pi*f2*t) + 7*sin(2*pi*f3*t - pi)  # signal with no noise
noise = runif(length(t), 0, 15)*cos(2*pi*rnorm(length(t), 13, 5)*t)  # Noise amplitude isn't trivial
s <- s_clean + noise # s got noise

Plot the signals

# Plot the signals
par(mfrow=c(3,1))
plot(t, s_clean, type="l")
plot(t, noise, type="l")
plot(t, s, type="l")

Plot some temporal relationships

par(mfrow=c(3,1))
acf(s_clean)
acf(noise)  # looks good. No discernible periodic pattern 
acf(s)

par(mfrow=c(2,1))
ccf(noise, s)
ccf(s_clean, s)

Perform Fourier Transform on s. We will look into s in frequency domain.

# create frequency index
fb <- 0:length(t)/T  # Up to sampling frequency. 2*folding_freq

# Perform Fourier Transform on s
Fs <- fft(s)  # Fs = DFT(s)

# For demo ; DFT <-> IDFT relationship
Fs_inv = fft(Fs, inverse=T)/length(Fs)  # Fs_inv = IDFT(Fs) should equal s
e = (s - Fs_inv)
sum(sqrt(Re(e)^2+Im(e)^2)/length(e))  # Of course. IDFT(DFT(s)) is s  

## [1] 9.259311e-15

# Did you know that R lets one do complex number based math so easy. 
i = 0 + 1i  # Yup, i is that old i, a complex number whose square is -1. (ie. i*i = -1)
i*i  # see

## [1] -1+0i

exp(i*pi/4)  # R support this. remember, exp(i*x) = cos(x) + i*sin(x) : Euler's formula

## [1] 0.7071068+0.7071068i

cos(pi/4) + i*sin(pi/4) # see. But we don't have to do it this way.  

## [1] 0.7071068+0.7071068i

create power spectrum

# create power spectrum of s
mag <- sqrt(Re(Fs)^2+Im(Fs)^2)*2/length(Fs)
phase <- atan(Im(Fs)/Re(Fs))
par(mfrow=c(2,1))
plot(t,s,type="l",xlim=c(0,T)) ; abline(h=0,lty=3)
plot(fb, mag[1:length(fb)], type="h")

Simple filtering example

par(mfrow=c(3,1))
plot(t,s_clean,type="l",xlim=c(0,T), main="Clean signal") ; abline(h=0,lty=3)  
plot(t,s,type="l",xlim=c(0,T), main="Noisy Signal") ; abline(h=0,lty=3)
plot(fb, mag[1:length(fb)], type="h", main="Power spectrum of noisy signal")

Naive freq. based filtering. We will filter out freq. components with power spectrum below 1.5 which is primarily due to noise.

par(mfrow=c(2,1))
# simple temporal smoothing (averaging) on s
s_smooth = filter(s, rep(1/7, 7), method="convolution", sides=2)
plot(t, s_smooth, type="l",xlim=c(0,T), main="Noisy Signal smoothed on temporal domain") ; abline(h=0,lty=3)

# Simplest frequency analysis based filtering one can think of 
Fss = Fs
Fss[mag<1.5] = 0
Fss_inv = fft(Fss, inverse=T)/length(Fss)
plot(t,Re(Fss_inv),type="l",xlim=c(0,T), main="Noisy signal after frequency domain digital noise filtering")
abline(h=0,lty=3)