2019年12月17日火曜日

CORDIC sample

/
/* Portlist
CORDIC_SC cordic
(.Clock(Clock),.Reset(Reset),.Start(),.Theta(),
.SinOut(),.CosOut(),.Busy(),.End());
*/

module CORDIC_SC
#(
////////////////////////////////////////////////////
// Parameters
parameter bw_theta = 9,
parameter bw_out = 16
)
(
////////////////////////////////////////////////////
// Ports
input Clock,Reset,Start,
input [bw_theta-1:0] Theta,
output reg [bw_out-1:0] CosOut,SinOut,
output reg Busy,End
);
localparam fxp = 4; //小数点位置
localparam convert_vector = 19898; //変換ベクトル
localparam bw_tromaddr = 4;
localparam bw_tromq = 13;

////////////////////////////////////////////////////
// Registers
reg signed [bw_out:0] rX,rY; //符号付
reg [bw_theta-1:0] rInTheta; //入力Thetaを格納するだけ
reg signed [bw_tromq+1:0] rTheta; //Thetaを計算した結果を格納する
reg [bw_tromaddr-1:0] rTROMAddr; //Startを遅延させる
reg rDStart;

////////////////////////////////////////////////////
// Net
wire signed [bw_out:0] wResAsX,wResAsY; //AddSubの結果を出力
wire signed [bw_out:0] wSX,wSY; //X,Yをシフトした数値
assign wSX[bw_out] = wSX[bw_out-1];
assign wSY[bw_out] = wSY[bw_out-1];

wire signed [bw_tromq+1:0] wResAsTheta;
wire [bw_tromq+1:0] wTROMQ;
assign wTROMQ[bw_tromq+1:bw_tromq] = 2'b0;

