This document contains code to perform the FFT of stimulation waveform used in the R21 project. Inputs : CSV file that contains the stimulation waveform as [time, channelA, channelC] exported from the Picoscope.

library(waved)
## 
## Attaching package: 'waved'
## 
## The following object is masked from 'package:base':
## 
##     scale
library(fftw)
library(pracma)
## Warning: package 'pracma' was built under R version 3.2.3
## 
## Attaching package: 'pracma'
## 
## The following objects are masked from 'package:waved':
## 
##     fftshift, rot90
file <- file.choose()
data <- read.csv(file,skip=1,header=T)
colnames(data) <- c('Time','ChannelA','ChannelC')

Initial plot of the recorded waveform:

plot(data$Time, data$ChannelA, type='l', xlab='Time (ms)', ylab='Stimulation Voltage (V)',col = 'blue')

The following code calculates the FFT and the frequence base from the obtained timebase. The total number of datapoints recorded in the timebase is 992067.

data_ft <- fftshift(FFT(data$ChannelA))
data$ft <- data_ft
data$ft_abs <- abs(data$ft)
timebase <- data$Time
ampdata <- data$ChannelA
omega_max <- 1/(2*(timebase[2]-timebase[1]))
omega_min <- -omega_max
freqbase <- linspace(omega_min,omega_max,992067)
len <- length(freqbase)

We use the length of the frequency base to control the display of the FFT.

plot(freqbase[(len/2-1000):(len/2+1000)]*1000, data$ft_abs[(len/2-1000):(len/2+1000)],type='l', xlab = 'Frequency in Hz', ylab = 'Magnitude FFT',col = 'brown')

A zoomed in version of the plot shows the frequency with the maximum power:

plot(freqbase[(len/2-25):(len/2+25)]*1000, data$ft_abs[(len/2-25):(len/2+25)],type='l', xlab = 'Frequency in Hz', ylab = 'Magnitude FFT', col = 'green')