π型フィルタの計算

ラインインピーダンスＺとカットオフ周波数fを決めたら
以下の式にしたがってＬ、Ｃを計算する。
Ｚ＝２πfＬ＝１／２πfＣ

エラーアンプ付きスイッチングレギュレータ メモ

○出力電圧の計算
エラーアンプのー端子（ＦＢ端子）に戻す電圧で決まる。
具体的には出力電圧を抵抗分圧した電圧を戻してやる。で、エラーアンプはスイッチングレギュレータの内部で作られた参照電圧とＦＢ端子の電圧を比較して次段にあるコンパレータの参照電圧を作る。コンパレータのもう一方の端子には三角波が入力されていて、エラーアンプの出力電圧が変わればデューティも変わって出力電圧が変化する。

フィードバックループ概略
出力電圧が下がる→ＦＢ端子（ー端子）電圧下がる→エラーアンプ出力上がる→コンパレータの参照電圧上がる→デューティ比が増加　or 減少（外付けのパワトラがPNPかNPNかＰ型かＮ形かで増加するか減少するかは変わる）→出力電圧上がる

外部メモリのアドレス線接続についてメモ

バス幅が16bitのメモリを接続する場合は、
以下のようにマイコンのアドレス線を1bitシフトダウンする。

＜16bitメモリとマイコンのアドレス線接続＞

マイコンのアドレス端子　　　　外部メモリのアドレス端子
A1-------------------------------A0
A2-------------------------------A1
A3-------------------------------A2


バス幅が32bitメモリを接続するんなら2bitシフトダウンして以下のようになる

＜16bitメモリとマイコンのアドレス線接続＞

マイコンのアドレス端子　　　　外部メモリのアドレス端子
A2-------------------------------A0
A3-------------------------------A1
A4-------------------------------A2

なんでシフトするかといえば、たぶんマイコンはバイトアクセスするつもりでアドレッシングするからだと思う。（＝バス幅が8bitのメモリがくっついているつもりでアドレッシングする）
だから、アドレス線が16本あるマイコンが０バイト目をアクセスしたいときには
A[0:16]=0x0000
１バイト目だったら
A[0:16]=0x0001
２バイト目だったら
A[0:16]=0x0002
とアドレス出力する。


○16bitバスのメモリがくっついている場合
マイコンが０バイト目のデータが欲しくて
A[0:16]=0x0000
というアドレスを出力すると
メモリは０バイト目から16bit分をデータバスに吐き出す。
つまり、０バイト目と１バイト目を吐き出す。
この場合は余分な１バイト目がくっついているけど、マイコンが欲しい０バイト目も出てきているから問題ない。余分な分は無視すればいいのだから。
では、マイコンが１バイト目のデータが欲しいときどうなるか？
A[0:16]=0x0001
というアドレスをメモリがもらうと、1 x 16bit目から16bit分を吐き出す。
つまり、2byte目と3byte目がでてきてしまう。
これにはマイコンが欲しかった1バイト目は含まれていない。
こんな問題が起きてしまうので、1bitシフトダウンしておく。
こうすれば、ぴったりマイコンが欲しいと思ってアドレッシングしたデータが出てくる。

16bitメモリだから1bitシフトだが、32bitならもう1bitシフトしないといけないので結局2bitシフトすることになる。

Pyscoで高速化

Pyscoってのを使うと早くなるといううわさがあるのでメモ
import psyco
psyco.full()
と冒頭に入れるらしい。

Super bound unboundについて

superという関数について調べる。


class super(object)
| super(type) -> unbound super object
| super(type, obj) -> bound super object; requires isinstance(obj, type)
| super(type, type2) -> bound super object; requires issubclass(type2, type)
| Typical use to call a cooperative superclass method:
| class C(B):
| def meth(self, arg):
| super(C, self).meth(arg)

super(type)っての要するにtypeの基底クラスを与えるのか？
では、
super(type,obj) ->bound super object.ってどういうこと？


