これは、1年くらい前に、Rを使って作った。
この時、メインで使用していたのは、Scilabだったが、文字のプロットする機能がScilabで見当たらなかったので、Rを使った。
Rは、統計解析向けのプログラミング言語及びその開発実行環境だけれども、文字をプロットできそうだという理由だけで使った。
scilabのプロット関連の関数から、textを表示するものが見当たらなかったので、
コレを、Matplotlibを使って、もう一度やろうという試み。
以前は、
1 変な数字の配列を使って、一個一個の数字が分かる、畳み込み動画
2 実際に101タップほどのFIRフィルタを畳み込んで、1,10,100Hzを足した波形から、10HzだけをBPフィルタを取り出す
の順で説明して、何をやっているか理解してもらっていた。
式は、よく分からないので、飛ばしてしまっていた。
これは、BPフィルタを畳み込んで、ノイズ処理をしている様子だが、「これも畳み込みで実現しています。」と言わなければ、分かってもらえない様な細かいプロットになってしまっている。
かなり頑張って作ったのに。
今回は、
2 全ての四則演算を、人が判る速度でアニメする。
3 さらに、公式の理解への準備等を載せる。
全部、大変なので、出来なかった。
いまはまだ、テキストプロットがおぼつかない。
1については、下の動画を見ながら、考えてほしいんだけど、動いてるやつが左右対称になっている所を利用して、掛け算の回数をなるべく減らそうっていう考え方で、掛け算減らすといいことづくめらしいよという。。。
なんで掛け算減らすといいのかというと、コンピュータの中の積算回路は、足し算の回路よりも数が少なくて、節約しないとすぐに一杯になるからで、リアルタイム処理するDSPとかFPGAとかマイコンとかでは、これがボトルネックになって、性能を決めてしまうらしい。
自分の作ったFIR&FFTの、tap数って所をいじってもらうと分かると思うけど、
http://soft.sakuraweb.com/firandfft/
tap数は高ければ高いほど、周波数特性が綺麗になる。
周波数特性を見ると、同じカットオフ周波数でも、タップ数の高い場合の方が、カットオフ周波数の周辺のスロープが急峻になる。
つまり、任意の周波数の部分だけを取り出すことが出来る。
じゃあ、タップ数あげまくってみればと思いきや、tap上げると、計算に時間がかかる。
タップ数(オーダー数ともいう)は、掛け算する関数の長さだから、当然計算が多くなっている。
FIRでノイズフィルタやるときは知っておきたい豆知識。元ネタは、EDN japanか何かで見た、が忘れた。
左右対称だと何でいいのか。
フィルタ対象のオリジナル信号 1 2 3 4 5 6 7
左右対称のフィルタ関数 1 3 5 7 5 3 1
これを全て、掛け算すると、
フィルタ対象のオリジナル信号 1 2 3 4 5 6 7
* * * * * * *
左右対称のフィルタ関数 1 3 5 7 5 3 1
の7回。7個も積算回路を占有してしまう。
下の、「左右対称のフィルタ関数」の3だけに注目すると、3は2個あるので、
3*2+3*6
掛け算2回。
これは、積和の公式というので整理できる。
教科書の最初過ぎて、当たり前過ぎて普段意識しない。
意識しないから、盲点になる系のアレ。
3*(2+6)
掛け算1回に減ってる。
積算回路使わなくて済む。
2倍早くするのか、2倍信号品質(音質)をよくするのかは、人次第だけど、とにかく、性能アップするコツ。
左右対称なので、約1/2に減る。
なぜ約なのかというと、タップ数は、絶対奇数じゃないとだめだから。
python的に書くなら、
floor(tap/2)+1
個分の、タップ数になる。
実はそこら辺よく分かってない。
単に使っている関数が、奇数じゃないとエラーを返してくるから程度の認識しかない。
フィルタ関数を作ってプロット。
戻ってきた。
その上に、文字列として数字をプロット。
視認性とか考えると、50px2桁で10個が限度ではないだろうか。
いや、計画しても無駄だ、手を動かそう。
import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
n = 21
tap = 11
nyquist_frq = n/2
lower_cutoff_frq = 5.0
higher_cutoff_frq = 8.0
window_type = "hamming"
fir_filtre = scipy.signal.firwin(tap, [lower_cutoff_frq, higher_cutoff_frq], pass_zero = False, nyq = nyquist_frq, window = window_type)
str_fir_filtre=[(str(np.round(fir_filtre[i],2))) for i in range(tap)]
f = 1.5
#Ideal sin curve
x1 = np.linspace(0, 1, 2001)
theta1 = np.linspace(0, 2*np.pi*f, 2001)
sin_curve1 = np.sin(theta1)
#Discrete sin curve
x2 = np.linspace(0, 1, n)
theta2 = np.linspace(0, 2*np.pi*f, n)
sin_curve2 = np.sin(theta2)
start_for_x3 = (tap/2.)/(n)
for i in range(n):
f = plt.figure()
ax = plt.subplot(111)
plt.plot(x1, sin_curve1, alpha = 0.5)
plt.plot(x2, sin_curve2, "o")
x3 = np.linspace((-start_for_x3+(i)/float(n-1)), (start_for_x3+(i)/float(n-1)), tap)
for m in range(tap):
plt.text(x3[m], fir_filtre[m], str_fir_filtre[m], fontsize = 4, backgroundcolor=(1,1,1,0.5))
ax.spines['bottom'].set_position(('data',0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))
ax2 = plt.twinx(ax)
ax2.tick_params(right='off', labelright='off')
plt.xlim(-0.3, 1.3)
plt.ylim(-1,1)
for label in ax.get_xticklabels() + ax.get_yticklabels():
label.set_fontsize(16)
label.set_bbox(dict(facecolor='white', edgecolor='None', alpha=0.65 ))
plt.tight_layout()
f.set_size_inches(7.2,4.05)
f.savefig("string plot"+"%04.f"%(i)+".png", dpi=266, facecolor='w', edgecolor='w',
orientation='portrait', papertype=None, format=None,
transparent=False, bbox_inches=None, pad_inches=0.1,
frameon=None)
plt.clf()
【pythonの最新記事】
- 中速フーリエ変換 ~離散フーリエ変換より..
- 断面二次モーメントを、座標点の配列から計..
- 断面二次モーメントを、座標点の配列から計..
- fontファイルの文字データ(グリフ)を..
- matplotlibのpyplot.pl..
- 計算力学技術者試験の問題集 自炊(裁断→..
- pythonで、ホワイトノイズやピンクノ..
- 脳ドッグに行ってきた。→MRIの画像デー..
- matplotlibのimshowで円を..
- matplotlibの、cmapを、徐々..
- matplotlibのmake_axes..
- matplotlib floatinga..
- matplotlib plotの色を、値..
- Pythonで、「二次元フーリエ変換した..
- matplotlibのlinestyle..
- どちらが正しいRGBか。(matplot..
- matplotlibのannotateの..
- matplotlibで、x軸とy軸の数字..
- VBAで、pythonのrangeとか、..
- matplotlibのaxes3Dで、a..