wire [fxp-1:0] wFxp = 1'b0<<fxp;
wire signed [bw_tromq+2:0] wInTheta = {2'b0,rInTheta,wFxp};

wire wSeqEn = (rTROMAddr < bw_tromq) && Busy;
wire wEnd = Busy && !wSeqEn;

////////////////////////////////////////////////////
// Instantiations

//Address Depth:13 Q BitWidth:13 4bit FixPoint
DISTROM_CORDIC_THETA theta_rom (.Address(rTROMAddr), .OutClock(Clock),
.OutClockEn(1'b1), .Reset(Reset), .Q(wTROMQ[bw_tromq-1:0]));

//算術右シフトレジスタ
wire wRSStart = rDStart | wRSEnd;
SRightShifter2 #(.bw_in(bw_out)) rs2(
.Clock(Clock),
.Reset(Reset|Start), //Startでリセット
.Start(wRSStart),
.IN1(wResAsX[bw_out:1]),
.IN2(wResAsY[bw_out:1]),
.Amount(rTROMAddr),
.OUT1(wSX[bw_out-1:0]),
.OUT2(wSY[bw_out-1:0]),
.Busy(wRSBusy),
.End(wRSEnd)
);

wire wCompRes = rTheta > wInTheta; //Thetaを比較計算したθが大きければ１
AddSub #(.width(bw_out+1)) asX(
.SubEn(~wCompRes),
.A(rX),
.B(wSY),
.S({1'bz,wResAsX})
);
AddSub #(.width(bw_out+1)) asY(
.SubEn(wCompRes),
.A(rY),
.B(wSX),
.S({1'bz,wResAsY})
);
AddSub #(.width(bw_tromq+2)) asT(
.SubEn(wCompRes),
.A(rTheta),
.B(wTROMQ),
.S({1'bz,wResAsTheta})
);

always@(posedge Clock or posedge Reset) begin
if(Reset)begin
rX <= 0;
rY <= 0;
CosOut <= 0;
SinOut <= 0;
Busy <= 0;
End <= 0;
rTheta <= 0;
rInTheta <= 0;
rTROMAddr <= 0;
rDStart <= 0;
end else begin
rDStart <= Start;
if(Start) begin
rX <= convert_vector;
rY <= convert_vector;
end else if (wRSEnd) begin
rX <= wResAsX;
rY <= wResAsY;
end

if(Start) begin
rTheta <= wTROMQ;
end else if (wRSEnd) begin
rTheta <= wResAsTheta;
end

if(Start) begin
rInTheta <= Theta;
Busy <= 1'b1;
end else if (wEnd) begin
CosOut <= wResAsX[bw_out-1:0];
SinOut <= wResAsY[bw_out-1:0];
Busy <= 1'b0;
end

if((Busy || rDStart) && !wRSBusy)
rTROMAddr <= rTROMAddr + 1;
else if (wEnd)
rTROMAddr <= 0;

if(End)
End <= 1'b0;
else if (wEnd)
End <= 1'b1;

end
end
endmodule

/* Portlist
RightShifter2 rs2(.Clock(Clock),.Reset(Reset),.Start(),.IN1(),.IN2(),
.Amount(),.OUT1(),.OUT2(),.Busy(),.End());
*/
//符号付
module SRightShifter2
#(
parameter bw_in = 15
)
(
input Clock,Reset,Start,
input [bw_in-1:0] IN1,IN2,
input [3:0] Amount,

output reg [bw_in-1:0] OUT1,OUT2,
output reg Busy,End
);

//シフトレジスタを使って任意桁の右ビットシフト
reg [3:0] rCount,rAmount;
reg [bw_in-1:0] rSR1,rSR2;

wire wSeqEn = (rCount < rAmount);

always@(posedge Clock or posedge Reset) begin
if(Reset)begin
rCount <= 0;
rSR1 <= 0;
rSR2 <= 0;
End <= 0;
Busy <= 0;
rAmount <= 0;
OUT1 <= 0;
OUT2 <= 0;
end else begin
if(Start) begin
rSR1 <= IN1;
rSR2 <= IN2;
rCount <= 0;
rAmount <= Amount;
end else if(wSeqEn) begin
rSR1 <= {rSR1[bw_in-1],rSR1[bw_in-1:1]};
rSR2 <= {rSR2[bw_in-1],rSR2[bw_in-1:1]};
rCount <= rCount + 1;
end else if (!wSeqEn && Busy) begin
OUT1 <= rSR1;
OUT2 <= rSR2;
end

if(Start)
Busy <= 1'b1;
else if (!wSeqEn)
Busy <= 1'b0;

if(End)
End <= 1'b0;
else if (!wSeqEn && Busy)
End <= 1'b1;
end
end
endmodule

CORDIC １

https://scrapbox.io/craftmemo/Python_%E3%81%A7_CORDIC
抜粋

Python で CORDIC

CORDICは関数電卓等で古くから使われてきたそうで，その一番の特徴はアルゴリズム内に乗算を必要としない点である．
それはそうと，何に使えるのかというと，以下の近似計算ができる．
三角関数([$ \cos{z}],[$ \sin{z}])
極座標変換([$ \sqrt{x^2+y^2}],[$ \tan^{-1}{\frac{x}{y}}])
しかもこの2項目の2つの計算についてそれそぞれの計算は同時に計算可能である．
つまり，FPGAなど並列演算が可能なデバイスを用いれば完全に同じタイミングで数値を算出でき，複素計算をするとき等に役立つ．
アルゴリズムの解説は，「CORDIC」で検索すれば丁寧な解説が出てくるのでここでは割愛して，整数で計算するときの値の変換方法に重点を置きたい．
ただ，私の理解としては，角度を二分探索的に収束させていく手法だと理解している．

[* 定数について]
定数について，ひとつなぜかどこにも説明がないもの([$ K])があったため，自分で導出した．
CORDICで N bits の近似値を得たい場合，N個の[$ \tan^{-1}{\frac{1}{2^n}},n=0\cdots (N-1)]と定数 [$ K]が必要である．この[$ K]は以下の式で表される．
[$ K=k_{N-1}, k_{n+1}=\sqrt{{k_{n}}^2+\left(k_{n}\frac{1}{2^n}\right)^2},k_{0}=1, n=0\cdots N-1 ]
いきなりわけわからんと思うが，直角三角形の一辺が[$ k_0=1]でそれに直角な辺が[$ k_0\frac{1}{2^{n=0}}=k_0=1]なら，
斜辺は
[$ \sqrt{{k_0}^2+k_0^2}=\sqrt{1^2+1^2}=\sqrt{1+1}=\sqrt{2}\approx1.41421356\cdots =k_{1}]
である．
要はピタゴラスの定理である．この三角形の斜辺に対して[$ \frac{1}{2^n}]倍の辺をもつ直角三角形をかく．
今度の斜辺は，
[$ \sqrt{{k_1}^2+\left(k_1\frac{1}{2}\right)^2}=\sqrt{2+\left(\sqrt{2}\times\frac{1}{2}\right)^2}=\sqrt{2+2\times\frac{1}{4}}=\sqrt{2.5}\approx1.581138830\cdots =k_{2}]
となる．
これを繰り返すと，[$ K=k_{N-1}\approx1.6467602, in N=12]が求まる．[$ N]が大きければ，だいたい似た数字になる．ここで説明した直角三角形の直角を構成する長辺側の斜辺との角度が[$ \tan^{-1}{\frac{1}{2^n}},n=0\cdots (N-1)]というわけである．
[WolframAlpha https://www.wolframalpha.com/input/?i=x(n%2B1)%3D(x(n)^2%2B((x(n)*(1%2F(2^n)))^2))^(½)]で数列のグラフや表が出てくる（便利）．以下は自分が表計算で算出したものである．
table:k_n
n k[n]
0 1.4142135623731
1 1.58113883008419
2 1.62980060130066
3 1.64248406575224
4 1.64568891575725
5 1.64649227871248
6 1.64669325427364
7 1.6467435065969
8 1.64675607020488
9 1.64675921113982
10 1.64675999637562
11 1.6467601926847
12 1.64676024176197
13 1.64676025403129
14 1.64676025709862
15 1.64676025786545
16 1.64676025805716
17 1.64676025810509

全く関係ないが，夏目漱石の「こゝろ」という小説にKという登場人物が出てくる．『精神的に向上心のない者は馬鹿だ』はとても印象的である．しかしながら，時には「精神的な向上心」から解き放つことにより得られる「向上」もあること，あまり根詰めても仕方がなく，本能に任せることも大切であろう．あまりオチを語りたくないが，実際Kはこの言葉によって自分の首を絞めることになる．．．
閑話休題，CORDICによってNbit精度に対して上で説明したN+1個の定数をメモリに格納して，アルゴリズムをN回実行することによって求められるのである．
ここで，FPGAで計算することを考える．
[** 角度の表現]
象限毎に上位2bitを割り当てる
具体的には例えば12bitなら
[$ 0 \mathrm{rad} = \mathrm{0x000}]
[$ \pi \mathrm{rad} = \mathrm{0x400}]
[$ 2\pi \mathrm{rad} = \mathrm{0x800}]
[$ 3\pi \mathrm{rad} = \mathrm{0xC00}]
[* 定数Kの整数値換算]

2019年10月6日日曜日

https://yolo.love/matplotlib/mplot3d/　　を参照

3Dグラフをmatplotlibのmplot3dで作成する方法

2018/10/01

今回はmatplotlibで3Dグラフを作成する方法を紹介していきます。matplotlibでは3Dグラフを生成することもできます。mplot3dが用意されており、3Dの散布図や等高線プロットなどを描画できます。

公式のチュートリアルがよくまとまっているので、詳細はこちらでご覧ください。本記事も公式を参考に進めていきます。

参考mplot3d tutorialMatplotlib

もくじ [非表示]

mplot3dの下準備

まずはインポート部分から見ていきましょう。一般的なmatplotlibのインポート（import matplotlib.pyplot as plt）に加えて、Axes3Dをインポートします。

空間だけ描画すると以下のようになります。ここにこれから描画していきましょう！

Python

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

注意

Axes3Dを使用していないように見えますが、ここでインポートされていないとKeyError: ‘3D’が出力されるので、要注意です。

Line plotsを描画してみる

先ほどの空間にline plotsを描画してみると、以下のようになります。parametric curveというみたいですね。

Python

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl

mpl.rcParams['legend.fontsize'] = 10

fig = plt.figure()
ax = fig.gca(projection='3d')
theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
z = np.linspace(-2, 2, 100)
r = z**2 + 1
x = r * np.sin(theta)
y = r * np.cos(theta)
ax.plot(x, y, z, label='parametric curve')
ax.legend()

plt.show()

Scatter plotsを描画してみる

散布図を3D上に表示することもできます。

Python

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def randrange(n, vmin, vmax):
return (vmax - vmin)*np.random.rand(n) + vmin

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

n = 100
for c, m, zlow, zhigh in [('r', 'o', -50, -25), ('b', '^', -30, -5)]:
xs = randrange(n, 23, 32)
ys = randrange(n, 0, 100)
zs = randrange(n, zlow, zhigh)
ax.scatter(xs, ys, zs, c=c, marker=m)

ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()

Surface plotsを描画してみる

surface plotsは以下のようなグラフになります。

Python

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

fig = plt.figure()
ax = fig.gca(projection='3d')

X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)

surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

ax.set_zlim(-1.01, 1.01)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

まとめ

様々な3Dプロットを見てきましたが、他にもたくさんの3Dグラフを描画できます。公式のページにはサンプルコード付きですべて載っているので、ぜひ見てみてください！

参考mplot3d tutorialMatplotlib

また、Numpyの知識も必要になると思うので、関連記事をこちらに載せておきます。ぜひチェックしてみてください！

numpyのrandomで生成できる乱数を総まとめ

２D　リアルタイムグラフ描画　python その２

matplotlibでリアルタイム描画

この記事は最終更新日から3年以上が経過しています。

モチベーション

センサーとかの値を確認したり，機械学習の学習状況とかそんなんをリアルタイムにグラフにプロットしながら見たいときがある．
pythonにはmatplotlibっていう便利なグラフ作成ライブラリがあるのでこれを使いたい．

既知の問題とか

普通matplotlibでグラフを作るときは

import matplotlib.pyplot as plt

"""
plotしたいデータを計算する処理
...
"""
plt.plot(data)
plt.show()

とかやる．
だけど，show()を呼ぶとプログラムがそこでブロックされてリアルタイムな描画ができない．
調べてみるとshow()ではなくてdraw()を使おう!みたいな記事とか，Tkinter(tkinter)とかQtとかのGUIライブラリと組み合わせた方法とかがでてくるんだけど，Windows環境だとグラフに一度でも触るとGUIが固まったりしてとてもじゃないけど使えなかったりした．
(たしかMacとかのUnix環境だったら使えたような気がする．)

pause()を使おう!!!

例の如くstackoverflowにそれ関係のやりとりがあった．
参考: real-time plotting in while loop with matplotlib

結論はplt.show()の代わりにplt.pause(interval)を使うということ．
下に簡単なsin関数を永遠にplotする例を載っける．

注意点などはコメントに書いたので，読んでください．

matplotlib_realtime_plot_example.py

# -*- coding: utf-8 -*-
"""
matplotlibでリアルタイムプロットする例

無限にsin関数をplotし続ける
"""
from __future__ import unicode_literals, print_function

import numpy as np
import matplotlib.pyplot as plt


def pause_plot():
    fig, ax = plt.subplots(1, 1)
    x = np.arange(-np.pi, np.pi, 0.1)
    y = np.sin(x)
    # 初期化的に一度plotしなければならない
    # そのときplotしたオブジェクトを受け取る受け取る必要がある．
    # listが返ってくるので，注意
    lines, = ax.plot(x, y)

    # ここから無限にplotする
    while True:
        # plotデータの更新
        x += 0.1
        y = np.sin(x)

        # 描画データを更新するときにplot関数を使うと
        # lineオブジェクトが都度増えてしまうので，注意．
        #
        # 一番楽なのは上記で受け取ったlinesに対して
        # set_data()メソッドで描画データを更新する方法．
        lines.set_data(x, y)

        # set_data()を使うと軸とかは自動設定されないっぽいので，
        # 今回の例だとあっという間にsinカーブが描画範囲からいなくなる．
        # そのためx軸の範囲は適宜修正してやる必要がある．
        ax.set_xlim((x.min(), x.max()))

        # 一番のポイント
        # - plt.show() ブロッキングされてリアルタイムに描写できない
        # - plt.ion() + plt.draw() グラフウインドウが固まってプログラムが止まるから使えない
        # ----> plt.pause(interval) これを使う!!! 引数はsleep時間
        plt.pause(.01)

if __name__ == "__main__":
    pause_plot()

これで色々捗る!!!

おまけ

上で載っけたstackoverflowにdrawnowというパッケージを使った方法も載っているので紹介する．
drawnowはpipで入るので

$ pip install drawnow

でOK．

使い方

matplotlib_drawnow_realtime_plot_example.py

# -*- coding: utf-8 -*-
"""
matplotlibでリアルタイムプロットする例

無限にsin関数をplotし続ける
"""
from __future__ import unicode_literals, print_function

import numpy as np
import matplotlib.pyplot as plt
from drawnow import drawnow, figure

def drawnow_plot():
    # drawnowパッケージに含まれているfigure関数を使って
    # figureオブジェクトを作る
    fig = figure()
    x = np.arange(-np.pi, np.pi, .1)
    y = np.sin(x)

    def draw():
        """
        drawnow関数に渡すplotに必要な処理を記述
        plotに必要な変数はこの関数がアクセスできるスコープの変数にしておく
        """
        plt.plot(x, y)

    while True:
        x += .1
        y = np.sin(x)
        drawnow(draw)

...書いてからソースコード見てみたが，drawnowは結局中でdraw()を呼んでるみたいで，Windowsでやると上手く行かなかった(GUIが固まる現象)...

結論

リアルタイムにプロットしたいときはplt.pause()を使おう!!!

２D　グラフ　リアルタイム表示

@nv-h

2019年07月10日に更新

PythonでGUI画面とリアルタイムグラフ表示する

Python matplotlib Tkinter

前提

急にGUIを作れと言われた。ボタンとかでいろいろ制御できるリアルタイムグラフ描画をしたいが、性能は問われない。
gnuplot触ったことないし、C#もよくわからんが、Pythonは触ったことあるのでmatplotlibを使ってみる。

難しいこととか変な(?)ライブラリへの依存はしたくないのでPython3系のtkinterを使う。

表示したいのは以下。

折れ線グラフ複数
カラーマップと対応するカラーバー
制御用のボタンとか

事前にpip install numpy matplotlibしておく必要がある。
(Ubuntuの場合はその前にsudo apt install python3-pip python3-tkしておく)

詳細

最初の一歩

ボタンとグラフを同時に表示する必要があるので、matplotlib.orgの検索窓にでtkinterと入れてみる。
するとEmbedding in Tkが一番上にある。

それから、リアルタイム描画ということでアニメーションする必要があるので、同様にanimationと入れてみる。
すると、今度は少し下のほうにAnimated line plotが見つかる。

今回はボタンなどを配置しつつアニメーションしたいので上記二つを組み合わせる必要がある。

Animated plot in Tk

Embedding in TkとAnimated line plotを合体して、少し整理した。

#!/usr/bin/env python3


# https://matplotlib.org/gallery/user_interfaces/embedding_in_tk_sgskip.html
# https://matplotlib.org/gallery/animation/simple_anim.html


import tkinter

from matplotlib.backends.backend_tkagg import (
    FigureCanvasTkAgg, NavigationToolbar2Tk)
from matplotlib.figure import Figure
import matplotlib.animation as animation

import numpy as np


def _quit():
    root.quit()     # stops mainloop
    root.destroy()  # this is necessary on Windows to prevent
                    # Fatal Python Error: PyEval_RestoreThread: NULL tstate


def init():  # only required for blitting to give a clean slate.
    line.set_ydata(np.sin(x))
    return line,


def animate(i):
    line.set_ydata(np.sin(x + i))  # update the data.
    return line,


root = tkinter.Tk()
root.wm_title("Embedding in Tk anim")

fig = Figure()
# FuncAnimationより前に呼ぶ必要がある
canvas = FigureCanvasTkAgg(fig, master=root)  # A tk.DrawingArea.

x = np.arange(0, 3, 0.01)  # x軸(固定の値)
l = np.arange(0, 8, 0.01)  # 表示期間(FuncAnimationで指定する関数の引数になる)
plt = fig.add_subplot(111)
plt.set_ylim([-1.1, 1.1])
line, = plt.plot(x, np.sin(x))

ani = animation.FuncAnimation(fig, animate, l,
    init_func=init, interval=10, blit=True,
    )

toolbar = NavigationToolbar2Tk(canvas, root)
canvas.get_tk_widget().pack()

button = tkinter.Button(master=root, text="Quit", command=_quit)
button.pack()

tkinter.mainloop()

ボタンとかは、普通にtkinterで配置できる。

見た目

注意点

animation.FuncAnimation()で返ってくるオブジェクト(上記コードではani)の参照がない(?)とグラフが更新されないことがある。
クラスにした場合は、self.aniとかに格納して大事に持っておく必要がある。
(上記コードでは、ani = ...がないとだめかもしれない。)
set_ydata()とかset_data()では、データ以外の表示は更新されない。
Y軸範囲を自動追従したい場合とかは、自前でset_ylim()する必要がある。
blit=Trueだとグラフが更新されない(軸ラベルだけの更新とかが反映されない)ことがある。
ボタンでX軸表示モード変えたいのに変わらない時はblit=Falseにしてみるといいかもしれない。

登録: 投稿 (Atom)

2019年12月17日火曜日

CORDIC １

2019年10月6日日曜日

mplot3dの下準備

Line plotsを描画してみる

Scatter plotsを描画してみる

Surface plotsを描画してみる

まとめ

２D リアルタイムグラフ描画 python その２

matplotlibでリアルタイム描画

モチベーション

既知の問題とか

pause()を使おう!!!

おまけ

使い方

結論

２D グラフ リアルタイム表示

PythonでGUI画面とリアルタイムグラフ表示する

前提

詳細

最初の一歩

Animated plot in Tk

見た目

注意点

２D　リアルタイムグラフ描画　python その２

２D　グラフ　リアルタイム表示