[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