[R] Fourier Transfrom (FFT) Example

delic kranius at optonline.net
Sat Sep 26 02:53:45 CEST 2009


LOL Rolf. Yes I am sure it isn't homework. I am working on an aeroacoustics
problem and was trying to figure out how to implement a fourier transform in
R. I normally don't work in this field so this stuff was new to me at the
time of writing. I have since figured it out.

Unfortunately I don't have my actual code where I am now but here is an
older version, it might have some bugs in it since I never verified this
version. Anyway I hope it helps someone, even if it's your homework!
Apparently some don't realize that there are different ways of learning,
learning by example being one of those ways.



func<-function(x)
{
   mag2<<-mag^2
   f<<-f
	approx(f,mag2,x)$y
}


layout(matrix(c(1,2,3,4), 4, 1, byrow = TRUE))
#SETUP
   T    <- 5. #time 0 -> T 
   dt   <- 0.01 #s
   n    <- T/dt
   F    <- 1/dt # freq domain -F/2 -> F/2
   df   <- 1/T
   t    <- seq(0,T,by=dt)  
   freq <- 5 #Hz

#SIGNAL FUNCTION
   y     <- 10*sin(2*pi*freq*t) 

#FREQ ARRAY
   f <- 1:length(t)/T 

#FOURIER WORK
   Y     <- fft(y)
   mag   <- sqrt(Re(Y)^2+Im(Y)^2)*2/n #Amplitude
   phase <- Arg(Y)*180/pi 
   Yr    <- Re(Y)
   Yi    <- Im(Y)

#PLOT SIGNALS
   plot(t,y,type="l",xlim=c(0,T)) 
   grid(NULL,NULL, col = "lightgray", lty = "dotted",lwd = 1)
   par(mar=c(5, 4, 0, 2) + 0.1)
   
   plot(f[1:length(f)/2],phase[1:length(f)/2],type="l",xlab="Frequency,
Hz",ylab="Phase,deg")
   grid(NULL,NULL, col = "lightgray", lty = "dotted",lwd = 1)
   
   plot(f[1:length(f)/2],mag[1:length(f)/2],type="l",xlab="Frequency,
Hz",ylab="Amplitude")
   grid(NULL,NULL, col = "lightgray", lty = "dotted",lwd = 1) 
   
   plot(f[1:length(f)/2],(mag^2)[1:length(f)/2],type="l",xlab="Frequency,
Hz",ylab="Power, Amp^2",log="xy",ylim=c(10^-6,100))
	
	pref<-20E-6 #pa
	p<-integrate(func,f[1],f[length(f)/2])
   pwrDB<-10*log10(p$value/pref^2)
	cat("Area under power curve: ",p$value,"Pa ",pwrDB," dB\n")









Rolf Turner-3 wrote:
> 
> 
> On 17/09/2009, at 3:39 AM, delic wrote:
> 
>>
>> I wrote a script that I anticipating seeing a spike at  10Hz with the
>> function 10* sin(2*pi*10*t).
>> I can't figure out why my plots  do not show spikes at the  
>> frequencies I
>> expect. Am I doing something wrong or is my expectations wrong?
> 
> (a) Is this a homework question?
> 
> (b) Have you figured it out yet?
> 
> (c) Hint:  You have spikes at +/- 40 in a range from -50 to 50. You  
> *want* spikes
> at 10 and 90 Hz. Could it be that you haven't set your frequency  
> vector ``f''
> quite right? :-)
> 
> 	cheers,
> 
> 		Rolf Turner
> 
> P. S. You won't get spikes bang on at 10 and 90 Hz. because these are  
> *not*
> Fourier frequencies when n = 256.  If you want spikes in your  
> periodogram
> at bang on 10 and 90 Hz use a value of n that is divisible by 10,  
> e.g. n=500.
> Why would you want a power of 2 anyhow?  (Well, the fft goes faster  
> when n
> is a power of 2, but who cares?)
> 
> 		R. T.
> 
> ######################################################################
> Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 

-- 
View this message in context: http://www.nabble.com/Fourier-Transform-fft-help-tp25475063p25621211.html
Sent from the R help mailing list archive at Nabble.com.




More information about the R-help mailing list