■技術
AD変換回路を使ったRaspberry Piのセンサ信号読み込み
LTsipiceを使えるようにする(LDのAPC回路シミュレート)
小規模発注システムの市場調査(Webアプリ開発向け)
Interfaceを読みつつ画像認識による物仕分けモデル計画確認
JavaScriptによるリアルタイムデータ表示スクリプト作成
MatlabでTektronicsのオシロスコープからUSBで信号取り込む方法確認
サーボ制御の基本を勉強
PDとアンプを組み合わせた光検出回路の基礎を勉強
スペクトラムアナライザの基本を勉強
USB機器の通信基礎の勉強
arduinoとPCのシリアル通信
温度計設計検討
フィードバック制御の復習と実験系の検討

■その他の学び
レアジョブ
ロシア語
バイオリン
キーボード
硬筆
デッサン

■日常
部屋の掃除
観葉植物探し
会社用ニットの長袖かベスト探し
カジュアル用白シャツ探し
ダニエルボブの鞄でいいものがないか探す
香水探し
音読
macのノート購入検討

2017-03-05

Matlabによる信号処理【フィルタの作り方】

Matlab 情報処理

DSP(Digital Signal Processing)でフィルタを使う訳ですが、今の知識だとFIRとIIRどちらを選べばいいのかわからないのでお勉強します。

【記事1】のNI様のところに分かりやすく書いてあります。

【記事1】
IIRフィルタとFIRフィルタ - DIAdem 2012ヘルプ - National Instruments

とりあえず、今後の基準は以下の通りと考えておきます。

■IIR
メリット: 係数が少ない/計算ボリュームが小さい
デメリット: 非線形位相応答
得意とする情報: 振幅

■FIR
メリット: 非線形位相応答
デメリット: 係数が多い
得意とする情報: 位相

こんなところです。

Matlabのドキュメンテーションに【記事2】があるので、基本的には記事にそって進めていきます。

【記事2】
jp.mathworks.com

まずは信号を設定します。信号はMatlabのサンプルデータとして用意されているchirpを使います。信号は以下の通りです。

f:id:cocosuzu:20170304124209p:plain

さて、フィルタを作ります。

blo = fir1(34,0.48,chebwin(35,30));

がフィルタの係数を生成する関数です。fir1はウィンドウを指定してフィルタの係数を生成する関数で、引数の最後にはウィンドウを定義するchebwinがいます。この辺はひとまずwindow関数が選べますよという程度にとどめておきます。

続いてフィルタを適用するためにfilter関数に処理したいデータとフィルタの係数セットを渡します。

すると時間波形は以下の通りとなります。

f:id:cocosuzu:20170305133126p:plain

大きな信号が消えました。ぱっと見で変化は分かりますが、より分かりやすくするために周波数スペクトルを見てみましょう。

f:id:cocosuzu:20170305164639p:plain

上段がフィルタ無しの生データ、下段がフィルタを適用した後の信号です。2000Hz以上の主ピークたちがごっそりいなくなっています。

これでフィルタが適用できていることがよく分かります。

ここで一歩戻って今回どんなフィルタを作ったのか見てみましょう。freqz関数でフィルタ係数の特性を見ることが出来ます。先ほど作ったbloを見てみると…

freqz(blo,1)

f:id:cocosuzu:20170305164351p:plain

どこにも2000Hzとは書いてません。正規化周波数で書かれています。ただし、ここで言う正規化周波数は書籍などでよく見るサンプリング周波数を1とする正規化周波数ではなく、ナイキスト周波数を1としています。それで正規化周波数の単位には2 $\pi$ rad/Samplingではなく $\pi$ rad/Samplingと書かれているのですね。

サンプリング周波数は8192Hzであることから、その半分である4096Hzを1とするとfir関数で引数としていた0.48は1966Hzに相当します。この1966Hzがローパスフィルタの周波数となっているわけですね。

さて、大体【記事2】の内容が分かってきたところでフィルタを自由自在に作る方法に移っていきます。

ここで使うのがFDA(Filter Design and Analisit) toolです。これはフィルタの形状を見ながらちょちょいとフィルタを作るには便利なツールです。

先ずは、ここまでで扱ってきたchirpデータの周波数スペクトルを見直してみます。

f:id:cocosuzu:20170305170509p:plain

