[R] Fourier Transfrom (FFT) Example

2009-09-25 Thread delic

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@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.

__
R-help@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.


Re: [R] Fourier Transform fft help

2009-09-16 Thread delic

After extensive research I found some information hinting at my freq array
being wrong. I thought that the data spanned -F/2 to F/2. It seems that I
should only plot the first half of the fft vector? Can anyone confirm or
shed some light on this?

One other thing I thought the magnitude should be 10 for 10* sin(2*pi*10*t),
my magnitudes are much higher.



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?
> 
> 
> require(stats)
> layout(matrix(c(1,2,3), 3, 1, byrow = TRUE))
> #SETUP
>n<- 256 #nextn(1001) gives next power 2
>F<- 100 #Hz -50 to 50 Hz
>dt   <- 1/F
>T<- n*dt
>df   <- 1/T
>t<- seq(0,T,by=dt)  #also try ts function
>t<-t[1:length(t)-1]
>freq <- 5 #Hz
> 
> #SIGNAL FUNCTION
>y <- 10* sin(2*pi*10*t)   #10*sin(2*pi*freq*t) 
> #FREQ ARRAY
>#f <- seq(-F/2,F/2,by=df)
>   f <- F/2*seq(-1,1,length.out=n+1)
>f<-f[1:length(f)-1]
> #FOURIER WORK
>Y <- fft(y)
>mag   <- sqrt(Re(Y)^2+Im(Y)^2)
>phase <- atan(Im(Y)/Re(Y))
>Yr<- Re(Y)
>Yi<- Im(Y)
> 
> #PLOT SIGNALS
>plot(t,y,type="l",xlim=c(0,T)) 
>plot(f,Re(Y),type="l")
>plot(f,Im(Y),type="l")
> 

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

__
R-help@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.


[R] Fourier Transform fft help

2009-09-16 Thread delic

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?


require(stats)
layout(matrix(c(1,2,3), 3, 1, byrow = TRUE))
#SETUP
   n<- 256 #nextn(1001) gives next power 2
   F<- 100 #Hz -50 to 50 Hz
   dt   <- 1/F
   T<- n*dt
   df   <- 1/T
   t<- seq(0,T,by=dt)  #also try ts function
   t<-t[1:length(t)-1]
   freq <- 5 #Hz

#SIGNAL FUNCTION
   y <- 10* sin(2*pi*10*t)   #10*sin(2*pi*freq*t) 
#FREQ ARRAY
   #f <- seq(-F/2,F/2,by=df)
f <- F/2*seq(-1,1,length.out=n+1)
   f<-f[1:length(f)-1]
#FOURIER WORK
   Y <- fft(y)
   mag   <- sqrt(Re(Y)^2+Im(Y)^2)
   phase <- atan(Im(Y)/Re(Y))
   Yr<- Re(Y)
   Yi<- Im(Y)

#PLOT SIGNALS
   plot(t,y,type="l",xlim=c(0,T)) 
   plot(f,Re(Y),type="l")
   plot(f,Im(Y),type="l")
-- 
View this message in context: 
http://www.nabble.com/Fourier-Transform-fft-help-tp25475063p25475063.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@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.