離散三角フーリエ変換とその逆変換

  • 資料
  • 資料のソースを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")
#plot(Is,ylim=c(min(Is,Im(y)),max(Is,Im(y))))
#par(new=TRUE)
#plot(Im(y),ylim=c(min(Is,Im(y)),max(Is,Im(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))