バンドパスフィルタで主ピークの両端を削るようなフィルタを作ってみましょう。となるとバンドパスフィルタの周波数を2700～3500Hzにします。

f:id:cocosuzu:20170305170750p:plain

こんな感じの画面が出てきますので、以下のように設定します。細かい数値は適当です。

f:id:cocosuzu:20170305171614p:plain

エクスポートしてワークスペースに係数ベクトルをNumとして出力します。これを先のbloの代わりに代入してあげればバンドパスフィルタとして処理されるはずです。

ただ、メインのスクリプトがclear allで変数を都度リセットするようにしているので、都度Numを生成してくれるように作成したフィルタをmファイルとして出力しておきます。Matlabコードを生成を選ぶとmファイルを作ってくれます。

% All frequency values are in Hz.
Fs = 8192;  % Sampling Frequency

Fstop1 = 2500;            % First Stopband Frequency
Fpass1 = 2700;            % First Passband Frequency
Fpass2 = 3500;            % Second Passband Frequency
Fstop2 = 3700;            % Second Stopband Frequency
Dstop1 = 0.001;           % First Stopband Attenuation
Dpass  = 0.057501127785;  % Passband Ripple
Dstop2 = 0.0001;          % Second Stopband Attenuation
dens   = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 1 ...
                          0], [Dstop1 Dpass Dstop2]);

% Calculate the coefficients using the FIRPM function.
bpf  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

ここで、Fsはスクリプト内で共通なので削除しておいて、生成されたフィルタ係数をfreqzで見てみましょう。

f:id:cocosuzu:20170305172558p:plain

ちゃんと出来てますね。それでは先のbloの代わりにbpfを使いましょう。先ずは時間波形を見てみます。

f:id:cocosuzu:20170305172823p:plain

これだけだと生データと大きくは変わらないですね。続いて生データとBPF処理後の周波数スペクトルです。

f:id:cocosuzu:20170305174018p:plain

上段が生データ、下段がBPF後の周波数スペクトルです。

ちゃんと設計したとおり2500以下は-60dB、3700以上は-80dBになるようになっています。以上、Matlabにおけるフィルタ設計の基本でした。

結局FIRとIIRはまだどっちを選べばいいか聞かれたら答えられません。追々分かったら書いていきたいと思います。

2017-03-05

最近買って便利だったもの

日常

絵を描く用のイーゼルがホワイトボードを立てておくのに便利です。

f:id:cocosuzu:20170305153355j:plain

https://www.amazon.co.jp/gp/aw/d/B001WB3K7G/ref=yo_ii_img?ie=UTF8&psc=1

ちょっと剛性は弱いですが、うちの壁は画鋲とか刺さらないので、コンパクトで自立式となるとこの辺がいいかもしれないです。

2017-03-04

今週の成績

投資

f:id:cocosuzu:20170304094901p:plain

先週は以下の通りです。

f:id:cocosuzu:20170224203939p:plain

先週比+2.43%です。

ビットコインの相場はというとこんな感じです。

f:id:cocosuzu:20170304111749p:plain

1月の暴落は過熱からのトレンド戻し&オーバーシュートと考えられますので、16年8月くらいから上昇トレンドを維持しつつ、上昇率も徐々に増加傾向にあります。今も若干ながらトレンドより上にいるので、若干の押し目はあるかもしれません。10~20k円ほどの急落はあるかも知れませんが、このまま上昇していくのではないかという印象を受けます。このあたりはまだ数値的な売買条件を設けていないのですが、長期で上昇するものだとは思っています。それがいつまで、どのくらいになるのかを見積もるためにもお金の仕組みと歴史、ビットコインの市場動向を調べる必要があります。

ちなみに、現在保有しているビットコインは21日前に取引を始めて+17%となっています。ちなみに会社の先輩は平均買い付け価格9万円なので、167%とかなりの運用成績になっています。車買えるくらいはつぎ込んでると言ってたので500万とすると…833万の+333万！ワオ！

そして先輩は100、1000倍とかになるんじゃない？と言ってるので、5億円！50億！

と、取らぬ狸の皮算用な上に人の話というどうしようもない感じですが、そんなドリームもあるということです。

私はまだまだ投入資金の規模が小さいので、せっせと買い増ししていきたいところです。

そうそう。トランプさんの演説見てなかったので見てみました。

www.youtube.com