謹賀新年
（遅くなりましたが）あけましておめでとうございます
本年もよろしくお願いいたします。
今年の干支をウェーブレットで描かしてみました。
必要な環境は id:ikeuchihiroki:20090102 と同じです。
from pylab import * import iwavelets.pycwt as w #Generate Time series Fs = 1024.0*4 time = arange(0,3.0/4.0+1.0/8.0,1/Fs) #horizontal lines signal=sin(2*pi*200*time)+sin(2*pi*100*time) signal[0:Fs/2+Fs/8]+=sin(2*pi*300*time[0:Fs/2+Fs/8]) win=zeros(signal.size) win[Fs/8:Fs-Fs/4]=(ones(Fs-Fs/4-Fs/8)) signal=signal*win #vertical lines t0=time[0:Fs/8] sig1=sin(2*pi*(100+t0*200*4)*t0) signal[Fs/4:Fs/4+Fs/8]+=sig1*2 signal[Fs/2:Fs/2+Fs/8]+=sig1*2 #Frequency array freqs=arange(1.0,400.0,1.0) #Morlet fc fc=5 #calc r=w.cwt_f(signal,freqs,Fs,w.Morlet(fc)) rr=r.real**2+r.imag**2 #show clf() imshow(flipud(rr),aspect='auto',extent=(time[0],time[-1],freqs[0],freqs[-1])) xlabel('Time[s]') ylabel('Frequency[Hz]') show() #output Wave file import wave f=wave.open('test.wav','wb') f.setnchannels(1) f.setframerate(Fs) f.setsampwidth(2) f.writeframes(signal) f.close()
Pythonで連続ウェーブレット解析
swanという連続ウェーブレット解析のモジュールがあった。
swan-0.5.4a.tar.gzを試してみる。
必要なモジュール
Gtkは最新版2.12.9を入れた。
PyGtkはバイナリからインストール。
swanはアーカイブを展開して、
python setup.py install
とする。ただし、Windowsの場合はそのままではインストールできないので、
setup.cfgを以下のように書き換えてから実行する。
# very basic config [install] prefix=C:\Python25 optimize=2
インストール後にC:\Python25\Scriptsに「swan」というファイルができるので
move swan swan.py
としておく。
試す
2. Simple usageを見ながら
まずは信号をつくる。「ipython -pylab」で起動して
import math,numpy,scipy,matplotlib Fs = 65.0 time = arange(0,16,1/Fs) signal = sin(2*pi*1.2*time) + cos(2*pi*3.7*time) save('signal.dat', signal)
実行はコマンドプロンプトから
swan.py -s 65 signal.dat
として、Apply CWTを押すと実行される。
pycwtのみを使ったプログラム
swan: some visible examplesにあるように、GUIなしでプログラムを書いてもよい。
import pylab as p import iwavelets.pycwt as w Fs = 65.0 time = p.arange(0,16,1/Fs) signal = p.sin(2*p.pi*1.2*time) + p.cos(2*p.pi*3.7*time) freqs=p.arange(0.1,5,0.025) fc=1.5 r=w.cwt_f(signal,freqs,Fs,w.Morlet(fc)) rr=r.real**2+r.imag**2 p.imshow(p.flipud(rr),aspect='auto',extent=(time[0],time[-1],freqs[0],freqs[-1])) p.xlabel('Time[s]') p.ylabel('Frequency[Hz]') p.show()
つかえそうだ。
numpy、SciPy、matplotlib
Pythonは科学計算にもよくつかわれるようだ。
科学者に必要なpythonモジュールはなにか
科学者のための Python 入門
とりあえず、Python本体に加えて以下のモジュールをインストールしてみる。
- Python
- SetupTools モジュールを簡単インストールできる
- numpy 行列などを扱える
- SciPy 数値計算など
- Matplotlib: Python plotting — Matplotlib 3.0.3 documentation グラフプロットなど。easy_installでインストール
- PIL 画像処理ライブラリ。easy_installでインストール
- http://ipython.scipy.org/moin/:pyreadline IPythonで使用するモジュール
- IPython インタラクティブシェル。easy_installでインストール
準備
Pythonはインストーラーで簡単にインストールできる。最新はPython 3.0だが使いたいモジュールが対応していないので、Python 2.5.2を選ぶ。
次に、SetupToolsをインストールする。同じくインストーラーでインストールする。SetupToolsは、easy_installというモジュールのダウンロードからインストールまでを自動でするコマンドが含まれる。
SetupToolsインストール後に、PATHを通しておく必要がある。マイコンピューター右クリック→プロパティ→詳細設定→環境変数→PATHに「C:\Python25\Scripts;C:\Python25;」を追加する。
モジュールをインストール
- numpy
- SciPy
- pyreadline
は、EXEのインストーラーでインストールする。easy_installの場合コンパイルが必要になるからだ。
pyreadlineはIPythonのサイトからダウンロードできる。
- matplotlib
- Python Image Library
- IPython
は、easy_installでバイナリをダウンロード、インストールできた。
コマンドシェルにて
easy_install matplotlib easy_install PIL easy_install IPython
と入力する。
試してみる
コマンドシェルにて
ipython -pylab
import math,numpy,matplotlib,scipy t=linspace(0,1,1024); y=sin(2*pi*10*t); plot(t,y) title 'Sin Plot' xlabel 'Time[s]' ylabel 'Amplitude'
けっこう使えそうだ。
Scilab（Matlabも）は、言語仕様が独特すぎてプログラムが読みにくい書きにくい印象があったが、Pythonならば可読性がよいプログラムが書けそうだ。
追記
http://www.enthought.com/に、科学計算に使うモジュール一式がある。必要に応じてサポート付きのライセンスもあるようだ。プロユースにも使えそう。
Scilabで連続ウェーブレット解析
id:ikeuchihiroki:20070222 で書いたように、連続ウェーブレット解析ができるフリーの環境というのはなかなか良いものがない。
Scilab上で動作するTOOLBOXでは以下のものがある。
- Wavelab
- Fraclab
- Scilab Wavelet Toolbox
Scilab Wavelet Toolbox http://scwt.sourceforge.net/ は、小規模なウェーブレットツールボックスなのだが、地味に更新が続けられている。
連続ウェーブレット変換のコードはCで書かれていて、なかなか速い。
テストプログラム
fs=512; x=[0:2047]/fs; y=sin(2*%pi.*x.*(x/max(x)*64)); // increase memory size stacksize(10000000); // continous wavelet transform scal=2^[1:1/32:6]; coef=cwt(y,scal,'cmor'); // colormap obtaining cmap=jetcolormap(512); // matrix normalization n=wcodemat(abs(coef),512); scf(0); clf; f=gcf();f.color_map=cmap; xlabel Time[s] ylabel Frequency[Hz] getf scal2freq.sci a=gca(); Matplot(n);
なかなかいい感じ。
でも、このままだと縦軸の意味がわからないので、スケールを周波数に変えて表示するようにしたい。
MATLABには、スケールを周波数に変換する関数が存在したので、同じような機能をする関数を書いてみる。
function [frq]=scal2freq(A,wname,fs) y=zeros(512,1); y($/2)=1; coef=cwt(y,A,wname); frq=[]; for k=1:length(A) fftcoef=abs(fft(coef(k,:))); [m ind]=max(fftcoef,'c'); tmp=ind-1; if tmp > length(y)/2 tmp=length(y)-tmp; end frq=[frq;tmp/length(y)*fs]; end endfunction
Aがスケール、wnameがウェーブレットの名前、fsはサンプリング周波数。
理屈は簡単で、インパルスをウェーブレット変換して、各スケールの中心周波数を見ているだけだ。たぶん理屈はあっているはずだ。CWTの関数の制約上、スケールは長さ２以上のベクトルを入れないといけないので注意。
次に、縦軸横軸を調整する関数。
function coefs=cwt_plot(data,scal,wname,fs) //plot settings y_ticks_distance=16; // increase memory size stacksize(10000000); coef=cwt(data,scal,'cmor'); // colormap obtaining cmap=jetcolormap(512); // matrix normalization n=wcodemat(abs(coef),512); scf; clf; f=gcf();f.color_map=cmap; xlabel Time[s] ylabel Frequency[Hz] getf scal2freq.sci a=gca(); ay=a.y_ticks; ay.locations=[0:y_ticks_distance:length(scal)]'; frq=scal2freq(scal($:-y_ticks_distance:1),'cmor',fs); ay.labels=string(frq); a.y_ticks=ay; Matplot(n); ax=a.x_ticks; ax.labels=string(ax.locations/fs); a.x_ticks=ax; endfunction
dataは解析対象信号,scalはスケール,wnameはウェーブレット関数の名前,fsはサンプリング周波数。
ウェーブレット関数の名前は、cwt.cを見るとわかる。今のところ以下の関数が使える。
sinus,poisson,mexh,morl,DOG, gaus1,gaus2,gaus3,gaus4,gaus5,gaus6,gaus7,gaus8, cmor,shan,fbsp,cauchy, cgau1,cgau2,cgau3,cgau4,cgau5,cgau6,cgau7,cgau8
使用例：
fs=512; x=[0:2047]/fs; y=sin(2*%pi.*x.*(x/max(x)*64)); cwt_plot(y,2^[1:1/32:8],'cmor',fs);
縦軸、横軸がサンプリング周波数に合わせて適正な表示になった。
意外にScilab Wavelet Toolboxが良い感じだ。
[Scilab] Scilabの変数スコープについて
Scilabは、変数のスコープの扱いが独特だ。
関数を作るときには、注意しないとややこしいことになるので、知っておかないとまずいと思い、ドキュメントを訳してみた。
http://wiki.scilab.org/howto/global_and_local_variables
より
In what follows level_0 is the level of the interactive shell; level_1 is that of functions called from level_0, level_2 that of subfunctions of level_1, and so on.
local variables and scoping
- local variables, which are not defined at level_N, are inherited from upper levels
- local variables defined at level_N remain local at that level. If a level_N variable has the same name of an upper level variable, a transient duplicate variable is created at level_N
- local variables inherited from upper levels can be read with no overhead; but if they are written, a transient local with the same name is created (not returned back unless specifically required). Note that the transient variable is first initialized to []; thus if this variable is implicitly grown (by assigning array elements beyond the current size), the local transient overrides over the outer variable.
- the appearance of the variable on an empty line of a function (which would cause immediate display at level_0) is sufficient to create the local duplicate
- input values of the function are assigned to the corresponding local variables of the function
- output arguments of a function are assigned upon return, at the upper level, to those variables specified in the function call
- argument passing is thus by value, in common programming terms; however, as in 3., the value is not really copied unless needed; for reading only, passing is by reference. So to say, it is a "lazy" passing by value, as copying the value is delayed to when really needed in the progam flow.
- only for adequately written primitives and interfaces of primitives, there is an option of by value/by reference (see intppty)
- local variables can be returned to their immediate upper level with resume. The difference between resume and output arguments is that resume specifies the name of the upper variable which is assigned, output values are assigned to whatever name choosen at the caller level.
- the command for undefining local variables is clear
Global variables
- global variables can have the same name of local variables, and coexist (the two namespaces are separate)
- global variables are shared by all scopes in which they are declared global, respectless of their calling depth
- declaring a variable global in a function does not help in making an upper local variable with the same name visible in the function: the upper local is already visible, and a global declaration can have unexpected side effects
- only when a variable is declared global, assignments affect its global version
- if a variable is declared global, and a local variable with the same name already exists, the value of the local is copied to the newly created global
- if a variable is declared global, and a local variable with the same name does not exist, an empty matrix [] is created in the global namespace and a pointer to it is created in the local namespace.
- the command for undefining global variables is clearglobal
適当に訳してみた。よくわからないところは意訳になっているので、内容は保障しません。
level_0はインタラクティブシェルのレベルである。
level_1はlevel_0から呼び出された関数であり、level_2はlevel_1から呼び出されたサブ関数である。それ以降も同様。ローカル変数とスコープ
- level_Nで定義されていないローカル変数は、上位のレベルから引き継がれる
- level_Nで定義されているローカル変数は、そのレベルに留まる。もし、level_Nの変数が上位レベルの変数と同じ名前である場合、一次的な複製がlevel_Nで作成される
- 上位レベルから引き継がれた変数は、オーバーヘッド無しに読み込むことができる。しかし、書き込みがある場合、一時的に同じ名前を持つ変数が作成される（特に指定されることがない限り上位レベルの値は更新されない）。この際に作成される一時的な変数の初期値は[]である点に注意すべきである。もし、この変数の領域が拡大する場合（配列の要素サイズが再定義されるような場合）、上位レベルの値は引き継がれない。
- 変数単体を記述したライン（level_0に内容が表示される）は、ローカルに複製を生成する
- 関数の引数は、ローカル変数になる
- 関数の戻り値は、関数を呼び出した上位レベルの変数になる
- 引数は値渡しでおこなわれるが、「3.」のような読み込みのみで複製の必要がない値は、参照渡しで行われる。値の複製が行われる「のろのろ」した値渡しは、プログラムを遅くする原因であるからである。
- 基本的には、値渡し・参照私のオプションは存在しない（intpptyを見よ）
- resumeを使うとローカル変数は上位レベルに戻すことができる。resumeと戻り値の違いは、resumeは上位レベルでの変数名を指定できるのに対し、戻り値は呼び出し側で名前を指定する
- ローカル変数を消去するコマンドはclearである
- グローバル変数は、ローカル変数と同じ名前を持つことができ、同時に存在できる（２つの名前空間は分けられている）
- グローバル変数は、コールの深さによらず、すべてのスコープで共有される
- 関数内でグローバル変数を定義しても、上位レベルで同じ名前のローカル変数は生成されない。上位レベルに同名のローカル変数が存在する場合に予期せぬ影響が考えられるからである。
- 値がグローバルのみで定義される場合、値の割り当てはグローバルにのみ影響をあたえる
- グローバル変数を定義する時、同じ名前のローカル変数がすでにある場合には、ローカル変数の値がグローバル変数にコピーされる
- グローバル変数を定義する時、同じ名前のローカル変数がない場合には、空行列[]がグローバル名前空間に生成され、そこへのポインターがローカル名前空間に生成される
- グローバル変数を消去するコマンドはclearglobalである
え？そうだったの？ってところがいくつもある。
- ローカル変数はサブ関数内にもスコープを持つ
- ただし読み込み専用
- 書き込むと同じ名前のローカル変数が作られる
- その際空行列[]で初期化されるため、行列の一部を書き換えるような書き方はできない
- 書き換えたとしても、上位レベルの値は変わらない
これを知らないと、変に悩むことがあるだろう。注意が必要だ。
- 追記
後で気がついたが、これらのルールはFortran風のルールのようだ。Fortranは使ったことがないので気がつかなかった。MATLABも同様なのだろうか。
遺伝的アルゴリズム in Scilab
遺伝的アルゴリズム(Genetic Algorithm)を勉強中。
Scilab 5から、遺伝的アルゴリズムのライブラリが標準で付いてくるようになった。
http://wiki.scilab.org/Genetic_algorithms
ネット上に情報はほとんどないので、Scilabのヘルプ
help optim_ga
と、サンプルプログラム(
ヘルプにのっているサンプルプログラムには些細なバグがある。
修正したうえで、グラフを表示させてみた。
（Scilab 5.0.1で確認)
deff('y=f(x)','y = sum(x.^2)');
PopSize = 100;
Proba_cross = 0.7;
Proba_mut = 0.1;
NbGen = 10;
NbCouples = 110;
Log = %T;
nb_disp = 10; // Nb point to display from the optimal population
pressure = 0.05;
ga_params = init_param();
// Parameters to adapt to the shape of the optimization problem
ga_params = add_param(ga_params,'minbound',[-2; -2]);
ga_params = add_param(ga_params,'maxbound',[2; 2]);
ga_params = add_param(ga_params,'dimension',2);
ga_params = add_param(ga_params,'beta',0);
ga_params = add_param(ga_params,'delta',0.1);
// Parameters to fine tune the Genetic algorithm. All these parameters are optional for continuous optimization
// If you need to adapt the GA to a special problem, you
ga_params = add_param(ga_params,'init_func',init_ga_default);
ga_params = add_param(ga_params,'crossover_func',crossover_ga_default);
ga_params = add_param(ga_params,'mutation_func',mutation_ga_default);
ga_params = add_param(ga_params,'codage_func',coding_ga_identity);
ga_params = add_param(ga_params,'selection_func',selection_ga_elitist);
//ga_params = add_param(ga_params,'selection_func',selection_ga_random);
ga_params = add_param(ga_params,'nb_couples',NbCouples);
ga_params = add_param(ga_params,'pressure',pressure);
Min = get_param(ga_params,'minbound');
Max = get_param(ga_params,'maxbound');
x0 = (Max - Min) .* rand(size(Min,1),size(Min,2)) + Min;
[pop_opt, fobj_pop_opt, pop_init, fobj_pop_init] = optim_ga(f, PopSize, NbGen, Proba_mut, Proba_cross, Log, ga_params);
disp(pop_opt(1));
// Plot 3D Surface and Optimized Point
x=Min(1):0.1:Max(1);
y=Min(2):0.1:Max(2);
z=[];
for k=1:length(x)
for m=1:length(y)
z(m,k)=f([x(k);y(m)]);
end
end
clf
plot3d(x,y,z);
p=pop_opt(1);
param3d(p(1),p(2),f(p));
e=gce(); //the handle on the 3D polyline
e.mark_mode='on';
e.mark_size=3;
e.mark_foreground=color('red');
実行結果
optim_ga: Initialization of the population optim_ga: iteration 1 / 10 - min / max value found = 0.005286 / 0.842114 optim_ga: iteration 2 / 10 - min / max value found = 0.000859 / 0.189952 optim_ga: iteration 3 / 10 - min / max value found = 0.000115 / 0.043110 optim_ga: iteration 4 / 10 - min / max value found = 0.000076 / 0.014128 optim_ga: iteration 5 / 10 - min / max value found = 0.000020 / 0.003757 optim_ga: iteration 6 / 10 - min / max value found = 0.000011 / 0.001156 optim_ga: iteration 7 / 10 - min / max value found = 0.000000 / 0.000396 optim_ga: iteration 8 / 10 - min / max value found = 0.000000 / 0.000134 optim_ga: iteration 9 / 10 - min / max value found = 0.000000 / 0.000041 optim_ga: iteration 10 / 10 - min / max value found = 0.000000 / 0.000011 0.0001984 0.0000989
ちなみに解析解は[0 0]。
Wavelab for scilab
WaveLabをScilabで試してみた。
インストールは特に必要ない。
Wavelab for scilab
http://www.scilab.org/contrib/index_contrib.php?page=displayContribution&fileID=934
をダウンロードし、展開。
Scilabを実行し、
exec loader.sce
とすると、Wavelabが組み込まれる。
実行例
//デモ用信号生成 sig = MakeSignal('Doppler',256); //CWT実行 nvoice=12; oct = 2; scale = 4; obj1 = CWT( sig, nvoice, 'Morlet',oct,scale); //軸のラベルを生成 [n nscale] = size(obj1); xtix = linspace(0,n,n); ytix = linspace(1+(oct-floor(log2(scale))),log2(n)-oct,nscale); skipx = floor(length(xtix)/8); skipy = floor(length(ytix)/nvoice); ha = gca(); ha.axes_visible="on"; def=format(); format('v',5); ha.x_ticks=tlist(['ticks' 'locations' 'labels'],1:skipx:length(xtix),string(xtix(1:skipx:$))); ha.y_ticks=tlist(['ticks' 'locations' 'labels'],1:skipy:length(ytix),string(ytix(1:skipy:$))); tmp=['v','e']; format(tmp(def(1)),def(2)); //画像化 subplot(211); ImageCWT(obj1, 'Overall', 'jetcolormap(256)') subplot(212) plot2d(sig);
ImageCWTがあまり良くないので、いずれ自作しようと思う。
