- 資料
- 資料のソースをRに書き換えるとfft()の結果と同様の結果になる
x<-seq(from=0,to=2*pi,length.out=100)
z<-cos(2*x)+2*cos(5*x)
par(mfcol=c(2,2))
plot(z,type="l")
y<-fft(z)
plot(Re(y),type="l")
plot(Re(y[1:10]),type="l")
a<-fft(y, inverse = TRUE) / length(y)
plot(Re(a),type="l")
x<-seq(from=0,to=2*pi,length.out=100)
z<-cos(2*x)+2*cos(5*x)+rnorm(length(x))
Rs<-c()
Is<-c()
for(i in 1:length(z)){
tmpR<-0
tmpI<-0
for(j in 1:length(z)){
tmpR<-tmpR+z[j]*cos(2*pi*(j-1)*(i-1)/length(z))
tmpI<-tmpI+z[j]*sin(2*pi*(j-1)*(i-1)/length(z))
}
Rs<-c(Rs,tmpR)
Is<-c(Is,tmpI)
}
y<-fft(z)
a<-fft(y, inverse = TRUE) / length(y)
par(mfcol=c(1,2))
plot(Rs,ylim=c(min(Rs,Re(y)),max(Rs,Re(y))))
par(new=TRUE)
plot(Re(y),ylim=c(min(Rs,Re(y)),max(Rs,Re(y))),col=2,type="l")
z2<-rep(0,length(z))
for(i in 1:length(Rs)){
for(j in 1:length(Rs)){
z2[i]<-z2[i]+Rs[j]*cos(2*pi*(j-1)*(i-1)/length(z))+Is[j]*sin(2*pi*(j-1)*(i-1)/length(z))
}
}
z2<-z2/(length(z2))
plot(z,ylim=c(min(z,z2),max(z,z2)))
par(new=TRUE)
plot(z2,type="l",col=2,lwd=3,ylim=c(min(z,z2),max(z,z2)))
par(new=TRUE)
a<-fft(y, inverse = TRUE) / length(y)
plot(Re(a),type="l",col="dark blue",lwd=1,ylim=c(min(z,z2),max(z,z2)))
par(mfcol=c(1,1))