高速フーリエ変換

20140910_01
#合成波仕様
a<-c(40,2,4,6,8,10,12) #amplitude
frequency<-c(10,2,4,6,8,50,60)
omega<-2*pi*frequency #Angular Velocity
alpha<-pi/180*c(35,90,-15,-60,80,30,-45) #phase
#サンプリング仕様
SamplingPoint<-4096
n<-0:(SamplingPoint-1)
SamplingFrequency<-1000 #(Hz)
t<-n/SamplingFrequency #(s)
f<-n*SamplingFrequency/SamplingPoint
#合成波
SyntheticWave<-numeric(SamplingPoint)
for(iii in 0:(length(n)-1)){
wave<-a*sin(omega*t[iii]+alpha)
SyntheticWave[iii]<-sum(wave)
}
tmp<-fft(SyntheticWave)
Spectrum<-Re(tmp)^2+Im(tmp)^2 #or abs(fft(SyntheticWave))^2
par(mfrow=c(3,1),mar=c(3,3,3,3),ps=15,cex.axis=1,cex.lab=1,cex.main=1,cex.sub=1)
plot(t,SyntheticWave,type="l")
plot(f,Spectrum,type="l",xlim=c(0.1,SamplingFrequency/2),log="xy") #NyquistFrequency
plot(1/f,Spectrum,type="l",xlim=c(0.1,1/(SamplingFrequency/2)),log="xy")
title(paste("\n","SamplingPoint=",SamplingPoint," SamplingFrequency=",SamplingFrequency,
" t=",1/SamplingFrequency," f=",SamplingFrequency/SamplingPoint
," Mean=",mean(SyntheticWave),sep=""),outer=T)