FFTの練習

FFT についてわかりやすいサイトA DFT and FFT TUTORIAL


FFT のバタフライダイアグラムのコードをみても、いまいち何をやって
いるのかわからない。なので、自分で書いてみました。
自分のコードとほかの人のコードはだいたい同じになったので、ここまで
やってやっと人のが読めそうな気がしてきました。手を動かすのは大事。

前述のサイトのバタフライダイアグラム をみながらコードにしてみました。
演算の中身はないです。画面に "回転因子はこれ" で "xとyをバタフライ"
というメッセージを出すだけ。

N=8 でのバタフライ演算(中身なし)のコード


output

0
twiddle(0, 2)
butterfly(0, 1)
2
twiddle(0, 2)
butterfly(2, 3)
4
twiddle(0, 2)
butterfly(4, 5)
6
twiddle(0, 2)
butterfly(6, 7)
0
twiddle(0, 4)
butterfly(0, 2)
twiddle(1, 4)
butterfly(1, 3)
4
twiddle(0, 4)
butterfly(4, 6)
twiddle(1, 4)
butterfly(5, 7)
0
twiddle(0, 8)
butterfly(0, 4)
twiddle(1, 8)
butterfly(1, 5)
twiddle(2, 8)
butterfly(2, 6)
twiddle(3, 8)
butterfly(3, 7)

VisualStudio2015 で glew-1.13.0 をビルド

1. glew のサイトからソースコード( glew-1.13.0.zip )をダウンロード
2. build\vc12\glew.sln を VisualStudio Community 2015 で開く
3. VC12 のソリューションなので変換
4. "VC++ コンパイラとライブラリをアップグレードします"ダイアログでOK
5. ビルド->バッチビルドを開いてすべてチェック
6. ビルド

であっさり完了。

bin, include, lib を置きたい場所にコピーすればOK

wxWidget を VS2015 でビルド

wxWidget を使ってみようと思いダウンロードすると VS2015 用の
プレビルドがなかったのでビルドしてみました。

まだ使っていないので、あっているかどうかわからないけど、ビルドは
通ったのでとりあえずメモ。 nmake のオプションは build\msw\config.vc
に説明が書いてあるのを参考にしました。

1. ソースコードダウンロード
git clone https://github.com/wxWidgets/wxWidgets.git

2. ビルドの設定が反映されるヘッダー？を作成
include/wx/msw/setup0.h を include/wx/msw/setup.h にリネーム

3. x86 版をビルド
VS2015 x86 Native Tools コマンドプロンプトで以下を実行



4. x64 版をビルド
VS2015 x64 Native Tools コマンドプロンプトで以下を実行

boost::property_tree の使い方メモ

Sony Camera Remote API の getEvent で取れる json のパースに
boost::property_tree を使うときに、いまいちよくわかってなかったので
メモ。

