2020-11-08

熱雑音の導出

A：抵抗の熱雑音のパワースペクトル密度は $4kTR$ っていいますけど，導出についての論文を見ると2つの抵抗を伝送線路で接続している系を考えるみたいなんです．なんで伝送線路が急に登場するのかがわからなくて，そこから先が頭に入ってこないんです．

B：確かに唐突な印象がありますね．でも，伝送線路はエネルギーを考えるのに結構都合がいいんですよ．

A：というと？

B：理由は2つあります．１つ目は，伝送線路の両端を終端されているときに伝送線路に閉じ込められているエネルギーが，定在波の数から計算できることです．

A：確かに，論文では伝送線路に抵抗をつないだ状態から，急に抵抗を切り離して伝送線路を終端するようなことが書いてありました．

B：両方が固定端になるので，定在波のモードがわかりますね．

A：周波数間隔が一定になるやつですね．伝送線路をエネルギーが伝搬するのにかかる時間を $\tau$ とすると，周波数間隔 $\frac{1}{2\tau}$ で定在波のモードが存在することになる，であってます？

B：正解．それを伝送線路の自由度とみると，エネルギー等分配則（1自由度当たりのエネルギーが $\frac{kT}{2}$ ）から，伝送線路に閉じ込められるエネルギーは全部で $4\tau (f_1-f_0) \frac{kT}{2}$ と分かります．

A： $f_1$ と $f_0$ は，対象とする周波数の上端と下端ですね．で，定在波のモード数を自由度とみるわけだ．でも，2倍大きい気がするんですが．

B：一つの周波数で，磁界と電界の2自由度があるので2倍するみたいです．

A：みたい？

B：そこはよくわからないです．それはおいといて，もう一つの理由にいきましょう．２つ目の理由は，検討対象のエネルギーを伝送線路内に通せるってことです．

A：？

B：ある抵抗があって，その一方が接地されていて，他方が伝送線路を介して別の抵抗で終端されている系を考えてみてください．2つの抵抗と伝送線路の特性インピーダンスが一致していれば，熱雑音電圧の半分，つまり熱雑音電力の 1/4 が伝送線路を通って反対側の抵抗に向かっていくことになります．その状態で伝送線路を切り離すと，熱雑音電力が伝送線路に閉じ込められますよね？

A：なるほど，伝送線路に閉じ込めてしまいさえすれば，１つ目の話に近づきそうです．でも，こっちは熱雑音「電力」で，さっきはエネルギーだったので，次元が合いませんよ．時間積分が必要ですがどうするんですか？

B：伝送線路をエネルギーが伝搬する時間 $\tau$ を使います．定常状態では電力は時刻によらず一定なので，電力に $\tau$ をかければエネルギーにできます．式にしてみましょう．

A：えっと，熱雑音パワースペクトルを $Sn(f)$ とすると，これをさっきと同じ周波数範囲で積分するので，伝送線路を通る熱雑音のエネルギーは $\tau\frac{1}{4R}\int^{f_1}_{f_0}Sn(f)df$ ですか？

B：惜しい．抵抗が2つあって，双方から同じだけ無相関な熱雑音電力が伝送線路に入ってくるので，その2倍になります． $2\tau\frac{1}{4R}\int^{f_1}_{f_0}Sn(f)df$ ですね．

A：また2倍か・・・．

B：まぁまぁ．で，1つ目と2つ目が等しいはずだから， $4\tau (f_1-f_0) \frac{kT}{2} = 2\tau\frac{1}{4R}\int^{f_1}_{f_0}Sn(f)df$ となります．もうすぐゴールです．

A： $f_0 = 0$ でもいいわけだから， $4\tau f_1 \frac{kT}{2} = 2\frac{1}{4R}\int^{f_1}_{0}Sn(f)df$ ．さらに，両辺を $f_1$ で微分すると， $\tau kT = \frac{1}{4R}Sn(f)$ となって，整理すると．．． $Sn(f) = 4kTR$ ！！

B：めでたしめでたし．

A：伝送線路を持ち込んだ発想がすごいですね・・・．

B：本当にそうですね．実は，偉そうしてましたが California Institute of Technology の講義動画がとても分かりやすくて，その受け売りなんです．ありがたいですね．

2017-05-22

juliaで可変ワード長の演算を試す

juliaが結構便利なので、勉強がてら有限ワード長のモジュールを作っています。その過程をメモ代わりにそのまま吐き出してみようと思います。

概要を考える

ワード長は可変にしたいので、値とワード長の二つのフィールドが必要になりそうです。値の方は、簡単のため UInt で済ませます。演算をするのであれば、コンストラクタと加減算、あとは、連結はほしいところです。

というざっくりしたイメージでスタートしました。

コンストラクタ

こちらは、特に悩む必要もなくtypeを使うだけで終わります。任意ワード長の RegVarのほかに、8bit固定の Reg8 も用意しておこうかと思い、これらをまとめる抽象型 AbstractReg も作りました。