しらべてみるとboundってのはメソッドがどっかのインスタンスに属している。unboundは属していないってことらしい。

boundメソッドは、そのまま関数として使える。unboundメソッドはそのまま関数として使おうとするとだめ。

次のコードを例にとって考えてみる。

class Bassclass:

    def __init__(self,str1):

        self.a = str1+"bass class"

    def func1(self):

        return self.a



print Bassclass().func1()

C++だと、コンストラクタにselfなんて渡さなかった（たぶん）けど、pythonだとselfを渡さないとだめ。また、何も引数がいらなさそうなfunc1()もselfを渡すような形で定義してやら無いとダメ。
上のコードだとBassclassっていうスコープの中だからいちいちself.aなんて書かなくてもいいのかと思ってたけどそうじゃないみたい。


疑問点が１つ。
class Bassclass(object)
という文があるが、これを
class Bassclass()
にするとなんか知らんけどだめ。

新スタイルクラスとクラシックスタイルクラスってのがあるらしい。
新スタイルクラスではobjectの派生クラスとして定義しないとダメらしい。
んで、新スタイルクラスではスタティックメソッドとクラスメソッドが使えるそうな。
スタティックメソッドは、インスタンスなしで使用できるメソッド。普通のクラスに属さないメソッドみたいなもんか。
クラスメソッドは、クラスを渡すメソッド。こっちは全く分からん。

class Bassclass(object):

    def __init__(self,str1):

        self.a = str1+"bass class"



    def func1(self):

        return self.a



class Derivedclass(Bassclass):

    def __init__(self, str1, str2):

        super(self.__class__, self).__init__(str1)

        self.b = str2

    def func1(self):

        return self.a

    def func2(self):

        return self.b



d = Derivedclass("ouou",12)

print d.func1()

print d.func2()

Gnuplot-py のサンプル(マンデルブロー集合) Gnuplot-py sample code (for Mandelbrot set)

Gnuplot-pyライブラリをインストールすると、どこかにdemo.pyができているのでこいつを参考にしてコードを書いた。

import psyco

psyco.full()

import Gnuplot, Gnuplot.funcutils

from scipy import *



STOP_VALUE = 4



#CX_AREA=[-1.495,-1.475]

#CY_AREA=[-0.02,0.045]

CX_AREA=[0,0.4]

CY_AREA=[-0.65,-0.35]

SIZE = 300

YSTEP = (CY_AREA[1]-CY_AREA[0])/SIZE



XSTEP = (CX_AREA[1]-CX_AREA[0])/SIZE



def mandelbrot_seq(x,y):

    mandelbrot_seq = 0.+0.j

    c=x+y*1j

    for n in range(0,100):

        mandelbrot_seq=mandelbrot_seq**2+c

        if abs(mandelbrot_seq)>STOP_VALUE:#avoid overflow(mandelbrot_seq is too large number)

            break

        else:

            pass

    return abs(mandelbrot_seq**2+c-mandelbrot_seq)





g = Gnuplot.Gnuplot(debug=1)

#g('set parametric')

#g('set data style lines')

g('set hidden')

g('set pm3d map')

#g('set contour base')

#g.title('An example of a surface plot')

g.xlabel('Re axis')

g.ylabel('Im axis')



x = arange(CX_AREA[0],CX_AREA[1],XSTEP)

y = arange(CY_AREA[0],CY_AREA[1],YSTEP)



g.splot(Gnuplot.funcutils.compute_GridData(x,y, mandelbrot_seq, binary=0))

g('set pm3d')をg('set pm3d map')に変更すると底面のマップのみが表示されるのでこっちのほうが見やすい。

マンデルブロー集合の描画 Rendering Mandelbrot set

pythonの練習でマンデルブロー集合を描画するコードを書いた。

Fig:実行結果　The result of execution

from scipy import *

from pylab import *

STOP_VALUE = 4

CX_AREA=[-1,0.5]

CY_AREA=[-1,1]