この api で返ってくる json は {{result:[{...},"null","",[...]}}
みたいなので、null になったり 空文字になったり状態によって変化する曲者。
なのでちゃんと処理しないと例外が出る。例外はデバッグのときに紛らわしい
ので、なるべく出さないようにしたいので、事前に判定したい。

ノードがあるかわからないときは optional で取得すれば例外を出さずに
判定できる。

ただの文字列("null" or "") か何か値が入っている ("{x:y}" or "[...]") かどうか
は child.size() をみればOK(たぶん)

code


output

e.first.data()=
e.second.data()=
child.get_optional("")=
child.size()=2
----
e.first.data()=
e.second.data()=null
child.get_optional("")= null
child.size()=0
----
e.first.data()=
e.second.data()=
child.get_optional("")=
child.size()=0
----
e.first.data()=
e.second.data()=
child.get_optional("")=
child.size()=1
----

const array メンバの初期化

メンバ関数のポインタを、配列でメンバ変数に持ちたい。
その配列は固定長なので const にしておきたい。
たまにしか使わないので、毎回つまづくからメモ。

boost::asio::streambuf にバイナリファイルを読み込む

boost1.59 を VS2013 でビルド

マルチスレッドDLL とかの設定がほかの lib と合わずエラーがでるので、boost を
ビルドしなおしました。

毎回面倒なので全種類作っておこうと思い、
- x86/x64
- static lib or dll
- static runtime link or dll runtime link

で 8 種類(debug/releaseあわせると16種類)いるのか！と思いましたが、
dll の boost に static の runtime は link できないみたい。
確かに別の dll や本体アプリで違う種類の runtime がリンクされてたら
インスタンスが2つ要るのでおかしなことになる。

リリース/デバッグビルドは "variant=debug,release"と書けば両方ビルドされるみたい。
あと "link=static,shared" もいけるとどこかにあったけど未確認。
使ったVSは2013 , boost1.59 です。

rem x86, static lib, static rt(MT/MTd) : xxx-vc120-mt-x_xx.lib, xxx-vc120-mt-gd-x_xx.lib
b2.exe --toolset=msvc-12.0 variant=debug,release link=static threading=multi runtime-link=static --stagedir=stage\x86\vc12\static\mt --without-python --without-mpi --build-type=complete address-model=32
rem x86, static lib, dynamic rt(MD/MDd) : xxxx-vc120-mt-x_xx.lib, xxxx-vc120-mt-gd-x_xx.lib
b2.exe --toolset=msvc-12.0 variant=debug,release link=static threading=multi runtime-link=shared --stagedir=stage\x86\vc12\static\md --without-python --without-mpi --build-type=complete address-model=32
rem x86, dynamic lib, static rt(MT/MTd) : この組み合わせはNG
b2.exe --toolset=msvc-12.0 variant=debug,release link=shared threading=multi runtime-link=static --stagedir=stage\x86\vc12\dynamic\mt --without-python --without-mpi --build-type=complete address-model=32
rem x86, dynamic lib, dynamic rt(MD/MDd) : xxx-vc120-mt-x_xx.{dll,lib}, xxx-vc120-mt-gd-x_xx.{dll,lib}
b2.exe --toolset=msvc-12.0 variant=debug,release link=shared threading=multi runtime-link=shared --stagedir=stage\x86\vc12\dynamic\md --without-python --without-mpi --build-type=complete address-model=32


rem x64, static lib, static rt(MT/MTd)
b2.exe --toolset=msvc-12.0 variant=debug,release link=static threading=multi runtime-link=static --stagedir=stage\x64\vc12\static\mt --without-python --without-mpi --build-type=complete address-model=64
rem x64, static lib, dynamic rt(MD/MDd)
b2.exe --toolset=msvc-12.0 variant=debug,release link=static threading=multi runtime-link=shared --stagedir=stage\x64\vc12\static\md --without-python --without-mpi --build-type=complete address-model=64
rem x64, dynamic lib, static rt(MT/MTd) : この組み合わせは NG
b2.exe --toolset=msvc-12.0 variant=debug,release link=shared threading=multi runtime-link=static --stagedir=stage\x64\vc12\dynamic\mt --without-python --without-mpi --build-type=complete address-model=64
rem x64, dynamic lib, dynamic rt(MD/MDd)
b2.exe --toolset=msvc-12.0 variant=debug,release link=shared threading=multi runtime-link=shared --stagedir=stage\x64\vc12\dynamic\md --without-python --without-mpi --build-type=complete address-model=64

hstファイルをcsvに変換

別にFXをやっている訳ではないんですが、いやちょっとはやってみたいとも思いますけど、そもそも元手が(10万とか)ないんで無理なのが悲しいところ。でもちょっと世界をのぞいてみたいと思いまして、為替相場の動きを観察しようと思います。

そのためにはまず過去のデータが必要です。2005年から1分ごとのデータが、FXDD という証券会社のメタトレーダーのヒストリカルデータというところからダウンロードできました。

さて、このデータ。メタトレーダー4というFX取引用のアプリのフォーマットのようです。
フォーマットの中身はMT4を使ってFXに解説されていました。

探せばごろごろあるのでしょうが、向学のために csv に変換するプログラムを書いてみました。本当は C# でサクッと行きたかったのですが、諸事情により VisualStudio をインストールしていないので、 xcode で c++ で書きました。

ソースコードは GitHub からダウンロードできます。

#include 
#include 
#include 
#include 

using namespace std;

typedef struct _HistoryHeader
{
    uint32_t    version;
    char        copyright[64];
    char        symbol[12];
    int32_t     period;
    int32_t     digits;
    uint32_t    timesign;
    uint32_t    last_sync;
    uint32_t    reserved[13];
} HistoryHeader;

#pragma pack(push,1)
typedef struct _RateInfo
{
    uint32_t    ctm;
    double      open;
    double      low;
    double      high;
    double      close;
    double      vol;
} RateInfo;
#pragma pack(pop)


void printHistoryHeader(const HistoryHeader& header)
{
    cout << "version   : " << header.version << endl;
    cout << "copyright : " << header.copyright << endl;
    cout << "symbol    : " << header.symbol << endl;
    cout << "period    : " << header.period << endl;
    cout << "digits    : " << header.period << endl;
    time_t tmp = header.timesign;
    cout << "timesign  : " << put_time(localtime(&tmp), "%F %T") << endl;
    tmp = header.last_sync;
    cout << "last_sync : " << put_time(localtime(&tmp), "%F %T") << endl;
}

void printRate(const RateInfo& rate, ostream& ost)
{
    time_t tmp = rate.ctm;
    ost << put_time(std::localtime(&tmp), "%F %T");
    ost << "," << rate.open;
    ost << "," << rate.low;
    ost << "," << rate.high;
    ost << "," << rate.close;
    ost << "," << rate.vol << endl;
}

int main(int argc, const char * argv[])
{
    if(argc < 2) {
        cout << "need to give hst file and output file" << endl;
        cout << "ex) hstconverter [hst file] [output file]" << endl;
        return 1;
    }
    
    ifstream ifs;
    ifs.open(argv[1], ios::in | ios::binary);
    if(!ifs) {
        cout << "file open error" << endl;
        return 1;
    }
    
    HistoryHeader header;
    ifs.read((char*)&header, sizeof(header));
    if(ifs.bad()) {
        cout << "file read error" << endl;
        return 1;
    }
    
    printHistoryHeader(header);
    
    ofstream ofs;
    ofs.open(argv[2], ios::out | ios::trunc);
    if(!ofs.is_open()) {
        cout << "failed to create a file : " << argv[2] << endl;
        return 1;
    }
    
    RateInfo rate;
    while(!ifs.eof()) {
        ifs.read((char*)&rate, 44);
        if(ifs.bad()) {
            cout << "file read error" << endl;
            return 1;
        }
        printRate(rate, ofs);
    }
    
    cout << "done" << endl;
    
    return 0;
}

#pragma pack でバイトアラインを無効にして、構造体に一気にデータを読み込んでいます。
これってバイトオーダーが違う処理系だったら大丈夫なのでしょうか。

そしてフォーマットにある time_t はなんと 32bit time_t のことの様なので、
uint32_t に置き換えて構造体に取り込んで、使うときに time_t にキャストしています。

次のように実行して USDJPY.hst を csv に変換して output.csv に書き出します。

hstconverter USDJPY.hst output.csv

できたファイルの内容はこんな感じ。

2005-01-10 11:31:00,104.79,104.79,104.79,104.79,5
2005-01-10 11:32:00,104.79,104.78,104.79,104.78,6
2005-01-10 11:33:00,104.78,104.77,104.78,104.77,5
...
2014-08-09 07:59:00,102.063,102.06,102.067,102.064,34
2014-08-09 07:59:00,102.063,102.06,102.067,102.064,34

さてこれを使ってどうするか。

自由切断を実装してみた

メタルギア ライジング リベンジェンスの自由切断がかっこよすぎるので、
自分でもなんちゃって実装してみました。

テクスチャを張ったキューブを自由切断したところ


ちょっと複雑な形の物体を自由切断したところ


断面に張るポリゴンも生成します。ただし断面のテクスチャは考慮していません(座標0,0になっています)。
また、穴が開いた断面には対応していません(外枠か、中の穴かどちらかに張るだけ)。


Vertex を用意

const GLfloat data[][8] =
{
   	{ -0.5f, -0.5f,  0.5f,  0.0f,  0.0f,  1.0f, 0.0f, 0.0f},
   	{  0.5f, -0.5f,  0.5f,  0.0f,  0.0f,  1.0f, 1.0f, 0.0f},
   	{  0.5f,  0.5f,  0.5f,  0.0f,  0.0f,  1.0f, 1.0f, 1.0f},
	...
};

平面を作成

MODELVEC3D p = MODELVEC3DMake(0.0f, 1.0f, 0.0f); // any point on the plane
MODELVEC3D n = MODELVEC3DMake(1.0f, 0.0f, 0.0f); // normal vector
MODELPLANE plane = MODELPLANEMake(p,n); // make a plane

結果を受け取るバッファを用意

GLfloat bufN[MAX_CHOP_BUF][8] = {0.0f}; // vertices located in normal vector sides
GLfloat bufA[MAX_CHOP_BUF][8] = {0.0f}; // vertices located in anti-normal vector sides
int bufNCount = 0;
int bufACount = 0;

自由切断実行

Chop(plane, &(data[0][0]), sizeof(data)/8/sizeof(GLfloat), bufN, bufA, bufNCount, bufACount);

処理内容は、ひたすら面と一つ一つのポリゴン(とどのつまり直線)の計算をしているだけです。

コードはGitHub においてあります。