nabeの雑記帳

検索条件

タグで絞り込みプログラム::Perl(1) プログラム::TCP/IP(3) プログラム::配布(1) メモ(2)

日付で絞り込み 1月(1) 2月(2) 5月(1) 7月(2) 8月(1)

全7件 (1/2ページ) 記事本文を表示

キット各種完売。
msBerryDAC改造情報

速くて使いやすいCMS - adiary / 頒布物一覧 / 寄付

2007/08/24(金)実数FFT/IFFT関数

メモプログラム

長年理解するのを拒否していたFFTとここ何週間か格闘しています。楽々とアルゴリズムを導出・実装できる人はいいのでしょうが、FFTのようなアルゴリズムをきちんと理解して実装するのは思うほど容易ではありません。既に導出されたアルゴリズム（や雛形サンプル）を何も考えず実装すれば簡単ですが、そういうことが生理的にできない場合、改良とかを考え出して簡単には実装できなくなるんですよね。

前置きはともかく、今日やっと実数専用のFFT/IFFTが作成できました。実数専用にすることで理論的には通常のFFTの半分の時間で処理できます。

rfft.cpp

見ての通り、実数FFTの対称性を利用したものすごく面倒くさく複雑な作りになっているんですが、有名なFFTの概略と設計法のソースと速度比較をしたら倍近く遅い……。64K点、double型、VIA C3 500MHzで動作させて、自作が100ms、リンク先のソースが60ms。

いかに細かい最適化をしても実装アルゴリズムの時点で差があると非常に大きいですね……。あとでアルゴリズムの差を検証予定。

ソース解読メモ

大浦氏のFFTソース(fft4g.c)解読メモ。つまりリバースエンジニアリングです。

サブルーチンの機能

cdft()	複素数FFT/IFFT。データサイズは N/2
rdft()	実数FFT/IFFT。内部的に複素FFTを呼び出し
makewt()	sin/cosテーブルを w[] に格納（ただし格納位置はビット反転）
makect()	cosテーブルを w[nw～] に格納（ただし格納位置はビット反転）
bitrv2()	ビット反転を実行
bitrv2conj()	ビット反転しつつ、読込データの複素共役を取る
cftfsub()	複素IFFT（f=forward, 回転子が正であるという意味）
cftbsub()	複素FFT（b=backward, 回転子が負であるという意味）
rftfsub()	実数IFFT
rftbsub()	実数FFT

IFFTをFFTで代用する方法

$C(k) = \sum_{n=0}^{N-1} c(n) W^{kn}$ 　　 $c(n) = \sum_{k=0}^{N-1} C(k) W^{-kn}$

なんだけども、複素共役を考えると

$\overline{c(n)} = \sum_{k=0}^{N-1} \overline{C(k) W^{-kn}} = \sum_{k=0}^{N-1} \overline{C(k)} W^{kn}$

つまり前処理と後処理としてデータの複素共役を取ってあげれば、同じFFTルーチンを使い回すことができる（効率を考えて、初段部分と終段部分だけ個別特別に実装すれば、中段部分は使い回せる）。

その他メモ

一般的なFFT(DFT)と比べ虚数部の符号が逆になっている。*1

周波数間引きFFTである

やっとこのプログラムの要が理解出来ました。

先にデータをスクランブルして（並べ替えて）いるため「時間間引きFFT」に見えますが、実際には周波数間引きFFTです。sin/cosテーブルもわざわざスクランブル位置に格納しています。

こうすることでFFT実行時のデータアクセスをシーケンシャルに行え*2、それが功を奏して実行速度が飛躍的に速くなっています。またsin/cosテーブル、偶数番地に cos、奇数番地に sin を格納し、全体としても π/2 しか用意しないことで、データアクセス量を減らしキャッシュが効きやすくなっています。