STEP = 0.01

def judge_conv_mandelbrot_seq(c):

    mandelbrot_seq = 0.+0.j

    for n in range(0,50):

        mandelbrot_seq=mandelbrot_seq**2+c

        if (abs(real(mandelbrot_seq))>STOP_VALUE)or(abs(imag(mandelbrot_seq))>STOP_VALUE):

            break

#    if abs(mandelbrot_seq)    if abs(mandelbrot_seq)>STOP_VALUE:

        return 0

    else:

        return 1

    

mandelbrot_set=[]

for cx in arange(CX_AREA[0],CX_AREA[1],STEP):

    for cy in arange(CY_AREA[0],CY_AREA[1],STEP):

        if judge_conv_mandelbrot_seq(complex(cx,cy)):

             mandelbrot_set.append(complex(cx,cy))

        else:

            pass

plot(real(mandelbrot_set),imag(mandelbrot_set),',')



show()

Gnuplot-pyを使った3次元グラフの書き方

3D plot をしたい。
matplotlibにplot3Dという関数があってこいつで3Dグラフをレンダリングできるらしいのだが、バージョン違いで上手く動かない。
そこでGnuplotを使う。

Sample Code

import Gnuplot

import numpy



gp = Gnuplot.Gnuplot(debug=1)

x=[0,1,0]

y=[0,1,0]

y2=[0,1,0]



d3=Gnuplot.Data(x,y,y2)

gp.splot(d3)

こんな感じでいける。
matplotlibと違って、いったんGnuplot.Data()でなんか知らないけど型を変換しないといけないらしい。
まだ良く分かっていないのでとりあえずメモ。

PythonでRC回路の過渡解析をしてみる その2 数値計算 by Scipy編

前回建てた微分方程式と初期条件
Vc'=1/(RC)(Vin-Vc)

Vin(0) = 0.0
Vc(0)= Vcc

前回、手計算で解いた結果
Vc(t)=Vcc×exp{-1/(RC)×t}

でも、定数変化法で解くので手間が掛かる。


手計算ではなく数値計算する方法。

Pythonの数値計算モジュールScipyを使う。図１に実行結果を示す。

図１：実行結果

from scipy import *

from pylab import *

deriv = lambda y,t : array([y[1],-y[0]])

# Initial parameters

R=10**3

C=10**-7

Vcc=5.





def Vin(t):

    return 0.



def Vcdash(Vc,t):

    return (Vin(t)-Vc[0])/(R*C)



start=0

end=10**-3

numsteps=10000

time=linspace(start,end,numsteps)



from scipy import integrate

Vc0=Vcc

Vc=integrate.odeint(Vcdash,Vc0,time)



#Plot Setting

figure(1,figsize=(6,4)).subplots_adjust(bottom=0.11)

xlabel('t')

grid(True)

hold(True)



plot(time,Vc)

#Setting Window title and Save Figure

title('Transient Analysis of RC Circuit')

savefig('transient analysis.png',dpi=72)



show()

integrate.odeint()関数には、微分方程式と初期値と時間に相当する行列を渡すと計算してくれます。
微分方程式はこの場合Vcdash()です。
なんで微分方程式が
def Vcdash(Vc,t):
return (Vin(t)-Vc[0])/(R*C)
って表現されるのか？その理由を以下に示します。

＜リストを用いて表現される理由＞
微分方程式を数値解析的に解くには差分方程式だったか漸化式だったかを用いる。（たしか）
差分方程式@wikipedia

解きたい微分方程式

Vc'=1/(RC)(Vin-Vc)

を前進差分で差分方程式にしてやると、
Vc'=(Vc(t+h)-Vc(t))/h=1/(RC)(Vin-Vc(t)) --(1)
ここで、hを数値計算するメッシュの分解能として考えると（＝メッシュがhずつ刻まれていると考えると）、
Vc(t) =V[n]
Vc(t+h)=v[n+1]
とおける。
すると(1)式は

(Vc[t+1]+V[t])/h=(Vin-Vc[t])/RC
V[t+1]=Vin/RC-Vc[t](h-1/RC)

となり、ここでhを0に近づければ、（計算機では０にはできないが）
Vc[t+1]≒1/(RC)(Vin-Vc[t]) --(2)


(1),(2)を見比べると、Vc'はVc[t+1]に、VcはVc[t]に置き換わっただけということが分かる。ということで、微分方程式はリストを使って表現できる。
2階微分方程式の場合、例えば、

y''=ay'+by+c

を解きたい場合には、

y'=x

とおくと

x'=ax+by+c

のように、1階常微分方程式になりscipy.integrate.odeint()関数で扱える形式になる。
ただし、今度は、yの初期値だけでなくy'=xの初期値も必要となってくるので注意が必要。
参考Scipy CookBook




ついでに、このプログラムのVin(t)という関数は電源電圧を表現している、これをsin(ωt)を返す関数にした例↓。



図２:Vin(t)=Vcc*sin(10**5*t)とした場合の実行結果

PythonでRC回路の過渡解析をしてみる その１ 手計算

↓の回路を対象に計算。最終的にはVcの式が知りたい。



まずは微分方程式を建てる。

Vc=1/C∫Idt --(1)
I=Vr/R --(2)
Vin-Vc=Vr --(3)
この３つの式をまとめる。
Vcをみたいので、VcとVinの関係式に持ち込むという方針で変形する。
(3)を(2)に代入すると、
I=(Vin-Vc)/R --(3')
(3')を(1)に代入すると
Vc=1/C∫(Vin+Vc)/Rdt
Rは定数なので移動して
Vc=1/(RC)∫(Vin-Vc)dt --(微分方程式)


<初期条件を建てる>
次に初期値を設定。
設定としては、Vinの電圧がVcc→0.0に落ちたということにします。あとは、コンデンサはフルチャージされていることとする。

Vin(-0) = Vcc
Vin(+0) = 0.0
Vc(0)= Vcc（コンデンサには電荷が完全チャージされている状態を想定。）

上で導出した微分方程式を並べる。

Vc=1/(RC)∫(Vin-Vc)dt --(微分方程式)
話が前後しますが、微分方程式を微分形式に直します。両辺を微分すると、

Vc'=1/(RC)(Vin-Vc) --(4)

微分したら積分定数が消える気がするが、積分定数とはすなわち初期条件のことなので（たぶん）ＯＫ。

<手計算で解く>
この微分方程式は、非斉次方程式ですので、定数変化法で解きます。

まずは、斉次方程式の一般解を求める。
Vc'=1/(RC)(Vin-Vc)
を斉次化するには
非斉次項Vinを
Vin=0とおけばよいので
Vc'=-1/(RC)Vc --(7)

次に、(5)に初期条件を代入してやると、
Vc(0)=Vcc=A×exp(0)
∴A=Vcc
よって
斉次方程式の一般解＝Vcc×exp{-1/(RC)×t} --(8)

(8)より
非斉次方程式の一般解＝B(t)×Vcc×exp{-1/(RC)×t}
となっているはずなので、B(t)を求める
初期条件を代入すると、
B(t)=1
となる。よって、

Vc(t)=Vcc×exp{-1/(RC)×t}
となる。
この式は、t=0でVc=Vccとなり、t>0の領域で指数関数的に減少するカーブを描く。
これはコンデンサの放電時の波形。

pythonワンパッケージ開発環境（@windows）

pythonの開発環境を構築するのは結構めんどくさい。

プレーンなPython開発環境を構築するのはあっさりできるけど、モジュール（ライブラリ）を追加していくのがメンドクセー。
ScipyとかNumpyとかMatplotlibとかVpythonとか便利そうなモジュールを使おうとするとだいたいわけ分からんことになって解決するのに時間がかかる。
私は、そんなことを解決するのには興味が無いのでこういう作業は苦痛でしかない。

私のような人にはPython(x,y)とかEnthought社が出しているワンパッケージの開発環境がお勧めだ。
これらは、pythonインタプリタ、数値計算モジュール(scipy,numpy,matplotlib)がワンパッケージになっている。ついでにPython(x,y)はEclipse+pyDevがついてくる。

めんどくさがりにはうってつけだ。

pywinautoインストール

python2.6にpywinautoをインストールした。

pywinautoはwindowsのGUIアプリを自動操作するのに便利そうなモジュール。

依存関係にあるモジュールをインストールして、いよいよpywinautoのモジュールをインストールしますかねというところでつまずいた。

pywinautoのsetup.pyを使ってコマンドプロンプト上で

setup.py install

って打ってみたら、ctypesモジュールないとのエラー。

よく分からないがpyScripterにsetup.pyを突っ込んで引数にinstallを指定して動かしたらインストール成功。

python help()関数

例えば、

print help(scipy.help.integrate.quad)

としてやれば、

Help on function quad in module scipy.integrate.quadpack:

quad(func, a, b, args=(), full_output=0, epsabs=1.4899999999999999e-008, epsrel=1.4899999999999999e-008, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)
Compute a definite integral.

Description:

Integrate func from a to b (possibly infinite interval) using a technique
from the Fortran library QUADPACK. Run scipy.integrate.quad_explain()
for more information on the more esoteric inputs and outputs.

Inputs:

func -- a Python function or method to integrate.
a -- lower limit of integration (use -scipy.integrate.Inf for -infinity).
b -- upper limit of integration (use scipy.integrate.Inf for +infinity).
args -- extra arguments to pass to func.
full_output -- non-zero to return a dictionary of integration information.
If non-zero, warning messages are also suppressed and the
message is appended to the output tuple.

Outputs: (y, abserr, {infodict, message, explain})

y -- the integral of func from a to b.
abserr -- an estimate of the absolute error in the result.

infodict -- a dictionary containing additional information.
Run scipy.integrate.quad_explain() for more information.
message -- a convergence message.
explain -- appended only with 'cos' or 'sin' weighting and infinite
integration limits, it contains an explanation of the codes in
infodict['ierlst']

Additional Inputs:

epsabs -- absolute error tolerance.
epsrel -- relative error tolerance.
limit -- an upper bound on the number of subintervals used in the adaptive
algorithm.
points -- a sequence of break points in the bounded integration interval
where local difficulties of the integrand may occur (e.g.,
singularities, discontinuities). The sequence does not have
to be sorted.

**
** Run scipy.integrate.quad_explain() for more information
** on the following inputs
**
weight -- string indicating weighting function.
wvar -- variables for use with weighting functions.
limlst -- Upper bound on the number of cylces (>=3) for use with a sinusoidal
weighting and an infinite end-point.
wopts -- Optional input for reusing Chebyshev moments.
maxp1 -- An upper bound on the number of Chebyshev moments.

See also:
dblquad, tplquad - double and triple integrals
fixed_quad - fixed-order Gaussian quadrature
quadrature - adaptive Gaussian quadrature
odeint, ode - ODE integrators
simps, trapz, romb - integrators for sampled data
scipy.special - for coefficients and roots of orthogonal polynomials

None

とかって表示される。

python,scipyメモ

小数を含むリストを自動生成したい。
pythonの組み込み関数range()だと

range(0.0 , 10.0 , 0.1)

みたいなことをすると、

step arguments must not be zero

エラーになる。step argumentsをzeroにしているつもりはないがダメらしい。内部的にintにキャストされる？

そこでscipyモジュールのarange()関数を使うとうまくいった。

from scipy import *

from scipy import integrate

from pylab import *



start=0

end=2*pi

step=0.001

value = []

for x in arange(start,end,step):

    value.append(integrate.quad(sin,start,x))



time=arange(start,end,step)

plot(time,value)

show()

みたいなことをやってもＯＫでした。不思議。