StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>The output of the <code>fft</code> (Fast Fourier Transform) function in R can be a little bit tedious to process. I wrote this <code>plotFFT</code> function in order to do a frequency vs power plot of the FFT. The <code>getFFTFreqs</code> function (used internally by <code>plotFFT</code>) returns the frequency associated to each FFT value. </p> <p>This was mostly based on the very interesting discussion at <a href="http://tolstoy.newcastle.edu.au/R/help/05/08/11236.html" rel="noreferrer">http://tolstoy.newcastle.edu.au/R/help/05/08/11236.html</a></p> <pre><code># Gets the frequencies returned by the FFT function getFFTFreqs <- function(Nyq.Freq, data) { if ((length(data) %% 2) == 1) # Odd number of samples { FFTFreqs <- c(seq(0, Nyq.Freq, length.out=(length(data)+1)/2), seq(-Nyq.Freq, 0, length.out=(length(data)-1)/2)) } else # Even number { FFTFreqs <- c(seq(0, Nyq.Freq, length.out=length(data)/2), seq(-Nyq.Freq, 0, length.out=length(data)/2)) } return (FFTFreqs) } # FFT plot # Params: # x,y -> the data for which we want to plot the FFT # samplingFreq -> the sampling frequency # shadeNyq -> if true the region in [0;Nyquist frequency] will be shaded # showPeriod -> if true the period will be shown on the top # Returns a list with: # freq -> the frequencies # FFT -> the FFT values # modFFT -> the modulus of the FFT plotFFT <- function(x, y, samplingFreq, shadeNyq=TRUE, showPeriod = TRUE) { Nyq.Freq <- samplingFreq/2 FFTFreqs <- getFFTFreqs(Nyq.Freq, y) FFT <- fft(y) modFFT <- Mod(FFT) FFTdata <- cbind(FFTFreqs, modFFT) plot(FFTdata[1:nrow(FFTdata)/2,], t="l", pch=20, lwd=2, cex=0.8, main="", xlab="Frequency (Hz)", ylab="Power") if (showPeriod == TRUE) { # Period axis on top a <- axis(3, lty=0, labels=FALSE) axis(3, cex.axis=0.6, labels=format(1/a, digits=2), at=a) } if (shadeNyq == TRUE) { # Gray out lower frequencies rect(0, 0, 2/max(x), max(FFTdata[,2])*2, col="gray", density=30) } ret <- list("freq"=FFTFreqs, "FFT"=FFT, "modFFT"=modFFT) return (ret) } </code></pre> <p>As an example you can try this</p> <pre><code># A sum of 3 sine waves + noise x <- seq(0, 8*pi, 0.01) sine <- sin(2*pi*5*x) + 0.5 * sin(2*pi*12*x) + 0.1*sin(2*pi*20*x) + 1.5*runif(length(x)) par(mfrow=c(2,1)) plot(x, sine, "l") res <- plotFFT(x, sine, 100) </code></pre> <p>or</p> <pre><code>linearChirp <- function(fr=0.01, k=0.01, len=100, samplingFreq=100) { x <- seq(0, len, 1/samplingFreq) chirp <- sin(2*pi*(fr+k/2*x)*x) ret <- list("x"=x, "y"=chirp) return(ret) } chirp <- linearChirp(1, .02, 100, 500) par(mfrow=c(2,1)) plot(chirp, t="l") res <- plotFFT(chirp$x, chirp$y, 500, xlim=c(0, 4)) </code></pre> <p>Which give</p> <p><a href="http://www.nicolaromano.net/misc/sine.jpg" rel="noreferrer">FFT plot of sine waves http://www.nicolaromano.net/misc/sine.jpg</a> <a href="http://www.nicolaromano.net/misc/chirp.jpg" rel="noreferrer">FFT plot of a linear chirp http://www.nicolaromano.net/misc/chirp.jpg</a></p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload