仮眠プログラマーのつぶやき

自分がプログラムやっていて、思いついたことをつぶやいていきます。

自作ゲームやツール、ソースなどを公開しております。
①ポンコツ自動車シュライシュラー
DOWNLOAD
②流体力学ソース付き
汚いほうDOWNLOAD
綺麗なほうDOWNLOAD
③ミニスタヲズ
DOWNLOAD
④地下鉄でGO
DOWNLOAD
⑤ババドン
DOWNLOAD
⑥圧縮拳（ツール）
DOWNLOAD
⑦複写拳
DOWNLOAD
⑧布シミュレーション
DOWNLOAD
⑨minecraft巨大電卓地形データ
DOWNLOAD
⑩フリュードランダー
デジゲー博頒布α版
DOWNLOAD
⑪パズドラルート解析GPGPU版
DOWNLOAD
⑫ゲーム「流体de月面着陸」
DOWNLOAD

2013年08月

2013年08月12日02:04
カテゴリ
GPGPU
１億桁×１億桁

円周率プロジェクト途中経過

なんかFFTよりも多倍長計算に向いている高速剰余変換（FMT）というものがあるらしい

で早速やってみた。

回転子であるωを4とし、すべての計算を257の剰余下で行う。
変換する全要素は8つであり、4^8を257で割った余りは1かつ、4^(8/2)を257で割った余りは-1(=256)であることがミソ

以下のソースでは4桁の10進数の数同士を乗算し結果を表示している



FFTを用いた計算よりもFMTが優れている点として、全て整数値同士の計算で完結するということがあげられる。
これは非常に重要だ。
FFTでは計算したい桁数を上げれば上げるほど、浮動小数点の誤差が蓄積してしまう。
doubleの精度を持ってしてでも1億桁×1億桁をやろうとすれば、double１要素につき4桁が限界だ。
これはdoubleつまり8byteに、10進数4桁≒1.66byteしか格納できなく、記憶するメモリ容量的にも無駄が多い。当然桁数あたりの必要な演算量も多くなってくる。
FMTでは4byte整数に約1.5～2byteは敷き詰められるしdoubleと違って演算スピードも早い。

誤差の話に戻るが、結果的にFFTを用いた多倍長の計算では1兆桁付近で精度の限界が生じてしまう。一方FMTは全て整数で計算するので理論上無限大まで誤差の蓄積なく計算することが可能だ。


というわけでFMT試作verのソースをUPする。
正直醜いスパゲッティなのであまり聞かないで・・・・
2014/2/3　YSRさんのご好意によりきれいなソースとなりました！ありがとうございます







#const beki 3
#const LIST_SIZE 1 << beki
#const LIST_SIZE2 LIST_SIZE >> 1
#const LIST_SIZE3 LIST_SIZE << 1
#const LIST_SIZE8 LIST_SIZE << 2
#const n LIST_SIZE
#const p 257
#const pn p/n
#const omega 4

#module
#deffunc fft array A, int inv
if(inv != 0){
dim c, LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =  cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic =  cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic = (iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a,c,LIST_SIZE8@
}
repeat beki@ ;fft
bekii = 1 << cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, t2*LIST_SIZE2@ / bekii, p@)
r4     = A(t1) * w \ p@
A(t1)  = (A(t0) + p@ - r4) \ p@
A(t0) += r4
A(t0)\ = p@
loop
loop
return

