マコリンとかまいんちゃんとかNHKは見るべきなのにスクランブルされてもいいとかイミフ

つーかマコリンかわいすぎ☆彡 スピンオフで番組始まって欲しいです.

あ, どーも僕です.

今クールのドラマはNHKが圧勝ですね.

・シングルマザー

・トトリ

・薄桜記

どれも面白いです. はい. ハエ女は保留ですかね.

で, こないだ離散フーリエ変換のレポートがでましたので Rでコード組んでみました.

まぁわざわざRでやる必要もないのですが・・・

つーかまぁRにはフーリエ変換の関数はあるのですが・・・

そのへんは気にせずにどんな感じだったのか載せたいとおもいます.

離散フーリエ変換の式はこちら

フーリエ変換って結局, パワースペクトルのグラフ作りが目的です. パワースペクトルを知れば どの周波数がどんだけの強さで混じっているのかわかるって話です. パワースペクトルを計算するにはフーリエ変換したものの実部と虚部があればいいのです. それらの式はこちら. まぁオイラー式を使うだけですが.

パワースペクトルの式はこちら. これは両側スペクトル.

片側スペクトル.

片側スペクトルは両側スペクトルの２倍ってだけですね.

レポートでは先生からもらった仮想データ(トレンド除去済み)の両側スペクトルと片側スペクトルを 書いてくるって話です. Rでやってみたら以下のようになりました. まずはデータと パワースペクトルを計算する関数の定義.

#****** データ　*********
N <- 36
n <- 0:35
k <- 0:36
dt <- 1
dat <- c( 9.6, -1.8, -10.9, -12.0,  -6.3,   2.0,  10.5,  20.8,  4.1,
          0.9,  6.9,  14.8,  -9.0, -16.0, -15.0, -12.1,  -8.9, -0.6,
         11.1, 17.5,   5.2,   6.1,   2.8,   0.7,  -9.4, -14.4, -6.6,
        -14.8, -0.8,  -1.9,   4.0,  18.5,  10.6,  -5.8,   2.0, -1.8)

# Nyquist frequency
Nyq.freq <- 1 / ( 2 * dt )

# fundamental frequency △f
T <- dt * N
delta.freq <- 1 / T


#******* 離散フーリエ変換の実部と虚部を計算　*********
for.real <- function(x)  dat * cos(2*x*pi*n/N)
Re       <- (colSums(sapply(k, for.real))/sqrt(N))**2
for.imag <- function(x)  dat * sin(2*x*pi*n/N)
Im       <- (colSums(sapply(k, for.imag))/sqrt(N))**2

# 2 - side
S <- function(x){
  x <- x + 1
  rtn <- (Re[x] + Im[x])/N
  return(rtn)
}
# 1 - side
G <- function(x){
  if(x %in% c(0, N/2)) return(S(x))
  return(2 * S(x))
}
#******************************************************

2-sideの結果とグラフ作り.

two_side.ps <- c(sapply(18:35, S), sapply(0:18, S))
def.par <- par(no.readonly = TRUE) # デフォルトの環境を保存
par(mai = c(1.5, 1.5, 0.5, 0.5))
plot(two_side.ps, type = "n", tcl = 0.5,
     ylab = "Power Spectrum S(k)", xlab = "k",
     cex.axis = 1.4, cex.lab = 1.4, xaxt = "n")
box()
axis(side = 1, at = seq(2,36,2), labels = seq(-17, 17, 2), tcl = 0.5)
points(two_side.ps, type = "b", lwd = 2.5, col = "#0068b7",
       pch = 21)
par(def.par)　　　　　　　　　　　 # デフォルトの環境を復帰

1-sideの結果とグラフ作り. あっそういえばこのグラフの片側スペクトルは周波数領域に変換した半連続 なものです. 簡単なので式は書きませんが,　つまりはGにTを乗じているだけです.

xaxis.label <- as.character(round(0:18*delta.freq, 2))
G.p <- function(x) G(x)*T　　　　　　　　　　　　　　#これで半連続に！
one_side.ps <- sapply(0:18, G.p)
def.par <- par(no.readonly = TRUE) # デフォルトの環境を保存
par(mai = c(1.5, 1.5, 0.5, 0.5))
plot(one_side.ps, type = "n", tcl = 0.5,
     ylab = expression(paste("Power Spectrum ", G^",","(",f[k],")")), 
     xlab = expression(paste("Frequency | ",f[k])), 
     cex.axis = 1.4, cex.lab = 1.4, xaxt = "n")
box()
axis(side = 1, at = seq(1,19,3),  tcl = 0.5,
     labels = xaxis.label[seq(1,19,3)])
points(one_side.ps, type = "l", lwd = 2.5, col = "#0068b7",
       pch = 21)
par(def.par)　　　　　　　　　　　 # デフォルトの環境を復帰

よかった!!!!それっぽい図になっていますね!!!!

フーリエ変換はちゃんと計算できるとデータをの自乗平均と片側スペクトルGの和が 同じになるので計算してみました.

mean(dat**2)
[1] 99.80944
sum(sapply(0:18, G))
[1] 99.80944

同じだぜ( • ̀ω•́ )✧

最後にこれはレポート関係ないのですが逆離散フーリエ変換して元に戻る のか試してみました. レポートと関係ないのでグラフはテキトーですが・・・.

逆離散フーリエ変換の式はこちら.(最初に書いたののデータが実数のときバージョン）. あと計算しやすいようにフーリエ変換の定義を変えてます・・・.

Rのコードと結果はこちら

#*********************************************************
#おまけ 逆離散フーリエ変換してみる
k <- 0:35
for.real  <- function(x)  dat * cos(2*x*pi*n/N)
Re        <- colSums(sapply(k, for.real))
Re.to.inv <- sapply(k, function(k) Re[k+1]*cos(2*k*pi*n/N))
for.imag  <- function(x)  dat * sin(2*x*pi*n/N)
Im        <- -colSums(sapply(k, for.imag))
Im.to.inv <- sapply(k, function(k) Im[k+1]*sin(2*k*pi*n/N))

inv.x <- rowMeans(Re.to.inv - Im.to.inv)
plot(dat, type = "b", lwd = 2, pch = 21)
points(inv.x, type = "h", col = "blue", lwd = 2)
#******************************************************

ばっちしやね！