Rで回路設計について調べるのは、検索性が悪いと思うのですよねー…
2018/04/04 追記
自分で作ってみました。
niszet.hatenablog.com
電気回路ではRという文字を抵抗を表すシンボルに使うため、うまいこと検索でひっかけられないのでした。
最近はR lang という形で"lang"と入れるようにして検索してます。そうしたら見つけてしまった次の記事：
LC LOW-PASS FILTER ANALYSIS IN R
paulorenato.com
これ、Rでやろうと思っていた小ネタの一つだったのですが、ggplot2使っている、complex型使っている、で完全に被ってました。むむむ。
ちょっとだけ解説というか、自分なりの（無駄）アレンジを入れてみます。
ちょっとだけ解説
Rには便利な型が色々とありますね。しかし、仲間はずれがいます。rawとcomplex型です。Rについては結構色々本を読んでいると思うものの、この二つについての解説はほとんど見たことがない気がします。曰く「あまり使わない」…
まぁ、普通の使い方をしていて使うシーンはあまりないとは思うのですが。
自分はraw型はWaveのファイルを読むために既に使っていますが、complex型はまだ使ってなかったです。
rawはバイナリファイルを読み書きするには必須ですが、complex型は複素数。統計分野で虚数なんて出てくるのでしょうか？使われなさそうですね。
電気電子回路や制御の分野だとフィルタの設計などでよく使うのです。最初にあげたサイトもLCフィルタについての記事となっています。
さて、ちょっと使ってみましょうか…。
まずは普通にplotします
plot(exp(-1i*seq(0,1,0.01)*2*pi))
plotはdata frameにしなくても表示できるので、ちゃちゃっと確認する際によく使いますね。
結果は下図のとおり。
これは横軸に実数部、縦軸に虚数部が描かれていますね。
これを明示的に行うと、
data.frame(Re(exp(-1i*seq(0,1,0.01)*2*pi)),Im(exp(-1i*seq(0,1,0.01)*2*pi))) %>% plot()
で、同じ図が得られる。
つづいて、ggplotで同じことをすると、
library(tidyverse)
data.frame(re_x=Re(exp(-1i*seq(0,1,0.01)*2*pi)),im_x=Im(exp(-1i*seq(0,1,0.01)*2*pi))) %>% ggplot(aes(x=re_x,y=im_x))+geom_point()
で表示できます。
さて、ここまでは実部と虚部をReとImを使用してそれぞれプロットしていたので、complex型をそのままy軸にプロットしているわけではない点に注意です。横軸が実数、縦軸が虚数という形。
y軸側をcomplexにしたらどうなるのでしょう？plotで描くなら下記のようになりますか
data_frame(x=seq(0,1,0.01),y=exp(-1i*seq(0,1,0.01)*2*pi)) %>% plot()
一応、下のように図が出ますが、なんか変ですね。
このプロットの実行時には下記のWarningが出ます。
Warning message:
In xy.coords(x, y, xlabel, ylabel, log) :
imaginary parts discarded in coercion
虚数部無視したよってことでしょうか。exp(ix)=cos(x+)isin(x)ってことで、余弦波となっているのですね。ふむ。
ggplot2の場合、同じようにすると
data.frame(x=seq(0,1,0.01),y=exp(-1i*seq(0,1,0.01)*2*pi)) %>% ggplot(aes(x=x,y=y))+geom_point()
これはエラーが出てプロットされません。エラーは下記のような内容。複素数のままだと大小関係わからないから確かにプロットできませんよね。エラーで止まるggplot2の方が親切かもしれません。
複素数は実部と虚部の2つをどうにかして表す必要があるので、最初にプロットしたように分けてプロットすることが考えられますが、もう一つ、絶対値と角度の形でも表現できます。極座標形式でしたっけ？それを使用したものがボード線図なのです。続く。
ようやくbode線図について。
ボード線図、あるいはボーデ線図と呼ばれる、周波数に対してのゲインと位相の関係をプロットした図については説明を省略…。
周波数応答系はわかりやすいというのは難しい気がするなぁ。一応参考になるかもなサイトを。
www.kairo-nyumon.com
さて、最初に貼ったサイトからソースを転載しますが、
fseries <- function(min, max, perdec) {
expmin = floor(log10(min))
expmax = floor(log10(max))
step = 1/perdec
return(10^(seq(expmin,expmax,by= step)))
}
Hf <- function(f, R, L ,C) {
s = as.complex(2i*pi*f)
Zser = as.complex( R + s*L)
Zpar = as.complex(1/(s*C))
return(Zpar/(Zpar+Zser))
}
ここまでは同じものを。先の関数は対数で見たときに等間隔になるようにする関数、後の方は伝達関数そのものを書いていますね（たぶん）
ここからはちょっと手を加えています。
ゲイン側は同じものを使用していますが、それに加えて位相分をArgによって取り出し、単位をradからdegへと変換しています。
library(ggplot2)
library(scales)
f = fseries(100,1e6,100)
R = 0.06; L = 1E-6; C = 47E-6
amplitude = 20*log10(Mod(Hf(f,R,L,C)))
df <- data.frame(x = f, y = amplitude)
p1 <- ggplot(df,aes(x,y)) + geom_line(lwd = 1, col = 'blue') +
scale_y_continuous(breaks=seq(-65,25,by = 5)) +
labs(x = "f [Hz]", y = "|H(f)| [dB]") +
scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+
annotation_logticks(sides = "bt")
phase = Arg(Hf(f,R,L,C))/pi*180
df <- data.frame(x = f, y = phase)
p2 <- ggplot(df,aes(x,y)) + geom_line(lwd = 1, col = 'blue') +
scale_y_continuous(breaks=seq(-180,180,by = 30)) +
labs(x = "f [Hz]", y = "Arg(H(f)) [deg]") +
scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x)))+
annotation_logticks(sides = "bt")
grid.arrange(p1,p2)
これで、（私）見慣れたbode線図の形となりました。
最後の、プロットを並べて表示についてはこちらを参考にさせていただきました。
id.fnshr.info
tidyverseからマスクされてますよってwarningが出ましたが、今回は使わないので良しとしました。
これ、geom化しようかと考えてはいたのですが、上記の通りある程度処理を加えればなんとかなってしまうので、わざわざgeomを作る必要はないのかな、と思いました。
ggplot2の中身の理解のために、これを題材に作るのは良いかもしれません。課題です。
終わりに
ということで、誰が得するのかわかりませんが、Rでも複素数を扱えて、bode線図も書くことが出来ますよ、という紹介でした。
これをRの得意な統計の分野とつなげられれば、使い道があるんですが…
しかし、リンク先にもあるように、
Well, this started as one of those “I'll do it just because I can" ...
です。Rは意外と色々なことができるので、色々試してみてください。
というわけで...
Enjoy!