#deffunc ufft array A, int inv
repeat beki@ ;fft
bekii = LIST_SIZE2@ >> cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, t2 * LIST_SIZE2@ / bekii, p@)
r4     = A(t1)
A(t1)  = (A(t0) - r4 + p@) * w \ p@
A(t0) += r4
A(t0) \= p@
loop
loop
if(inv != 0){
dim c, LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic = cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic=(iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a, c, LIST_SIZE8@
}
return

#deffunc ifft array A,int inv
if(inv != 0){
dim c, LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic = cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic=(iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a, c, LIST_SIZE8@
}
repeat beki@ ;fft
bekii = 1 << cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, LIST_SIZE@ - t2 * LIST_SIZE2@ / bekii, p@)
r4     = A(t1) * w \ p@
A(t1)  = (A(t0) + p@ - r4) \ p@
A(t0) += r4
A(t0) \= p@
loop
loop
return

#deffunc iufft array A,int inv
repeat beki@;fft
bekii=LIST_SIZE2@>>cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, LIST_SIZE@ - t2 * LIST_SIZE2@ / bekii, p@)
r4     = A(t1)
A(t1) =(A(t0) -r4 + p@) * w \ p@
A(t0) += r4
A(t0) \= p@
loop
loop
if(inv != 0){
dim c,LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic = cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic=(iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a, c, LIST_SIZE8@
}
return

#deffunc narasi array A,int inv
repeat LIST_SIZE@
arhgea  = a.cnt \ n@
a.cnt  /= n@
if(arhgea != 0) :a.cnt += (n@ - arhgea) * pn@ + 1
loop
return

#defcfunc jouyo int a, int inv, int cc
ninv=inv
m = 1.0 * cc
x = (1.0 * a) \ m
jouyoo = 1.0
repeat
/* n&1 はnが奇数(bitの1桁目が1)の時1、偶数(bitの1桁目が0)の時0を返す */
if(int(ninv) \ 2 == 1) :jouyoo = (jouyoo * x) \ m
ninv = ninv / 2
x = (x * x) \ m
if(ninv == 0) :break
loop
if(jouyoo < 0) :jouyoo += m
return int(jouyoo)
#global

dim aa, n
aa = 6,5,9,3,0,0,0,0
dim bb, n
bb = 7,8,6,5,0,0,0,0
dim cccc, n
dim dddd, n
memcpy cccc, aa, n*4
memcpy dddd, bb, n*4
;aa.1 *= -2
;aa.2 *= 4

ufft bb
ufft aa
repeat 8
bb.cnt *= aa.cnt
bb.cnt \= p
loop

ifft bb
narasi bb

repeat 7
bb.(cnt+1) += bb.cnt / 10
bb.cnt     \= 10
loop

pos 290, 0

repeat 4
mes cccc.cnt
pos ginfo_cx - 12, ginfo_cy - 18
loop

mes "*"
pos ginfo_cx-12, ginfo_cy-18

repeat 4
mes dddd.cnt
pos ginfo_cx - 12,ginfo_cy - 18
loop

mes "="
pos ginfo_cx - 12, ginfo_cy - 18

repeat 8
mes bb.cnt
pos ginfo_cx - 12, ginfo_cy - 18
loop

pos 0, 40
mes jouyo(400, 1024 * 64, 164 * 1024 * 128 + 1)
stop









ここで、FFTに時間間引きや周波数間引きの理論を取り入れることで、変換→逆変換の際2回あるビット逆順の計算を省くことが可能となる。
これも将来的にOpenCLへ落としこむ時、ランダムアクセスを減らしてパフォーマンス低下を防ぐ大きな意味を持つ。

さらにRADEON HD 7900シリーズ・・いわゆるGCNと言われるタイプのGPUコアは整数演算が非常に早い！
これでさらに面白くなってきそうだ！

もともとdouble精度の計算ができないGPUでも、mad24の機能で高速に整数演算ができる物が多いため、そういった点でもFFTよりFMTが有利だ！


というわけで次なるステップは、HSPCLへの移植だな


2019/10追記
FMTをPythonとCUDAで実装しました。ガウスルジャンドルで3億桁まで求めた記事をQiitaに公開しました。(この量だとlivedoorブログだと文字制限で圧倒的にオーバーする・・・)
https://qiita.com/Red_Black_GPGPU/items/e933e0d846b874b86c32
ソースはこちら
https://github.com/toropippi/FMT_QiitaSample
  • toropippitoropippi
  • コメント( 2 )
2013年08月11日04:20
カテゴリ
GPGPU
HSP

HSPCLでリアルタイムマンデルブロー集合描画

HSPCL ver1.11で大きな変更点はないが、カーネル内でグローバルアイテムを２次元方向に拡張できるようにした。

ここで面白くなってくるのがマンデルブロー集合の描画
ダウンロード 
無題


普通リアルタイムに描画できないくらい、何度も計算を繰り返さなければいけないのだが、その計算が案外単純で、あまり隣接要素へ相互間作用を及ぼさないアルゴリズムなのでOpenCL移植によって相当な高速化ができた。
core i7 3820とHD7970でfpsに約10～20倍の差がでた。


操作はgoogleマップと同じといえば分かるだろう。マウスホイールのない方は残念ですが・・・はい、拡大縮小出来ません・・・


計算精度がfloatのためあまり拡大すると粗が目立つが、ここはdouble-doubleの四倍精度計算のアルゴリズムを改良してfloat→float-floatの擬似二倍精度でやればほぼdoubleと同じ精度が得られると思われる。

逆にミドルレンジ以下GPUはまずdouble計算できない事が多いので、float-floatがベストだろう




とりあえずソースと実行ファイルをUP。
また、「砂遊び 」のサンプルを快くテスト実行してくださった通りすがりの皆様のお陰で、だいたいHSPCL.dllの動作環境がわかってきた。

恐らくこうであろう実行環境・・・・・
OS ： Windows XP以降
グラボ ： Geforce系なら9800GT以降、RADEONは不明、intel HD Graphics 2500/4000以降
CPU ： intelでかつcore iシリーズなら対応グラボなしでも「Intel SDK for OpenCL」をインストールで動作
Direct X ： 9.0c以降（？openclとは関係ないかも）


だが依然としてHSPSHADより要求スペックが高いのには変わりなかった
  • toropippitoropippi
  • コメント( 0 )
  • トラックバック( 0 )
2013年08月03日14:29
カテゴリ
HSP
アルゴリズム

「FFTモジュール活用例２巨大整数の乗算」のソース

以下のソースはHSPコンテストに出した作品「FFTモジュール活用例２巨大整数の乗算」のソースです。
下の方にあるfftモジュールは著作権フリーなので自由に使って下さい



「fft」の使い方ですが、まずこの命令を使うには変数
「beki」
「LIST_SIZE」
「LIST_SIZE2」
「LIST_SIZE8」
の作成をしておいて下さい。

beki には高速フーリエ変換するサンプル数(2^n)のnを入れます。以下のソースで言うとサンプル数を8192個作りたかったので
beki=13
になっています。

LIST_SIZEには2^beki を入れて下さい
LIST_SIZE2はLIST_SIZE\2を、LIST_SIZE8にはLIST_SIZE×8を代入して下さい。


高速フーリエ変換したい実数(虚数)配列は当然ddim で作り、配列数はLIST_SIZEを指定して下さい
ここで注意ですがbeki には整数しか指定出来ませんので、変換する配列数は2のべき乗に限られます



命令「fft」の引数は３つで
・実数部(浮動小数点配列)
・虚数部(浮動小数点配列)
・逆変換か否か(1 or -1       -1で逆変換)

です。命令を実行すると実数部、虚数部に変換後の値が入ります。












boxf
color 255,128,128
font "MS ゴシック",27,1
mes "Ａ"
color 255,255,255
pos 25,0
mes "×　　 ＝"
color 128,255,128
pos 60,0
mes "Ｂ"
pos 140,2
button goto "答えを計算",*jikkou
color 255,128,128
boxf 0,30,314,480
color 128,255,128
boxf 314,30,640,480
pos 10,66
sdim k1,32768/2;(2^24)
mesbox k1,290,404,1,32766/2-1
sdim k2,32768/2;(2^24)
pos 332,66
mesbox k2,290,404,1,32766/2-1



ddim d1,32768/4
ddim d2,32768/4
ddim chasr,32768/4
ddim chasi,32768/4
stop










;計算実行
*jikkou
clrobj
color 255,255,255
boxf
beki=13;beki 13では8192配列が作成され1配列あたり4桁入るので最終出力が32768桁までなら計算できることになる
LIST_SIZE=1<<beki
LIST_SIZE2=LIST_SIZE>>1
LIST_SIZE8=LIST_SIZE<<3
strk1=strlen(k1)
strk2=strlen(k2)
hot4=4
repeat 32768/4/2;入力の最大桁は32768桁の半分
ccnt=LIST_SIZE2-1-cnt
sitenk1=strk1-cnt*4-4
if sitenk1<=0:hot4=strk1-cnt*4:sitenk1=0
d1.ccnt=double(strmid(k1,sitenk1,hot4));
if sitenk1<=0:break
loop

hot4=4
repeat 32768/4/2
ccnt=LIST_SIZE2-1-cnt
sitenk1=strk2-cnt*4-4
if sitenk1<=0:hot4=strk2-cnt*4:sitenk1=0
d2.ccnt=double(strmid(k2,sitenk1,hot4));
if sitenk1<=0:break
loop

fft d1,d2,1;フーリエ変換

memcpy chasr,d1,8*LIST_SIZE,0,0
memcpy chasi,d2,8*LIST_SIZE,0,0


ReTm=0.0
imTm=0.0
ReTN=0.0
imTN=0.0
ReTm=chasr
ImTm=chasi
ReTN=chasr
ImTN=chasi
d1=0.5*(Retm*ImTm+ReTN*ImTN)
d2=0.25*( ((ReTN+ImTN)*(ReTN-ImTN)) + ( (ImTm+ReTm)*(ImTm-ReTm) ) )
itjn=LIST_SIZE-1
repeat LIST_SIZE-1,1;フーリエ変換後の各要素同士を掛け算。この部分はこちらのページを参考に作った→http://homepage2.nifty.com/m_kamada/math/fftmul.htm
ReTm=chasr.cnt
ImTm=chasi.cnt
ReTN=chasr.itjn
ImTN=chasi.itjn
d1.cnt=0.5*(Retm*ImTm+ReTN*ImTN)
d2.cnt=0.25*( (ReTN+ImTN)*(ReTN-ImTN) + (ImTm+ReTm)*(ImTm-ReTm) )
itjn-
loop


ddim chasr,1
ddim chasi,1

fft d1,d2,-1;逆変換
sdim k2,4
sdim k1,LIST_SIZE*4+80
dim kekka,LIST_SIZE


w=LIST_SIZE-1;繰り上げ
w1=LIST_SIZE-2
kisuu=10000
repeat LIST_SIZE-1;繰り上げ
d1.w+=0.5
d1.w-=d1.w\1
if d1.w>=kisuu{
kuriaged=int(d1.w/kisuu)
d1.w-=kisuu*kuriaged
d1.w1+=kuriaged
}
kekka.w=int(d1.w+0.5)
w-
w1-
loop
d1.w+=0.5
d1.w-=d1.w\1
kekka.w=int(d1.w+0.5)


onf=0
w=0
w1=0
sdim k2,12;文字列に出力
repeat LIST_SIZE-1
if onf=0{
if kekka.cnt!0{
onf=1
k2=str(kekka.cnt)
w1=strlen(k2)
k2=strf("%0"+w1+"d",kekka.cnt)
}
}else{
w1=4
k2=strf("%04d",kekka.cnt)
}
memcpy k1,k2,w1,w,0
w+=w1
loop

color 0,0,0
pos 0,0:mes "答え＝"
mesbox k1,630,430,0,LIST_SIZE*4+80

wait 5
sdim kekka,1;メモリ解放
sdim k2,1
sdim d1,1
sdim d2,1
sdim chasr,1
sdim chasi,1
stop







#deffunc fft array A,array B,int inv
ddim c,LIST_SIZE
ddim d,LIST_SIZE

repeat beki;fft
bekii=LIST_SIZE>>(cnt+1)
repeat LIST_SIZE2
t2 = cnt\bekii;
t0 = (cnt/bekii)*bekii*2 + t2;
t1 = t0+bekii;
r4=A(t1)
i4=B(t1);
rm0=A(t0)-r4
im0=B(t0)-i4;
pii=3.14159265358979323846264338328*t2/bekii*inv;
dsin=sin(pii)
dcos=cos(pii);
A(t0) += r4;
B(t0) += i4;
A(t1) = rm0*dcos-im0*dsin;
B(t1) = rm0*dsin+im0*dcos;
loop
loop


repeat LIST_SIZE;ビット逆順
ic=cnt/65536
ic = (ic&0x00005555)<<1 | (ic&0x0000AAAA)>>1;
ic = (ic&0x00003333)<<2 | (ic&0x0000CCCC)>>2;
ic = (ic&0x00000F0F)<<4 | (ic&0x0000F0F0)>>4;
ic = (ic&0x000000FF)<<8 | (ic&0x0000FF00)>>8;
iic=cnt\65536
iic = (iic&0x00005555)<<1 | (iic&0x0000AAAA)>>1;
iic = (iic&0x00003333)<<2 | (iic&0x0000CCCC)>>2;
iic = (iic&0x00000F0F)<<4 | (iic&0x0000F0F0)>>4;
iic = (iic&0x000000FF)<<8 | (iic&0x0000FF00)>>8;
ic=(iic*32768+ic/2)>>(31-beki);
C(cnt)=A(ic)
D(cnt)=B(ic)
loop

if inv=-1{
repeat LIST_SIZE
a.cnt=c.cnt/LIST_SIZE
b.cnt=d.cnt/LIST_SIZE
loop
}else{
memcpy a,c,LIST_SIZE8
memcpy b,d,LIST_SIZE8
}
return
  • toropippitoropippi
  • コメント( 0 )
  • トラックバック( 0 )
2013年08月01日02:21
カテゴリ
HSP
アルゴリズム

FFTモジュール活用例１スペクトラムアナライザーのソース

以下のソースはHSPコンテストに出した作品「FFTモジュール活用例１スペクトラムアナライザー」のソースです。
下の方にあるfftモジュールは著作権フリーなので自由に使って下さい



「fft」の使い方ですが、まずこの命令を使うには変数
「beki」
「LIST_SIZE」
「LIST_SIZE2」
「LIST_SIZE8」
の作成をしておいて下さい。

beki には高速フーリエ変換するサンプル数(2^n)のnを入れます。以下のソースで言うとサンプル数を32768個作りたかったので
beki=15
になっています。

LIST_SIZEには2^beki を入れて下さい
LIST_SIZE2はLIST_SIZE\2を、LIST_SIZE8にはLIST_SIZE×8を代入して下さい。


高速フーリエ変換したい実数(虚数)配列は当然ddim で作り、配列数はLIST_SIZEを指定して下さい
ここで注意ですがbeki には整数しか指定出来ませんので、変換する配列数は2のべき乗に限られます



命令「fft」の引数は３つで
・実数部(浮動小数点配列)
・虚数部(浮動小数点配列)
・逆変換か否か(1 or -1       -1で逆変換)

です。命令を実行すると実数部、虚数部に変換後の値が入っています。


#const beki 15
#const LIST_SIZE 1<<beki
#const LIST_SIZE2 LIST_SIZE>>1
#const LIST_SIZE3 LIST_SIZE<<1
#const LIST_SIZE8 LIST_SIZE<<3
#const  lstdim LIST_SIZE3/4+11
 _MSGothic = "ＭＳ ゴシック"
font _MSGothic,80,1
mes "※音量注意！"
font _MSGothic,30,1
mes "本ソフトは音声再生機能があります。\n注意して下さい"
wait 100

ddim onryo,5
randomize
onryo=570.0,50.0,2.0,0.012,0.00007
pos 140
datamode=0;0でロードしたデータ、1ブラウン2ピンク3ホワイト4ブルー5パアプル
gosub*syokika

sdim pinkname,64,5
pinkname="ブラウンノイズ","ピンクノイズ","ホワイトノイズ","ブルーノイズ","パープルノイズ"

objsize 140,40
pos 40,0
button gosub "波形読込",*load
pos 180,0
button gosub "波形連続再生or停止",*saisei
pos 320,0
button gosub "カラードノイズ試験生成",*colorz

objsize 52,24
pos 165,213
button gosub "wav解析",*fft0
pos 445,257
button gosub "wav復元",*fft1


sdim aaaaa,300
aaaaa="x1\nx10\nx100\nx10^3\nx10^4\nx10^5\nx10^6\nx10^7"
yokojiku1=0
tatejiku1=0
objsize 40,24;縦軸
pos 0,160
combox tatejiku1,,aaaaa
pos 55,193;横軸
combox yokojiku1,,aaaaa


repeat -1;メインルーチン
mbeki=beki0
await 15
if mbeki!beki0:gosub*syokika
redraw 0
gosub*byouga
redraw 1
loop






*syokika
ddim schasr,LIST_SIZE
ddim schasi,LIST_SIZE
ddim schasr2,LIST_SIZE
ddim schasi2,LIST_SIZE
ddim itjr,LIST_SIZE
ddim itji,LIST_SIZE
ddim hakei,LIST_SIZE
ddim bunkai,LIST_SIZE
if datamode=5{;読み込みデータの解析
if sampkz<LIST_SIZE{
dialog "波形データが短いです"
}else{
repeat LIST_SIZE
dataitj=wpeek(wavedata,sampt*cnt+wbit*hidari)\(1<<(wbit*8))
if wbit=2:if dataitj>32767:dataitj=dataitj-65536
if wbit=1:dataitj-=128:dataitj*=256
schasr.cnt=double(dataitj)
loop
hakeistk hakei,schasr,schasi
}
}else{
gosub*rando
mcis1
fft itjr,itji,1
repeat LIST_SIZE;周波数ごとにいじる
scare=onryo.datamode*expf(logf(sqrt(1.7+cnt))*(datamode-2))
itjr.cnt*=scare
itji.cnt*=scare
loop

hakeistk bunkai,itjr,itji
mcsi2
fft itjr,itji,-1
mcsi1
hakeistk hakei,schasr,schasi
}
return


*load
dialog "wav",16,"WAVEファイル"
exist refstr
dim wavedata,strsize/4+1
bload refstr,wavedata
wbit=wavedata.4/8
streo=wavedata.5/65536;1でモノラル2でステレオ
sampt=wbit*streo;サンプリング単位(byte)
sampkz=(strsize-44)/sampt
hidari=0
if streo=2:dialog "音声はステレオです\n左の音声を使用しますか？\n(いいえ選択で右)",2:if stat=7:hidari=1
datamode=5;0でロードしたデータ、1ブラウン2ピンク3ホワイト4ブルー5パアプル
gosub*syokika
return


*rando
gk=0.0
repeat LIST_SIZE;ランダムノイズ作る
r=double(rnd(32768)-16384)
schasr.cnt=r
gk+=r
loop
gk/=LIST_SIZE
repeat LIST_SIZE
schasr.cnt-=gk
loop
return

*colorz
datamode=rnd(5)
gosub*syokika
dialog pinkname.datamode+"を生成しました"
return
*saisei
if playon{
mmstop:playon-
}else{
dim wf,lstdim
wf=1179011410,36+LIST_SIZE3,1163280727,544501094,16,65537,44100,88200,1048578,1635017060,(LIST_SIZE3-82)

repeat LIST_SIZE
wpoke wf,44+cnt*2,limit((hakei.cnt+32768)\65536,0,65535)
loop

memfile wf ; ストリームの直前で指定してください
mmload "MEM:a.wav",0,1 ; 画像の拡張子識別のためダミー名a.jpgを使用
sdim wf,64 ; もったいないのでbufのメモリ領域を小さくしておきます
mmplay 0:playon+
}
return

*byouga
pget -1:boxf
font _MSGothic,100
color 128,128,128
pos 140,193
mes "▼"
pos 420,200
mes "▲"
font _MSGothic,10
pos 8,183:mes "縦軸\n　　／横軸"
tt="1"
repeat 5
pos 132.87712*cnt+25,450
mes tt+"Hz"
tt+="0"
loop


;グラフ形
repeat 2
color 200,200,200:boxf 40,47+255*cnt,640,191+255*cnt
loop



;波形
color ,165
yokoj600=LIST_SIZE>>yokojiku1
tatej600=expf(logf(10)*tatejiku1)
repeat 600
limi=limit(180-hakei(yokoj600*cnt/600)*tatej600*0.001983642578,50,180)
if cnt=0:pos cnt+40,limi
line cnt+40,limi
loop
;波形2
tatej600=1
bxs=40
clf=0
getkey key,1
repeat 600
mtatej600=tatej600
tatej600=int(expf(logf(LIST_SIZE-2)*(1.0*cnt/599.0)))
if mtatej600!tatej600{
limi=446-limit(logf(bunkai.tatej600)*12.5076816-144,0,144)
itjint=cnt+41
color clf*80,165
if abs(mousey-375)<73:if bxs<=mousex:if mousex<itjint:color :if key{
repeat tatej600-mtatej600,mtatej600+1
itjdl=expf(((446.0-mousey)*0.07995087+11.512925))
nanbai=itjdl/bunkai.cnt
schasr2.cnt*=nanbai
schasi2.cnt*=nanbai
bunkai.cnt=itjdl
loop
}
boxf bxs,limi,itjint-1,446
bxs=itjint
clf=1-clf
}
loop
return


*fft0
mcis1
fft itjr,itji,1
mcsi2
hakeistk bunkai,schasr2,schasi2
return

*fft1
mcis2
fft itjr,itji,-1
mcsi1
hakeistk hakei,schasr,schasi
return





;if 0 {
#deffunc fft array A,array B,int inv
ddim c,LIST_SIZE
ddim d,LIST_SIZE

repeat beki;fft
bekii=LIST_SIZE>>(cnt+1)
repeat LIST_SIZE2
t2 = cnt\bekii;
t0 = (cnt/bekii)*bekii*2 + t2;
t1 = t0+bekii;
r4=A(t1)
i4=B(t1);
rm0=A(t0)-r4
im0=B(t0)-i4;
pii=3.14159265358979323846264338328*t2/bekii*inv;
dsin=sin(pii)
dcos=cos(pii);
A(t0) += r4;
B(t0) += i4;
A(t1) = rm0*dcos-im0*dsin;
B(t1) = rm0*dsin+im0*dcos;
loop
loop


repeat LIST_SIZE;ビット逆順
ic=cnt/65536
ic = (ic&0x00005555)<<1 | (ic&0x0000AAAA)>>1;
ic = (ic&0x00003333)<<2 | (ic&0x0000CCCC)>>2;
ic = (ic&0x00000F0F)<<4 | (ic&0x0000F0F0)>>4;
ic = (ic&0x000000FF)<<8 | (ic&0x0000FF00)>>8;
iic=cnt\65536
iic = (iic&0x00005555)<<1 | (iic&0x0000AAAA)>>1;
iic = (iic&0x00003333)<<2 | (iic&0x0000CCCC)>>2;
iic = (iic&0x00000F0F)<<4 | (iic&0x0000F0F0)>>4;
iic = (iic&0x000000FF)<<8 | (iic&0x0000FF00)>>8;
ic=(iic*32768+ic/2)>>(31-beki);
C(cnt)=A(ic)
D(cnt)=B(ic)
loop

if inv=-1{
repeat LIST_SIZE
a.cnt=c.cnt/LIST_SIZE
b.cnt=d.cnt/LIST_SIZE
loop
}else{
memcpy a,c,LIST_SIZE8
memcpy b,d,LIST_SIZE8
}
return



#deffunc hakeistk array A,array B,array E
repeat LIST_SIZE
a.cnt=sqrt(B(cnt)*B(cnt)+E(cnt)*E(cnt))
loop
return


#deffunc mcis1
memcpy itjr,schasr,LIST_SIZE8
memcpy itji,schasi,LIST_SIZE8
return

#deffunc mcis2
memcpy itjr,schasr2,LIST_SIZE8
memcpy itji,schasi2,LIST_SIZE8
return

#deffunc mcsi1
memcpy schasr,itjr,LIST_SIZE8
memcpy schasi,itji,LIST_SIZE8
return

#deffunc mcsi2
memcpy schasr2,itjr,LIST_SIZE8
memcpy schasi2,itji,LIST_SIZE8
return
;}
  • toropippitoropippi
  • コメント( 0 )
  • トラックバック( 0 )
プロフィール

toropippi

記事検索


【お気に入りブログ】
┣アド吉ブログ
┣シープのオフィシャルブログ
┣「ファンの一人だよ。」
┣鋼鉄親子でゲーム漬け！
┣ゲームって何？
┗空のつくり方

アクセスカウンター

    QRコード
    QRコード
    最新記事(画像付)
    • ライブドアブログ

    Template by decoweb