abstract AbstractReg
type RegVar <: AbstractReg
  val :: UInt
  wl :: Int
end
type Reg8 <: AbstractReg
  val :: Uint
end

加算

Baseの+関数をimportして、ワード長を考慮しながら計算します。

import Base: +
export +
+(a::RegVar, b::RegVar) = addreg(a, b)
function addreg%(a::RegVar, b::RegVar)
  n = max(a.wl, b.wl)
  ub = 2^(n-1)-1
  lb = -2^(n-1)
  y = a.val + b.val
  if (y < lb)
    y = y + 2^n
  else (y > ub)
    y = y - 2^n
  end
  RegVar(y, n)
end

RegVar としてはこれでいいのですが、Reg8に適用しようとしてもReg8はwlをメンバに持っていないのでうまくいきません。かといって、Reg8にwlを追加するのも嫌だな……。と考えたところで、wl を総称関数にすればよいことに気づきました。

wl(::Reg8) = 8
wl(x::RegVar) = x.wl

これをつかって先ほどのaddregを書き換えると、こうなります。

import Base: +
export +
+(a::AbstractReg, b::AbstractReg) = addreg(a, b)
function addreg(a::AbstractReg, b::AbstractReg)
  n = max(wl(a), wl(b))
  ub = 2^(n-1)-1
  lb = -2^(n-1)
  y = a.val + b.val
  if (y < lb)
    y = y + 2^n
  else (y > ub)
    y = y - 2^n
  end
  RegVar(y, n)
end

a.wl などとしていたところを wl(a) にしてあげただけですが、これで Reg8 にも RegVar にも対応できるようになりました。

Juliaの場合、型のフィールドに対して自動的にアクセサが生成されるのですが、総称関数を定義するとこういうときに便利なのですね。納得。

減算

減算は、単項演算子のマイナスを定義して、先ほどの加算と組み合わせることにします。

import Base: -
export -
-(x::AbstractReg) = RegVar((1 << wl(x)) - x.val, wl(x))
-(a::AbstractReg, b::AbstractReg) = a + (-b)

ワード長が64の時にどうするんだという話はありますが、とくに問題ないようです。

julia> convert(UInt, 1) << 64 - 1
0xffffffffffffffff

連接

このあたりから自信がなくなってきます。演算子は ++ にすることにして、二項演算子としてはすぐに書けました。ただ、これを連続して使おうとした場合にどうするかで手が止まります。よくわからないので、julia本体から似たようなコードを探してみることに。sum なら可変長引数をとるかとおもって確認してみたところ、tuple.jlに記載がありそうです。

julia> methods(sum)
# 13 methods for generic function "sum":
sum(x::Tuple{Any,Vararg{Any,N<:Any}}) at tuple.jl:205
...

そこを確認したところ、ドット３つでタプルを展開してくれるような感じです。いわれてみれば、MATLABにそんなのがあったような気もします。

export ++
++(args...) = assoc(args)
assoc(x::Tuple{Any, Vararg{Any}}) = assoc(x...)
assoc(a::RegVar, args...) = assoc(a, assoc(args))
...

こんな感じで展開できました。

補足：表示、等号

書き始めた当初に、完全に忘れていたのが表示と等号。これらがないと REPL やテストで困ることに、後で気づきました。表示については、juliaの本体を検索すると、showという関数を定義すればよさそうなことがわかりました。

> ag "print\(" *.jl
...
version.jl
48:function print(io::IO, v::VersionNumber)
49:    v == typemax(VersionNumber) && return print(io, "∞")
50:    print(io, v.major)
51:    print(io, '.')
52:    print(io, v.minor)
53:    print(io, '.')
54:    print(io, v.patch)
...

というわけで、定義します。

import Base: show
function show(io::IO, x::AbstractReg)
  mroundup(x, n) = Int(ceil(x / n)) * n;
  print(io, wl(x));
  print(io, "'h");
  octmask = 1 << 4 -1;
  for i=mroundup(wl(x), 4)-4:-4:0
    z = octmask & (x.val >> i);
    if z < 10
      print(io, z);
    else
      print(io, 'a'+ (z-10));
    end
  end
end

これだけ定義すれば、REPL や print 関数などから自動的に呼んでもらえます。

julia> Reg8(10)
8'h0a

等号はこんな感じ。

==(a::AbstractReg, b::AbstractReg) = (a.val == b.val) && (wl(a) == wl(b))

できた。

2016-12-07

配列の配列から行列への変換

配列の配列を、行列に変換する方法がわからなかったのでメモ。
結果としては以下のようにすれば良いようです。

import Base.convert
convert{T}(::Array{T,2}, x :: Array{Array{T,1},1}) = hcat(x...)'

今回学んだのは２点。