256点の複素FFTを実行させたときの、cos/sinテーブル参照位置およびデータ参照位置は次のようになります。

 call cftmdl(512, 8)
  k1=1 (0.923880,0.382683), k2=2 (0.707107,0.707107)
    32 36 40 44
    33 37 41 45
    34 38 42 46
    35 39 43 47
  k1=2 (0.980785,0.195090), k2=4 (0.923880,0.382683)
    64 68 72 76
    65 69 73 77
    66 70 74 78
    67 71 75 79
  k1=3 (0.831470,0.555570), k2=6 (0.382683,0.923880)
    96 100 104 108
    97 101 105 109
    98 102 106 110
    99 103 107 111
  k1=4 (0.995185,0.098017), k2=8 (0.980785,0.195090)
    128 132 136 140
    129 133 137 141
    130 134 138 142
    131 135 139 143
  k1=5 (0.881921,0.471397), k2=10 (0.555570,0.831470)
    160 164 168 172
    161 165 169 173
    162 166 170 174
    163 167 171 175
  k1=6 (0.956940,0.290285), k2=12 (0.831470,0.555570)
    192 196 200 204
    193 197 201 205
    194 198 202 206
    195 199 203 207
  k1=7 (0.773010,0.634393), k2=14 (0.195090,0.980785)
    224 228 232 236
    225 229 233 237
    226 230 234 238
    227 231 235 239
 call cftmdl(512, 32)
  k1=1 (0.923880,0.382683), k2=2 (0.707107,0.707107)
    128 144 160 176
    129 145 161 177
（中略）
    142 158 174 190
    143 159 175 191

ほんと、この実装はすごいなぁ。

コメント（4件）

読み物

2007/07/08プログラム

コメント（0件）

2007/07/07(土)プラットホーム汎用で、固定長の整数型を使う

プログラムメモ

int型は整数型ですが、int や long や long long などは環境によってサイズが違ったりします。これらの型は、元もとサイズ（バイト長）を固定する目的で作られたものではないからです。

プラットホーム汎用でプログラムを書く際、固定長のデータを扱う時（バイナリデータ列など）は特に注意する必要があります。C99という規格で固定長整数型として int32t などが規定されましたが、すべての環境で使えるわけではありません。また Windows プラットホームならば、__int32 などが利用出来ますが、Windows以外では利用出来ません。

自分は、だいたいの環境でうまく動くマクロを作って、これを使っています。

#if defined(__C99__) || (defined(__GNUC__) && __GNUC__ >= 3)
#	include <inttypes.h>
#	include <stdint.h>
#else
#	if defined(__GNUC__)
		typedef short			 int16_t;
		typedef unsigned short		uint16_t;
		typedef int			 int32_t;
		typedef unsigned int		uint32_t;
		typedef long long		 int64_t;
		typedef unsigned long long	uint64_t;
#	elif (_MSC_VER || __BORLANDC__)
		typedef __int16			 int16_t;
		typedef unsigned __int16	uint16_t;
		typedef __int32			 int32_t;
		typedef unsigned __int32	uint32_t;
		typedef __int64			 int64_t;
		typedef unsigned __int64	uint64_t;
#	endif
#endif

コメント（1件）

2007/05/23(水)UDPをTCPにクローンするプログラム

プログラム::TCP/IP

UDPを受信してTCPにクローンするプログラム（複数クライアント対応）。本当はマルチスレッドで書いてあげるべきなのですが、そこまですると（実現する処理に対して）高級するぎる気がしたのでポーリングにしました。

というのも、TCPの通信エラーとかが出たときに（ブロックしないはずの）ソケットへの書き込み（データの送信）がバッファ一杯になるとストールして（止まって）しまうんですよね。もちろんソケットや書き込み時の関数はノンブロッキングに設定しておいても、です。C ではまだ確認していませんが、Perlで書いたときはそうでした。

TCP/IPはただ通信するだけなら簡単ですが、この手のエラー処理を考えはじめると非常に頭を使います。

プログラム

メイン部のみ。適当に main ルーチンを書けば Windows でも Unix系でも動きます。

int	max_connections=100;
char	buffer[0x10000];
int	connection_pool[max_connections];

int accept_client(int listen_sock) {
	int sock, len;
	struct sockaddr_in sin;

	// accept
	len = sizeof(sin);
	sock = accept(listen_sock, (struct sockaddr *)&sin, &len);
	if (sock<0) error_return("client accept error");

	// 接続者情報
	printf("[%02d] Connection from %s\n", sock, inet_ntoa(sin.sin_addr));
	// ソケットの設定
	set_non_blocking(sock);
	return sock;
}

//////////////////////////////////////////////////////////////////////////////
// server main
//////////////////////////////////////////////////////////////////////////////
int udp2tcp_server_main(int udp_sock, int tcp_sock) {
	int i;
	char buf[1024];
	// select用の設定
	fd_set fdbits;
	fd_set fdbits_org;
	FD_ZERO(&fdbits_org);
	FD_SET(udp_sock, &fdbits_org);
	FD_SET(tcp_sock, &fdbits_org);

	while(1) {
		// Select
		memcpy(&fdbits, &fdbits_org, sizeof(fdbits_org));
		select(FD_SETSIZE, &fdbits, NULL, NULL, NULL);

		// New connection from TCP
		if ( FD_ISSET(tcp_sock, &fdbits) ) {
			int newsock = accept_client(tcp_sock);
			if (newsock>0) {
				// search for free connection pool point
				for(i=0; i<max_connections; i++)
					if (!connection_pool[i]) break;

				if (i<max_connections) {	// ソケット番号をセーブ
					FD_SET(newsock, &fdbits_org);
					connection_pool[i] = newsock;
					dbg("[%02d] save to connection_pool[%d]\n", newsock, i);
				} else {		// 接続を切る
					printf("Connections max\n");
					close(newsock);
				}
			}
		}

		// Recieved from UDP
		int size = 0;
		if ( FD_ISSET(udp_sock, &fdbits) ) {
			size = recv(udp_sock, buffer, BUF_SIZE, MSG_DONTWAIT);
			if (size<0) error_exit2("UDP recv error(%d)", size);
			dbg("[%02d] Recieved UDP %d bytes\n", udp_sock, size);
		}

		// TCP接続クライアントにデータを送信
		for(i=0; i<max_connections; i++) {
			int sock = connection_pool[i];
			if (!sock) continue;
			// socket からデータ受信
			if ( FD_ISSET(sock, &fdbits) ) {
				int s = recv(sock, buf, 1024, MSG_DONTWAIT);
				if (s<0) {
					FD_CLR(sock, &fdbits_org);
					connection_pool[i] = 0;
					dbg("[%02d] clear to connection_pool[%d]\n", sock, i);
					printf("[%02d] Connection close\n", sock);
					close(sock);
					continue;
				}
			}
			// データ送信
			if (size>0) {
				//int s = write(sock, buf, size);
				int s = send(sock, buf, size, MSG_DONTWAIT);
				dbg("[%02d] write TCP %d bytes\n", sock, s);
			}
		}
	}
}

コメント（3件）

2007/02/04(日)Windows Media Service(WMS) 代替サーバ wmrelay

プログラム::配布プログラム::TCP/IP

wmrelayって？

Windows Media エンコーダ(WME9)には小規模なライブ中継機能（ストリーミングサーバ機能）を持っていますが、最大接続数に制限があり*1、中継元のアップ回線が細いといろいろと不都合が起こります。

Windows Media Services(WMS) などを用いて大規模な中継を行ったり、icecastのように複数の協力者で大規模中継を試みたいところですが、WMS は Windows Serverにしか付属せず、高い上に Windows を外部向けサーバとして運用させるのも気乗りしません。

wmrelay を利用すれば、Linux、FreeBSD、通常のWindowsマシンなどで誰でも手軽に中継サーバを実現できます（mmsプロトコルは非対応です）。協力者を集めれば大規模中継も簡単に行えます。

ダウンロード

ライセンス：GPL (Version 2 or later)

wmrelay.zip Version 2.13
- (Windows) wmrelay.exe を実行（ファイル単体でok）
- (Linux/FreeBSD) wmrelay.pl を Perl 5.8.1 以降で実行

変更履歴

Ver2.13 2010/04/10 -b オプションの追加。バッファディフォルトを 64 に変更。
Ver2.12 2009/08/24 [Windows] EXE化をActivePerl5.10ベースにしました。
Ver2.11 2007/03/08 [Windows] GUIで動作しない不具合を修正しました
Ver2.10 2007/02/13 [Windows] GUI操作画面を作成
Ver2.01 2007/02/05 Mac等で動かない問題を修正（Thanks to RXさん）
Ver2.00 2007/02/03 全面書き直し。pushサーバ機能追加
Ver1.00 2005/06/07 初公開

Linux/FreeBSD等での動作

Linux, FreeBSDなどではPerl 5.8.1以降でithreadが有効なものが必要です。.plをWindows環境(GUI)で使用する場合は、GUI-Loftのインストールが必要です（参考サイト）。

「perl -V」として実行したとき「useithreads=define」と出力されればOKです。