一つは、型自体で多重ディスパッチしたい時は ::型でいけるということ。
もう一つは、... の挙動が理解できていないこと。
後者は、現段階では「配列やタプルを分解してくれるもの」というぼんやりした認識です。そういえば MATLAB にも同じ文法がありますね。
詳細について、いずれドキュメントを確認しようと思います。

2016-11-28

Julia で scanl を使いたい

Julia で scanl に相当するものが見当たらなかったので、勉強がてら、考えてみました。

foldl が参考になるのではないかと思ってソースを見てみたところ、以下のようになっています。

foldl(op, v0, itr) = mapfoldl(identity, op, v0, itr)
mapfoldl(f, op, v0, itr) = mapfoldl_impl(f, op, v0, itr, start(itr))
function mapfoldl_impl(f, op, v0, itr, i)
  if done(itr, i)
    return r_promote(op, v0)
  else
    (x, i) = next(itr, i)
    v = op(r_promote(op, v0), f(x))
    while !done(itr, i)
      @inbounds (x, i) = next(itr, i)
      v = op(v, f(x))
    end
  end
end

色々呼び出していますが、f は map 用の関数で、op は v に結果を次々と付け加えていく関数と考えれば良さそう。
Haskell と（雰囲気だけですが）対応させると、こんな感じでしょうか。

itr :: (Ix i) => (b, i)
f :: b -> a
v :: a
op :: a -> a -> a  -- mappend 相当
mapfoldl f op v0 itr = foldl op v0 $ map f itr

scanl を作るためには、２つの引数を取る関数とmappend が必要ですが、 mapreduce には前者がないので、これで scanl を実現するのはちょっと難しそうです。

以上を踏まえ、素直に書くことにしました。こんな感じです。

 function scanl{T}(op, v0 :: T, xs)
   function scanList(op, vs, xs)
     isempty(xs) ? vs : scanList(op, [vs; op(vs[end], xs[1])], xs[2:end])
   end
   scanList(op, [v0], xs)
 end

再帰を使っていますが、foldl の実装を見ると素直にループを回した方が良いのかもしれません。この辺りも、Julia なら簡単に検討できるはずなので、いずれやってみようと思います。

2016-11-25

Julia で関数合成マクロをつくる

Julia は便利そうなのですが、Haskell から移ってくるとカッコの数が気になってきます。
せっかく Julia ではマクロが作れるので、関数合成を簡単に書けないかと考えてみました。一気に書くことができなかったので、関数列と引数を与えると合成した結果を返すマクロ @rc を用意し、さらに @compose で無名関数にするという２段構えです。
（もっとスマートな方法があれば教えてください）

macro rc(args ...)
  if length(args) == 1
    :($(args[1]))
  else
    :($(args[1])(@rc($(args[2:end]...))))
  end
end

macro compose(args ...)
  :(x -> @rc($(args ...), x))
end

rc は、recursive compose のつもり。使い方はこんな感じになります。

>>> @rc(asin, sin, 1)  # => 1
>>> @compose(asin, sin)  # =>(generic function with 1 method)
>>> @compose(asin, sin)(1)  # => 1

2016-11-01

MATLABでカリー化

一度関数型プログラムの考え方に触れると、どの言語でも同じようなことがしたくなってしまいます。
MATLAB は便利ですが、当然関数型ではないので勉強がてらちょっとずつ自作しています。
その過程で、カリー化をしようとしてハマったのでメモ。

結論からいうと、セル配列が使えました。セル配列は型の違うものをひとまとめに扱うしくみで、
取り出すときには小括弧と中括弧の２つの方法があります。
中括弧を使って範囲指定をして取り出すと、多値で答えが帰ってくるので、これを使って curry がつくれました。

function f = curry(f, varargin)
%curry do currying
%   feval(curry(@plus, 1), 10) => 11
%   feval(curry(@strcat, 'a', 'b', 'c'), 'z') => 'abcz'

if nargin >= 2
    xs = varargin;
    f = curry(@(varargin) f(xs{1}, varargin{:}), xs{2:end});
end

第１引数の f をそのまま返値にしているので、 nargin == 1 の場合はそのまま f が返ります。
この辺り、お行儀が良いのかはよくわかりません。明示的に分岐を書いた方が読みやすい気もします。

# 追記
と思ったら、まんま stack overflow に書いてあった。。。

2016-08-17

Stackの味見

他人事だと思っていた cabal hell についにはまったので、cabal から stack に乗り換えて見ることにしました。

brew uninstall ghc
brew install ghc
brew haskell-stack
stack new test-proj
cd ./test-proj
vim ./test-proj.cabal
# library の build-depends にカンマ区切りで使用するパッケージを追加
stack build
stack ghci

ソースファイルに import を書くのに加え、.cabal ファイルにも記載が必要なので二度手間のような感じです。
多分、ソースから .cabal に反映させる方法はあるのでしょうが、まだ見つけられていません。