~$ perl -V
Summary of my perl5 (revision 5 version 8 subversion 8) configuration:
＜中略＞
hint=recommended, useposix=true, d_sigaction=define
usethreads=define use5005threads=undef useithreads=define usemultiplicity=define

動作確認

エンコードソフト：Windows Media Encoder 9 (WME9)
再生ソフト：Windows Media Player 6.4*2/9.x, Winamp, MPlayer(Linux)
wmrelay動作環境：FreeBSD 6.2、Ubuntu 6.10(Linux 2.6)、Windows 2000

音声のみ、動画再生共に確認してます。

使い方

wmrelayはコマンドラインから起動します。使用例は次のようになります。

配信元 http://192.168.1.50:8888/
送信ポート標準値（8080）
```
wmrelay http://192.168.1.50:8888/
```
配信元 http://test-stream.dyndns.org:8080/

送信ポート 8000

wmrelay -p 8000 http://test-stream.dyndns.org:8080/

配信元 http://test-stream.dyndns.org:8080/path/file.asf
送信ポート 9000

最大接続数 200

wmrelay.pl -p 9000 -m 200 http://test-stream.dyndns.org:8080/path/file.asf

オプションの解説

-h	コマンドラインヘルプを表示します
-c	CUI(GUIをオフにして)使用します（Windowsでの動作時のみ）
-p [num]	送信に利用するポート番号を指定します（省略時:8080）
-m [num]	最大接続数を指定します（省略時:100）
-b [num]	バッファサイズを指定。16,32,64,128,256,512,1024 どれか（省略時:64）
-l [file]	ログファイル名を指定します。%p はポート番号で置き換わります
-f [file]	プッシュサーバモードでのパスワードファイルを指定します（後述）
-a [pass]	プッシュサーバモードでのパスワードを指定します（後述）
-s	リレー元からのパケット受信情報を表示しません
-d	デバッグ情報を表示します

最大接続数の目安

回線速度（アップ速度）をストリームのビットレートで割れば目安がわかります。スピードテストサイトなどで「Upload Speed」を計測してみてください*3。

「アップ速度：0.5Mbps　ストリーム：136kbps」の場合0.5*1024 / 136 = 3.7647... (connections max)端数は切り捨てて（必ず切り捨てること）3人が理論限界です。

「アップ速度：12.58Mbps　ストリーム：226kbps」の場合12.63*1024 / 226 = 57.226... (connections max) 57人が理論限界ですが50人程度にしておくほうが無難です。

push配信モード

ストリーム受信元URLを指定しないと、wmrelayはpushサーバモードで動作します。このモードでは中継クライアント（Windows Media Encoder9）から中継要求（サーバへのプッシュ要求）があった際に中継が開始されます。

中継先ポートには配信に利用するポート(-p port_number)と同じ番号を指定します。公開ポイントは適当で構いません（wmrelayは公開ポイントを無視します）。

また１つのプログラム（１つのポート）につき同時に１つのストリームのみが中継可能です。パスワードを指定しないと誰でも中継可能となりますのでご注意ください。サーバを止めるときは CTRL-C を何度か押してください。

共通パスワードモード

ユーザー名に関わらず、特定のパスワードを入力したクライアントの接続をすべて受け付けるモードです。パスワードは起動時に次のように指定します。

wmrelay -a password

ユーザー別パスワードモード

起動時にパスワードファイルを指定することで、ユーザー別のパスワードを設定します。

wmrelay -f password_file

パスワードファイルは「ユーザー名=パスワード」という書式になります。

# この行はコメントです
user1=pass1
user2=himitsu

パスワードの設定されているユーザーが居ない場合や、パスワードファイルがない場合には、誰でも中継可能になりますので注意してください。

クライアント（リスナー）の使いかた

例えば、10.11.22.33 というIPのマシンにおいて、標準状態の 8080 で中継した場合は、クライアントは次のようにアクセスします。

http://10.11.22.33:8080/

ブラウザから .asx ファイルを経由して表示させたい場合は、ブラウザで次のアドレスにアクセスするよう指示します。

http://10.11.22.33:8080/relay.asx

参考資料

コメント（13件）

		2007/08
日	月	火	水	木	金	土
			1	2	3	4
5	6	7	8	9	10	11
12	13	14	15	16	17	18
19	20	21	22	23	24	25
26	27	28	29	30	31