-
- [1]
改訂 十進 BASIC による プログレッシブ JPG の展開と画像化。
- 投稿者：SECOND
- 投稿日：2009年12月14日(月)20時36分38秒
-
-
!改訂 十進 BASIC による プログレッシブ JPG の展開と画像化。
!
!先投稿では、プログレッシブ JPG 再生過程の画像を、何枚も表示していたが、
!途中画像は、最初のＤＣ成分１枚だけ と、最終完成画、全２枚とした。
!Baseline JPG は、全１枚なので、画数で両者を区別できる。
!（必要なら、再生過程 全ての画像も、表示できるよう、SUB IZZRL0 に注釈行がある）
!
!色差成分 Cb Cr の間引き走査復元の塗潰しを、SUB IDDCT8X8 に統合し、簡素化した。
!その他、各所の手入れ多数。
!
!具体的、可視的なプログラムで、実行し画像化しますので、詳細事項の追跡と御参考に。
!再生できるファイルは、1000x1000 までの JPG だけで、
! baseline , spectral selection , successive approximation の３種類( web 上の、ほぼ全種)
! 参考資料：http://www.w3.org/Graphics/JPEG/itu-t81.pdf
!
!１）successive approximation AC.subsequence（Y,Cb,Cr 別々、１bit づつの処理）
!
! 0 1 1 0 0 0 0 0 0 ？
! 0 0 0 0 0 1 0 0 0 ？
! 0 1 0 0 0 0 0 1 0 ？
! 0 1 1 0 0 1 0 1 0 ？
! 0 1 1 0 0 1 0 1 0 ？
! --------------------------------------------------------------------------
! ±1 b1 b2 0 0 b3 0 b4 ±1 ？
! 前の終り RRRR RRRR RRRR extend. 次の始め
!
! huffman.
! RRRRssss extend. b1 b2 b3 b4 …
! 3 1 (0 or 1) bit_stream=？何個になるかは、上図で、上位桁 =0 の係数が
! -1 +1 (0 or 1) RRRR 個 になるまでに通過した上位桁 <>0 の個数。
! 0 ±1
! 0 → 無変化。
! 1 → ±符号は上位桁に合せて加算。(絶対値が＋1)
!
!２）successive approximation DC.subsequence（Y,Cb,Cr 別々、１bit づつの処理）
!
! ハフマン・コード RRRRssss 部は、存在せず、
! 頭からの bit_stream.で、１bit づつ、全てのblock の DC係数 に加える。
!
! AC と同様、0 → 無変化。1 → ±符号は上位桁に合せて加算。(絶対値が＋1)
!
!※上記、successive approximation AC, DC とも、加える１は、
! 2^Al 倍の point transfer. としてから、加えます。
!
DEBUG ON
!------------------------
!JPG.decoder
! Baseline
! Progressive( spectral selection )( successive approximation )
!------------------------
OPTION ARITHMETIC NATIVE
OPTION BASE 0
OPTION CHARACTER byte
SET TEXT background "OPAQUE"
ASK BITMAP SIZE bmx,bmy
SET WINDOW 0,bmx, bmy,0
SET ECHO "OFF"
SET COLOR MODE "NATIVE"
!
DIM D8(1000,1000) !MAIN65
DIM D2(1000,1000,2) !Y=D2(,,0) Cb=D2(,,1) Cr=D2(,,2)
DIM D1(1000,1000,2) !Y=D2(,,0) Cb=D2(,,1) Cr=D2(,,2)
DIM MH(2),MV(2) !R_BIN31 SOF0 MCU.Ybr.H()V()
DIM HDC(2),HAC(2) !R_BIN31 hT.table selection
DIM QS(2),CoID(255) !R_BIN31 qT.table selection
DIM M3(2)
!
DIM U(63),V(63) !zigzag
DIM DQ(7,7,3) !blk8x8 DQT
DIM DH(16,7),DV(255,7) !DHT
DIM B(255+1,7),L(255,7) !encorder & decorder's pre_table, length, ( MAKE_H2,MAKE_H0)
DIM A(2000,7) !decorder
DIM B2(2) !Ybr D.C.成分 starting & back_level for difference
DIM T(7,7),X(7),XO(7) !DDCT8X8, IDDCT8X8
!
LET BST=2 !huffman decorder's bit step 1=8.5s 2=6.5s 4=8.0s 8=50.0s
LET SHb=2^BST !huffman decorder n*SHb=(shl n,BST) n/SHb=(shr n,BST)
!
!---zigzag table
FOR V_=0 TO 7
FOR U_=0 TO 7
READ i
LET U(i)=U_
LET V(i)=V_
NEXT U_
NEXT V_
DATA 0, 1, 5, 6,14,15,27,28
DATA 2, 4, 7,13,16,26,29,42
DATA 3, 8,12,17,25,30,41,43
DATA 9,11,18,24,31,40,44,53
DATA 10,19,23,32,39,45,52,54
DATA 20,22,33,38,46,51,55,60
DATA 21,34,37,47,50,56,59,61
DATA 35,36,48,49,57,58,62,63
!
DO
FILE GETNAME FL$, "jpg"
IF FL$="" THEN
PRINT "入力ファイル名無し"
EXIT DO
END IF
PRINT "入力ファイル："& FL$
!---
CLEAR
CALL IZZRL0 ! D2()<-- decord JPG
PRINT "次のファイル[ 左クリック ]"
beep
DO
MOUSE POLL j,i,mlb,mrb !CHARACTER INPUT CLEAR: w$
WAIT DELAY 0
LOOP UNTIL 0< mlb OR 0< mrb
LOOP UNTIL 0< mrb
PRINT "終了。"
!-------- IZZRL0 call here for display D2()
SUB MAIN65
LET tester=TIME
PRINT "画像の準備中、";
CALL IDDCT8X8 ! D1()<-- iDCT<-- iDQT<-- D2()
!------ JPG 色空間 ----------------------------
! | Y | | 0.2990 +0.5870 +0.1140 | | R |
! |B-Y| = |-0.1687 -0.3313 +0.5000 | | G |
! |R-Y| | 0.5000 -0.4187 -0.0813 | | B |
!
! | R | | 1 0 +1.40200 | | Y |
! | G | = | 1 -0.34414 -0.71414 | |B-Y|
! | B | | 1 +1.77200 0 | |R-Y|
!----------------------------------------------
FOR V0=0 TO DY-1
FOR U0=0 TO DX-1
!--- RGB<-- Ybr
LET w1=IP(D1(U0,V0,0) +1.40200*D1(U0,V0,2)) !R
LET w2=IP(D1(U0,V0,0) -0.34414*D1(U0,V0,1) -0.71414*D1(U0,V0,2)) !G
LET w3=IP(D1(U0,V0,0) +1.77200*D1(U0,V0,1)) !B
IF w1< 0 THEN
LET w1=0
ELSEIF 255< w1 THEN
LET w1=255
END IF
IF w2< 0 THEN
LET w2=0
ELSEIF 255< w2 THEN
LET w2=255
END IF
IF w3< 0 THEN
LET w3=0
ELSEIF 255< w3 THEN
LET w3=255
END IF
LET D8(U0,V0)=w3*65536+w2*256+w1 !（逆）BGR
NEXT U0
NEXT V0
PRINT TRUNCATE(TIME-tester,2);"秒"
LET w=TRUNCATE(MIN( (bmx-1)/DX,(bmy-1)/DY),1)
IF 1< w THEN LET w=IP(w)
IF 4< w THEN LET w=4
PRINT "描画の倍率=";w
MAT PLOT CELLS,IN 1,1; DX*w, DY*w :D8
END SUB
!========================
!inverse haffman Transform.
SUB IZZRL0
LET byt=0 !!!
CALL ROPEN ! FL$
!---
CALL R_BIN31(0) !A() B(i,J)L(i,J)<-- DH(), return at img.top
PRINT right$("000"& BSTR$(byt,16),4) !!!
PRINT "(";STR$(DX);"x";STR$(DY);
!---
MAT D8=ZER(DX-1,DY-1) !MAIN65
LET i=8*MH(0) !MCU Y.Hsize
LET j=8*MV(0) !MCU Y.Vsize
LET DUM=CEIL(DX/i)*i !Uwidth=bound by MCU Y.Hsize
LET DVM=CEIL(DY/j)*j !Vwidth=bound by MCU Y.Vsize
MAT D1=ZER(DUM-1,DVM-1,2) !Y=D1(,,0) Cb=D1(,,1) Cr=D1(,,2)
MAT D2=ZER(DUM-1,DVM-1,2) !Y=D2(,,0) Cb=D2(,,1) Cr=D2(,,2)
LET MH_=MH(0)
LET MV_=MV(0)
LET DU =DUM !Uwidth=bound by MCU Y.Hsize
LET DV_=DVM !Vwidth=bound by MCU Y.Vsize
LET DU8=CEIL(DX/8)*8 !Uwidth=bound by block Y.Hsize
LET DV8=CEIL(DY/8)*8 !Vwidth=bound by block Y.Vsize
!---
PRINT "/ ";STR$(DU8);",";STR$(DV8);"/ ";STR$(DUM);",";STR$(DVM);")"
CALL frame
!---
PRINT "M3()=";M3(0);M3(1);M3(2)
CALL MAIN65 ! Baseline.最終、Progressive.1st.
!---
IF 0< M THEN PRINT " (";STR$(DX);"x";STR$(DY);"/";STR$(U0);",";STR$(V0);") ";debug$;" abort by ";BSTR$(M,16) !!!
PRINT right$("000"& BSTR$(byt-2*SGN(M),16),4) !!!
CALL R_BIN31(M) ! return at img.top, or EOI
!---
DO WHILE M=BVAL("DA",16) !SOS
IF 0<=HAC(0) THEN
LET MV(0)=1
LET MH(0)=1
LET DU=DU8
LET DV_=DV8
END IF
CALL frame
LET MV(0)=MV_
LET MH(0)=MH_
LET DU=DUM
LET DV_=DVM
!---
PRINT "M3()=";M3(0);M3(1);M3(2) !文末参照：M30<>99(balance), M30=99(un-balance)
IF M30=0 OR M30=64 THEN CALL MAIN65 !Progressive.最終スキャン後の画像
!IF M30<>99 THEN CALL MAIN65 !Progressive.各スキャン毎、Ybr 揃った画像のみ
!CALL MAIN65 !Progressive.各スキャン毎、全画像
!---
IF 0< M THEN PRINT " (";STR$(DX);"x";STR$(DY);"/";STR$(U0);",";STR$(V0);") ";debug$;" abort by ";BSTR$(M,16) !!!
PRINT right$("000"& BSTR$(byt-2*SGN(M),16),4)
CALL R_BIN31(M) ! return at img.top
LOOP
CLOSE #1 ! FL$
END SUB
SUB reset0
LET B2(0)=0 !ROUND( YDC0/DQ(0,0,QS(0)) ) !prediction YDC.( 1st.reference level)
LET B2(1)=0 !prediction CbDC.
LET B2(2)=0 !prediction CrDC.
LET Hx=0 !bits stream input buffer 0~(7+8)bits, use fraction
LET BC=0 !stored bits in Hx
LET NA=0 !nest adr. in A()
LET EOB=0 !counter( end_of_band)
LET M=0
LET ext=0
END SUB
!Page-2 へ続く
-
- [2]
Page-2
- 投稿者：SECOND
- 投稿日：2009年12月14日(月)20時39分14秒
-
-
!Page-2 の始め
SUB frame
PRINT " Ss Se AhAl: ";Ss_;Se_;STR$(Ah);STR$(Al)
PRINT " Y HDC HAC: ";IP(HDC(0)/2);IP(HAC(0)/2)
PRINT " Cb : ";IP(HDC(1)/2);IP(HAC(1)/2)
PRINT " Cr : ";IP(HDC(2)/2);IP(HAC(2)/2)
CALL reset0
!---
FOR V09=0 TO DV_-1 STEP 8*MV(0)
FOR U09=0 TO DU-1 STEP 8*MH(0)
IF rct=0 THEN
CALL R_BIN31(0) ! read marker
IF rct<>DRI THEN BREAK ! not RST0~7
CALL reset0 ! Restart
END IF
CALL MCUxx11 ! read picture data
LET rct=rct-1
!---
IF 0< ext THEN
IF ext=103001 THEN
PRINT "abort marker ";BSTR$(M,16)
IF BVAL("D0",16)<=M AND M<=BVAL("D7",16) THEN ! RST0~7( restart marker)
LET rct=DRI ! set counter
CALL reset0 ! Restart
ELSE
EXIT SUB ! others marker
END IF
ELSE
PRINT "file error. display fragment"
LET M=BVAL("D9",16) ! EOI
EXIT SUB
END IF
END IF
NEXT U09
NEXT V09
IF 0< EOB THEN PRINT "EOBn over frame";EOB !!!
END SUB
SUB MCUxx11
!---read MCU
FOR P=0 TO CMO
IF 0<=HDC(P) OR 0<=HAC(P) THEN
FOR V0=V09 TO V09+8*MV(P)-1 STEP 8
FOR U0=U09 TO U09+8*MH(P)-1 STEP 8
WHEN EXCEPTION IN
IF EOB=0 THEN CALL R_BLK0 ELSE LET EOB=EOB-1
USE
LET ext=EXTYPE
EXIT SUB
END WHEN
!---extend bitmap
IF 0< Ah AND 0< Se_ THEN
FOR i=A_ TO Se_
IF D2(U0+U(i),V0+V(i),P)<>0 THEN
LET L_=1
WHEN EXCEPTION IN
CALL DEC1_EX
USE
LET ext=EXTYPE
EXIT SUB
END WHEN
LET V_=SGN(D2(U0+U(i),V0+V(i),P))*V_*2^Al
LET D2(U0+U(i),V0+V(i),P)=D2(U0+U(i),V0+V(i),P) +V_
END IF
NEXT i
LET A_=Ss_
END IF
!---
NEXT U0
NEXT V0
END IF
NEXT P
END SUB
!------
SUB R_BLK0
IF Ss_=0 THEN
!---D.C.part
LET debug$="DC.huffman" !!!
IF Ah=0 THEN !-----baseline.progSS.progSA(1st.scan).
LET J=HDC(P) !huffman D.C.table selection P( 0=Y 1=Cb 2=Cr)
CALL DEC1_NS
LET EL=V_ !extent length
!---D.C.extent
LET debug$="DC.huffman extend" !!!
IF 0< EL THEN
LET L_=EL
CALL DEC1_EX !keep EL, V_=extent value( length EL bits)
LET W=2^(EL-1) !minimum in EL bits length
IF V_< W THEN LET V_=V_-(W*2-1) !restore signed value
LET B2(P)=B2(P)+V_*2^Al !point transform, integrate to D.C.
END IF
LET D2(U0+U(0),V0+V(0),P)=B2(P)
ELSE !-----progSA(2st.scan).
LET L_=1
CALL DEC1_EX
LET V_=SGN(D2(U0+U(0),V0+V(0),P))*V_
LET D2(U0+U(0),V0+V(0),P)=D2(U0+U(0),V0+V(0),P) +V_*2^Al
END IF
LET Sa_=1
ELSE
LET Sa_=Ss_
END IF
!---A.C.parts
IF Se_=0 THEN EXIT SUB !band Ss_~Se_
LET J=HAC(P) !huffman A.C.table selection P( 0=Y 1=Cb 2=Cr)
LET debug$="AC.huffman"
FOR A_=Sa_ TO Se_
CALL DEC1_NS
LET EL=MOD(V_,16) !extent length
LET RL= IP(V_/16) !run length
!---
IF RL<=14 AND EL=0 THEN !End Of Block(00). End Of Band(10,20,,E0)
!---EOBn extend
LET debug$="EOBn extend"& STR$(RL)
IF 0< RL THEN
LET L_=RL !extend= run_length
CALL DEC1_EX !keep RL, V_=run value( length RL bits)
LET EOB=V_+2^RL-1 !RL= End Of Band(10,20,,E0) run length
END IF
EXIT SUB
!---
END IF
!---RL=(0~15)EL=(1~10), RL=(15)EL=(0)
LET debug$="AC.huffman extend" !!!
IF Ah=0 THEN !-----baseline.progSS.progSA(1st.scan).
LET A_=A_+RL !skip zero_run_length 0~15
!---A.C.extent
IF 0< EL THEN !ZRL(16) only skip
LET L_=EL
CALL DEC1_EX !keep EL, V_=extent value( length EL bits)
LET w=2^(EL-1) !minimum in EL bits length
IF V_< w THEN LET V_=V_-(w*2-1) !restore signed value
!---
LET V_=V_*2^Al !point transform
LET D2(U0+U(A_),V0+V(A_),P)=V_
END IF
ELSE !-----progSA(2st.scan).
IF 0< EL THEN !ZRL(16) only skip
LET L_=EL
CALL DEC1_EX !keep EL, V_=extent value( length EL bits)
IF EL<>1 THEN PRINT "AC.2nd.=";EL;V_ !!!
LET V01=V_
END IF
FOR i=A_ TO Se_
IF D2(U0+U(i),V0+V(i),P)<>0 THEN !zz(k)=xxx_1?/0?
LET L_=1
CALL DEC1_EX
LET V_=SGN(D2(U0+U(i),V0+V(i),P))*V_
LET D2(U0+U(i),V0+V(i),P)=D2(U0+U(i),V0+V(i),P) +V_*2^Al
ELSEIF RL=0 THEN !zz(k)=000_V01
EXIT FOR
ELSE !zz(k)=000_0 ,zero run
LET RL=RL-1
END IF
NEXT i
IF 0< EL THEN !ZRL(16) skip
IF V01=0 THEN LET V01=-1
LET D2(U0+U(i),V0+V(i),P)=V01*2^Al
END IF
LET A_=i
END IF
NEXT A_
END SUB
!========================
! decorder
! J= huffman code table selection ( 0=YDC 1=YAC 2=CDC 3=CAC)
! V_= pickup RRRRssss <-- JPG.file
SUB DEC1_NS
DO
IF BC< BST THEN CALL DEC1_IN
LET W=IP(Hx) ! bits width BST
!----
LET W=A(NA+W,J)
IF 32768<=W THEN EXIT DO
LET NA=W ! nest adr. W=0 table end
LET BC=BC-BST
LET Hx=MOD(Hx*SHb,SHb)
LOOP
LET NA=0 ! DU0L LLLL VVVV VVVV
LET L_=MOD(IP(W/256),128) ! U0L LLLL
LET V_=MOD(W,256) ! VVVV VVVV
IF 16< L_ THEN PRINT "unused code" !BREAK !unused code ! LET V_=BVAL("8000",16)
!----
LET W=MOD(L_,BST)
IF 0< W THEN
LET BC=BC-W
LET Hx=MOD(Hx*2^W,SHb)
ELSE
LET BC=BC-BST
LET Hx=MOD(Hx*SHb,SHb)
END IF
END SUB
SUB DEC1_IN
CALL RED_D
LET W=ORD(D$)
IF W=255 THEN
CALL RED_D
LET M=ORD(D$)
IF M<>0 THEN LET w=1/0 ! EXTYPE=3001, ffxx marker, abnormally break
END IF
LET Hx=Hx+W*2^(BST-8-BC)
LET BC=BC+8
END SUB
!-------
SUB DEC1_EX
LET V_=0
DO
IF L_< 1 THEN EXIT SUB
IF BC< L_ THEN CALL DEC1_IN
LET W=IP(Hx)
!----
IF BST>=L_ THEN EXIT DO
LET V_=V_*SHb+W
LET L_=L_-BST
LET BC=BC-BST
LET Hx=MOD(Hx*SHb,SHb)
LOOP
LET V_=V_*2^L_+IP(W*2^(L_-BST))
!----
LET BC=BC-L_
LET Hx=MOD(Hx*2^L_,SHb)
END SUB
!Page-3 へ続く
-
- [3]
Page-3
- 投稿者：SECOND
- 投稿日：2009年12月14日(月)20時41分34秒
-
-
!Page-3 の始め
!============
! B(,J)L(,J)<-- DH(,J) for decorder table A(,J)
!
SUB makeH0(J)
LET i=0 ! コード生成 順番(短い順)
LET Hx=0
LET Tx=BVAL("8000",16)
FOR L_=1 TO 16
FOR P=1 TO DH(L_,J)
LET L(i,J)=L_
LET B(i,J)=Hx ! コード(生成順), 座標DV(頻度降順) と同順。
LET i=i+1
LET Hx=Hx+Tx
NEXT P
LET Tx=Tx/2
NEXT L_
LET B(256,J)=0
FOR i=i TO 255
LET L(i,J)=0
LET B(i,J)=0
NEXT i
END SUB
!============
!A(,J)=output decorder table<-- B(,J) L(,J) DH(,J) DV(,J)
!
SUB makeD0(J)
FOR LH=16 TO 1 STEP -1
IF DH(LH,J)<>0 THEN EXIT FOR
NEXT LH !length max. in huffman table
LET LM=CEIL(LH/BST)*BST !length max. bound by BST
!---
LET I=0 !start huffman table adr.
LET LA=0 !line adr.
LET P=BST !start Decord code width
LET U_=2^(16-BST) !start Decord code step
LET NC=0 !next start Decord code
DO
LET D_=NC !start Decord code
LET NC=-1
LET LB=LA+(65536-D_)/U_ !1st nest adr.
DO
CALL SERCH
IF 0< L_ THEN
LET A(LA,J)= BVAL("8000",16)+L_*256+DV(I,J) !b15=end. +L.+V.
ELSEIF P=LM THEN
LET A(LA,J)= BVAL("C000",16)+LH*256 !b15=end. b14=Unused. +L.
ELSE
IF NC=-1 THEN LET NC=D_
LET A(LA,J)=LB ! nest adr.
LET LB=LB+SHb ! next nest adr.
END IF
LET D_=D_+U_
LET LA=LA+1
LOOP UNTIL IP(D_)=65536
LET P=P+BST
LET U_=U_/SHb !shr(U_,BST)
LOOP UNTIL P>LM
!---
FOR LA=LA TO 255
LET A(LA,J)=0 !(0),table stop mark
NEXT LA
END SUB
SUB SERCH
FOR I=I TO DH(0,J)-1
LET L_=L(I,J)
IF L_<=P THEN LET w=IP(D_/2^(16-L_))*2^(16-L_) ELSE EXIT FOR
IF w<=B(I,J) THEN
IF w=B(I,J) THEN EXIT SUB ELSE EXIT FOR
END IF
NEXT I
LET L_=-1
END SUB
!===========
! Inverse Fast Cosin Transform.( 8x8, iDCT-2 ) ← Inverse Quantization.DQ()
SUB IDDCT8X8
FOR P=0 TO CMO ! (0=Y,1=Cb,2=Cr)
FOR V0=0 TO DV_-1 STEP 8*MV(0)/MV(P)
FOR U0=0 TO DU-1 STEP 8*MH(0)/MH(P)
!----
FOR V_=0 TO 7
FOR U_=0 TO 7
LET X(U_)=D2(U0+U_,V0+V_,P) *DQ(U_,V_,QS(P)) ! Inverse Quantization
NEXT U_
CALL IWANG ! inverse DCT horizontal_row_8
FOR X_=0 TO 7
LET T(X_,V_)=X(X_)
NEXT X_
NEXT V_
!----
FOR X_=0 TO 7
FOR V_=0 TO 7
LET X(V_)=T(X_,V_)
NEXT V_
CALL IWANG ! inverse DCT vertical_column_8
LET i=X_*MH(0) !expand pt.X
FOR Y_=0 TO 7
IF P=0 THEN
LET D1(U0+X_,V0+Y_,P)=X(Y_)+128 ! inverse level shift
ELSE
LET j=Y_*MV(0) !expand pt.Y
!----expand
FOR V_=j TO j+MV(0)-1
FOR U_=i TO i+MH(0)-1
LET D1(U0+U_,V0+V_,P)=X(Y_) ! CbCr to MCU scale
NEXT U_
NEXT V_
!----expand
END IF
NEXT Y_
NEXT X_
!----
NEXT U0
NEXT V0
NEXT P
END SUB
!----inverse Wang.( 8, iDCT-2 )
SUB IWANG
LET XO(0)=SQR(2/8)*X(0)
LET XO(1)=SQR(2/8)*X(4)
LET XO(2)=SQR(2/8)*X(2)
LET XO(3)=SQR(2/8)*X(6)
LET XO(4)=SQR(1/8)*X(1)
LET XO(5)=SQR(1/8)*X(5)
LET XO(6)=SQR(1/8)*X(3)
LET XO(7)=SQR(1/8)*X(7)
!
LET X(4)=(COS(PI /16)*XO(4)+SIN(PI /16)*XO(7))
LET X(5)=(COS(PI*5/16)*XO(5)+SIN(PI*5/16)*XO(6))
LET X(6)=(SIN(PI*5/16)*XO(5)-COS(PI*5/16)*XO(6))
LET X(7)=(SIN(PI /16)*XO(4)-COS(PI /16)*XO(7))
!
LET XO(4)= X(4)+X(5)
LET XO(5)= X(4)-X(5)
LET XO(6)=-X(6)+X(7)
LET XO(7)= X(6)+X(7)
!
LET X(0)=(COS(PI/4)*XO(0)+COS(PI /4)*XO(1))
LET X(1)=(COS(PI/4)*XO(0)-COS(PI /4)*XO(1))
LET X(2)=(SIN(PI/8)*XO(2)-SIN(PI*3/8)*XO(3))
LET X(3)=(COS(PI/8)*XO(2)+COS(PI*3/8)*XO(3))
LET X(4)=XO(4)
LET X(5)=XO(6)
LET X(6)=XO(5)
LET X(7)=XO(7)
!
LET XO(0)=X(0)+X(3)
LET XO(1)=X(1)+X(2)
LET XO(2)=X(1)-X(2)
LET XO(3)=X(0)-X(3)
LET XO(4)=X(7)*SQR(2)
LET XO(5)=X(6)-X(5)
LET XO(6)=X(6)+X(5)
LET XO(7)=X(4)*SQR(2)
!
LET X(0)=XO(0)+XO(7)
LET X(1)=XO(1)+XO(6)
LET X(2)=XO(2)+XO(5)
LET X(3)=XO(3)+XO(4)
LET X(4)=XO(3)-XO(4)
LET X(5)=XO(2)-XO(5)
LET X(6)=XO(1)-XO(6)
LET X(7)=XO(0)-XO(7)
END SUB
!=============
SUB R_BIN31(M) ! decord(M) before new.search(M)
DO
IF M=BVAL("D8",16) THEN ! SOI
MAT DH=ZER ! clear Huffman Table
LET DRI=0 ! clear Restart Interval.value for RST0~7(restart marker)
LET rct=-1 ! Interval.counter, valid (0<=rct), invalid (rct< 0)
MAT M3=ZER ! clear scan band sum
ELSEIF M=BVAL("D9",16) THEN ! EOI
EXIT DO ! close & end_sub
ELSEIF BVAL("D0",16)<=M AND M<=BVAL("D7",16) THEN ! RST0~7( restart marker)
LET rct=DRI ! set counter with Restart Interval
EXIT SUB
ELSEIF 0< M THEN !M=0 is data"FF" in picture area
CALL RED_D
LET N=ORD(D$)*256
CALL RED_D
LET N=N+ORD(D$)-2 ! N=remain size
!---
IF BVAL("E0",16)<=M AND M<=BVAL("EF",16) THEN ! APP0~APP15
CALL FFE0
ELSEIF M=BVAL("DD",16) THEN
CALL FFDD ! DRI load DRI & rct=DRI
ELSEIF M=BVAL("FE",16) THEN
CALL FFFE ! COMMENT
ELSEIF M=BVAL("C4",16) THEN
CALL FFC4 ! DHT
ELSEIF M=BVAL("DB",16) THEN
CALL FFDB ! DQT
ELSEIF M=BVAL("C0",16) OR M=BVAL("C2",16) THEN
PRINT right$("000"& BSTR$(byt-4,16),4)
!---
CALL FFC0 ! SOF0 SOF2
!---
PRINT " SOF";STR$(MOD(M,16));" MCU_HV Ybr ";STR$(MH(0));STR$(MV(0));
PRINT " ";STR$(MH(1));STR$(MV(1));" ";STR$(MH(2));STR$(MV(2))
ELSEIF M=BVAL("DA",16) THEN
CALL FFDA ! SOS
EXIT SUB ! without close
ELSE
BREAK ! new marker
END IF
END IF
!---
DO
LET M=BVAL("D9",16) ! EOI, 256 ! end of file
CHARACTER INPUT #1,IF MISSING THEN EXIT DO :D$
LET byt=byt+1 !!!
LET M=ORD(D$)
LOOP UNTIL M=255 ! 1st.mark
IF M<>255 THEN EXIT DO ! close & end_sub
CALL RED_D
LET M=ORD(D$)
LOOP
CLOSE #1
END SUB
!Page-4 へ続く
-
- [4]
Page-4
- 投稿者：SECOND
- 投稿日：2009年12月14日(月)20時43分29秒
-
-
!Page-4 の始め
!DRI
SUB FFDD
CALL RED_D
LET DRI=ORD(D$)*256
CALL RED_D
LET DRI=DRI+ORD(D$)
LET rct=DRI
END SUB
!APP0
SUB FFE0
FOR W=1 TO N
CALL RED_D
NEXT W
END SUB
!COMMENT
SUB FFFE
FOR W=1 TO N
CALL RED_D
NEXT W
END SUB
!DQT
SUB FFDB
DO WHILE 0< N
CALL RED_D
LET w= IP(ORD(D$)/16) !p=0(byte) p=1(word)
LET J=MOD(ORD(D$),16) !J=0~3 (QT.number)
FOR i=0 TO 63
CALL RED_D
LET DQ(U(i),V(i),J)=ORD(D$)
IF w=1 THEN
CALL RED_D
LET DQ(U(i),V(i),J)=DQ(U(i),V(i),J)*256+ORD(D$)
END IF
NEXT i
LET N=N-65-64*w ! remain size
LOOP
END SUB
!SOF0 SOF2
SUB FFC0
CALL RED_D
IF ORD(D$)<>8 THEN BREAK ! 8bit( 24bitColor ) at RGB.dimension
CALL RED_D
LET W=ORD(D$)*256
CALL RED_D
LET DY=W+ORD(D$) !V.pix.
CALL RED_D
LET W=ORD(D$)*256
CALL RED_D
LET DX=W+ORD(D$) !H.pix.
CALL RED_D
FOR i=0 TO ORD(D$)-1 !1~3 フレーム数、続くID 配置は、暗黙に Y,Cb,Cr の順
CALL RED_D
LET CoID(ORD(D$))=i ! CoID( ID=0~255) <--(Y=0, Cb=1, Cr=2)
CALL RED_D
LET MH(i)= IP(ORD(D$)/16) ! HV Y=11,12,21,22,41 Cb=11,11,11,11,11 Cr=11,11,11,11,11
LET MV(i)=MOD(ORD(D$),16)
CALL RED_D
LET QS(i)=ORD(D$) ! QT.number0~3 <-- QS( Y=0, Cb=1, Cr=2)
NEXT i
IF i=1 THEN LET CMO=0 ELSE LET CMO=2
END SUB
!DHT
SUB FFC4
DO WHILE 0< N
CALL RED_D
LET J=ORD(D$) ! (DC=0 AC=1 | ID0_0~3)
LET J=2*MOD(J,16)+IP(J/16) ! 0~1=ID0.DC~AC 2~3=ID1.DC~AC 4~5=ID2.…
LET DH(0,J)=0 !!!for 2nd.use for clear
FOR i=1 TO 16
CALL RED_D
LET DH(i,J)=ORD(D$)
LET DH(0,J)=DH(0,J)+DH(i,J)
NEXT i
FOR i=0 TO DH(0,J)-1
CALL RED_D
LET DV(i,J)=ORD(D$)
NEXT i
!---
FOR i=i TO 255
LET DV(i,J)=0
NEXT i
CALL makeH0(J) ! make Huffman Code table B() L()
CALL makeD0(J) ! make Huffman Decorder table A()
!---
LET N=N-1-16-DH(0,J) ! remain size
LOOP
END SUB
!SOS
SUB FFDA
CALL RED_D
LET M2=ORD(D$)
MAT HDC=(-2)*CON
MAT HAC=(-2)*CON
FOR i=1 TO M2
CALL RED_D ! ID=0~255( defined by SOFx)
LET w=ORD(D$)
CALL RED_D ! (DC_0~3|AC_0~3) huffman table selection
LET HDC(CoID(w))= IP(ORD(D$)/16)*2 !DC 0~3-->0,2,4,6
LET HAC(CoID(w))=MOD(ORD(D$),16)*2+1 !AC 0~3-->1,3,5,7
NEXT i
CALL RED_D
LET Ss_=ORD(D$) ! low of spectral selection
CALL RED_D
LET Se_=ORD(D$) ! high of spectral selection
CALL RED_D
LET Al=MOD(ORD(D$),16) !successive approximation bit position low ( point transform )
LET Ah=IP(ORD(D$)/16) !successive approximation bit position high ( preceding "Al" )
!--- balance monitor M30 M3() for display timing.
LET w=Se_-Ss_+1
IF Ah<>Al THEN LET w=w*(Ah-Al) ! prog.sa
FOR i=0 TO 2
IF 0<=HAC(i) THEN LET M3(i)=M3(i)+w ! M3()= scan band sum
NEXT i
IF CMO=0 OR M3(0)=M3(1) AND M3(1)=M3(2) THEN LET M30=M3(0) ELSE LET M30=99 ! Ybr.balance
!--- next image data top
END SUB
!プログレッシブ・シーケンスでは、段階的に画像ができるため、
!カラー画像表示のタイミング適正のため、上の様に、モニター変数 M30 を設けた。
!
! (Ah|Al) (バンド幅*ビット幅)の積算
!Ss Se Y Cb Cr ベース・ライン・シーケンス例 baseline
!00 3F 00 00 00 M3()= 64 64 64 …M30=64 完成画。
!Ss Se Y Cb Cr プログレッシブ・シーケンス例 spectral selection
!00 00 00 00 00 M3()= 1 1 1 …M30= 1
!01 05 -- 00 -- M3()= 1 6 1 = 99 99 は、バンド幅(Se-Ss+1)の
!01 05 -- -- 00 M3()= 1 6 6 = 99 累積が、
!01 05 00 -- -- M3()= 6 6 6 …M30= 6 Y,Cb,Cr で揃っていない時。
!06 3F -- 00 -- M3()= 6 64 6 = 99
!06 3F -- -- 00 M3()= 6 64 64 = 99
!06 3F 00 -- -- M3()= 64 64 64 …M30= 64 完成画。
!Ss Se Y Cb Cr プログレッシブ・シーケンス例 successive approximation
!00 00 01 01 01 M3()= -1 -1 -1 …M30= -1
!01 05 02 -- -- M3()= -11 -1 -1 = 99 99 は、バンド幅(Se-Ss+1)と、
!01 3F -- -- 01 M3()= -11 -1 -64 = 99 分割ビット長(Ah-Al)の積の
!01 3F -- 01 -- M3()= -11 -64 -64 = 99 累積が、
!06 3F 02 -- -- M3()=-127 -64 -64 = 99 Y,Cb,Cr で揃っていない時。
!01 3F 21 -- -- M3()= -64 -64 -64 …M30= -64
!00 00 10 10 10 M3()= -63 -63 -63 …M30= -63
!01 3F -- -- 10 M3()= -63 -63 0 = 99
!01 3F -- 10 -- M3()= -63 0 0 = 99
!01 3F 10 -- -- M3()= 0 0 0 …M30= 0 完成画。
SUB ROPEN
OPEN #1 :NAME FL$ ,ACCESS INPUT
END SUB
SUB RED_D
CHARACTER INPUT #1 :D$
LET byt=byt+1 !!!
END SUB
END
-
- [5]
改訂 十進 BASIC で、JPG ファイルを作る。
- 投稿者：SECOND
- 投稿日：2009年12月15日(火)04時53分6秒
-
-
!改訂 十進 BASIC で、JPG ファイルを作る。
!
!ハフマン符号木に、１つ空席を残さないと、開けないビューアの存在や、
!Linux 対策などが、先の投稿で、あちこち散乱していたのを、１つにし、
!Windows Linux 両用に、まとめた。
!
!このリスト最後の EXTERNAL SUB sample() で描く マンデルブロー の画像を、
!原画のサンプルにして、ベースラインJPG ファイルを作ります。
!DCT 変換、ハフマンコード化 …JPG ファイルまでを、BASIC だけで、実行。
!
!ハフマン・コードについては、標準テーブルを使わず、
!ランレングス・頻度の測定と、ハフマン・ツリーの作成を行い、それによる
!専用ハフマン・テーブルで、コード化します。（本来のハフマン・コード。）
!
DEBUG ON
!----------------------
! baseline JPG encoder
!----------------------
! テキスト・ウィンドウの、左上位置(x0,y0)と、幅(xw,yw)。
! CALL SetWindowPos( WinHandle("TEXT" ),0, 15,172,500,520, 0)
! SUB SetWindowPos( handle,C2, x0,y0,xw,yw, nFLG) !nFLG: 0=x0y0xwyw 1=x0y0 2=xwyw
! ASSIGN "user32.dll","SetWindowPos" !← Linux版 は使えない命令。
! END SUB
!----------------------------------------------------------
LET FL$="baseline.jpg" ! 削除すると、ダイアログ・ボックス入力。
SET ECHO "OFF"
ASK DIRECTORY s$
IF FL$>"" THEN PRINT "カレント DIR："& s$ ELSE file getname FL$, "jpg"
PRINT "出力ファイル："& FL$
IF FL$>"" THEN PRINT "上書き。又は作成されます。…"& "Ok？[Enter]"
IF FL$>"" THEN CHARACTER INPUT k$
IF FL$="" OR k$<>CHR$(13) THEN
PRINT "中止"
STOP
END IF
OPTION ARITHMETIC NATIVE
OPTION BASE 0
OPTION CHARACTER byte
SET TEXT background "OPAQUE"
ASK BITMAP SIZE bmx,bmy
SET WINDOW 0,bmx, bmy,0
!
DIM D8(1000,1000) !sample picture
DIM D2(1000,1000,2) !Y=D2(,,0) Cb=D2(,,1) Cr=D2(,,2)
DIM MH(2),MV(2) !MCU.Ybr.H()V()
DIM HDC(2),HAC(2) !hT.table selection
DIM QS(2),CoID(255) !qT.table selection
!
DIM U(63),V(63) !zigzag
DIM DQ(7,7,3) !DQT
DIM DH(16,7),DV(255,7) !DHT
DIM B(255+1,7),L(255,7) !huffman.code & length ( MAKE_H2 )
DIM B2(2) !Ybr D.C.成分 starting
DIM T(7,7),X(7),XO(7) !DDCT8X8
!
!---encorder
DIM SV(255,3) !ZFRE0 頻度SV( 座標)
DIM S_(255,3) !MAKE_DHT 頻度SV( 座標)--> 頻度S_( 降順No.) 座標DV( 降順No.)
DIM F_(255),Tr(16,255,3) !TREE3
!---lister
LET crlf$=CHR$(13)& CHR$(10)
DIM W_$(7)
LET W_$(0)="Y.DC"
LET W_$(1)="Y.AC"
LET W_$(2)="C.DC"
LET W_$(3)="C.AC"
!
SET VIEWPORT 0,0.4,0.6,1
CALL sample(DX$,DY$) !サンプルの描画。
LET DX=VAL(DX$)
LET DY=VAL(DY$)
MAT D8=ZER(DX-1,DY-1)
SET VIEWPORT 0, 1, 0, 1
SET WINDOW 0,bmx, bmy,0
! SET COLOR MODE "NATIVE" !← Linux版 は使えない命令。
ASK PIXEL ARRAY (0,0) D8
! MAT PLOT CELLS,IN 250,250; 450,450: D8 ! sample 原画 D8 のチェック
!
LET CMO=2 !CMO=0(mono.) CMO=2(color)
LET SD=1 !encoder 量子化テーブル調整 1/SD
CALL DQTINI !DQ(,,0)~DQ(,,1)~{zigzag U()V()},MH(),MV()
LET DU =CEIL(DX/(8*MH(0)))*8*MH(0) !Uwidth= (8X8)*2 bound by MCU size
LET DV_=CEIL(DY/(8*MV(0)))*8*MV(0) !Vwidth= (8X8)*1
MAT D2=ZER(DU-1,DV_-1,2) !Y=D2(,,0) Cb=D2(,,1) Cr=D2(,,2)
!
! CALL YbrRGB_n ! Ybr D2()<--RGB D8() !← color mode "NATIVE" の場合。
CALL YbrRGB_r ! Ybr D2()<--RGB D8() !← color mode "REGULAR" の場合。
LET MHT=1 ! flag. uncondition making huff.table
CALL DDCT8X8 ! D2() -->DCT -->Quantization
CALL ZZRL0 ! encoder
!---
PRINT "-------------------------"
PRINT "Encoder huffman Code"
FOR J=0 TO CMO+1
CALL list_HT(W_$(J),J) ! value Sort !H2.LST
NEXT J
beep
PRINT "終了"
!--------- JPG 色空間 -------------------------
! | Y | | 0.2990 +0.5870 +0.1140 | | R |
! |B-Y| = |-0.1687 -0.3313 +0.5000 | | G |
! |R-Y| | 0.5000 -0.4187 -0.0813 | | B |
!
! | R | | 1 0 +1.40200 | | Y |
! | G | = | 1 -0.34414 -0.71414 | |B-Y|
! | B | | 1 +1.77200 0 | |R-Y|
!----------------------------------------------
! Windows版 color mode "NATIVE" の場合。
!----------------------------------------
SUB YbrRGB_n !Ybr D2()<--RGB D8() 24bit Color
FOR V0=0 TO DY-1
FOR U0=0 TO DX-1
LET w1= MOD(D8(U0,V0),256) !R
LET w2=MOD(IP(D8(U0,V0)/256),256) !G
LET w3= IP(D8(U0,V0)/65536) !B
LET D2(U0,V0,0)= 0.2990*w1+0.5870*w2+0.1140*w3 !Y
LET D2(U0,V0,1)=-0.1687*w1-0.3313*w2+0.5000*w3 !Cb
LET D2(U0,V0,2)= 0.5000*w1-0.4187*w2-0.0813*w3 !Cr
NEXT U0
NEXT V0
END SUB
!---------------------------------------------------
! Windows版 or Linux版 color mode "REGULAR" の場合。
!---------------------------------------------------
SUB YbrRGB_r !Ybr D2()<--RGB D8() 色指標 Color
FOR V0=0 TO DY-1
FOR U0=0 TO DX-1
ASK COLOR MIX( D8(U0,V0)) w1,w2,w3 ! R,G,B (0~1)
LET D2(U0,V0,0)= 255*( 0.2990*w1+0.5870*w2+0.1140*w3) !Y
LET D2(U0,V0,1)= 255*(-0.1687*w1-0.3313*w2+0.5000*w3) !Cb
LET D2(U0,V0,2)= 255*( 0.5000*w1-0.4187*w2-0.0813*w3) !Cr
NEXT U0
NEXT V0
END SUB
!-----------
SUB DQTINI
RESTORE
!---DQT quantization
FOR j=0 TO 1
FOR V_=0 TO 7
FOR U_=0 TO 7
READ W
LET DQ(U_,V_,j)=CEIL(W/SD) ! inhibit 0
NEXT U_
NEXT V_
NEXT j
!---zigzag-U()V()
FOR V_=0 TO 7
FOR U_=0 TO 7
READ i
LET U(i)=U_
LET V(i)=V_
NEXT U_
NEXT V_
!---HT selection
MAT READ HDC !Y Cb Cr
MAT READ HAC !Y Cb Cr
!---QT selection
MAT READ QS !Y Cb Cr
!---MCU size
MAT MH=CON !Y Cb Cr =1
MAT MV=CON !Y Cb Cr =1
IF 0< CMO THEN
LET MH(0)=2 !Y
LET MV(0)=2 !Y
END IF
END SUB
!---quantization table
!輝度( SMPTE 370M ).Y
!*DQTY
DATA 32, 16, 17, 18, 18, 19, 42, 44
DATA 16, 17, 18, 18, 19, 38, 43, 45
DATA 17, 18, 19, 19, 40, 41, 45, 48
DATA 18, 18, 19, 40, 41, 42, 46, 49
DATA 18, 19, 40, 41, 42, 43, 48,101
DATA 19, 38, 41, 42, 43, 44, 98,104
DATA 42, 43, 45, 46, 48, 98,109,116
DATA 44, 45, 48, 49,101,104,116,123
!---quantization table
!色差( SMPTE 370M ).Cb Cr
!*DQTC
DATA 32, 16, 17, 25, 26, 26, 42, 44
DATA 16, 17, 25, 25, 26, 38, 43, 91
DATA 17, 25, 26, 27, 40, 41, 91, 96
DATA 25, 25, 27, 40, 41, 84, 93,197
DATA 26, 26, 40, 41, 84, 86,191,203
DATA 26, 38, 41, 84, 86,177,197,209
DATA 42, 43, 91, 93,191,197,219,232
DATA 44, 91, 96,197,203,209,232,246
!---Zigzag table
DATA 0, 1, 5, 6,14,15,27,28
DATA 2, 4, 7,13,16,26,29,42
DATA 3, 8,12,17,25,30,41,43
DATA 9,11,18,24,31,40,44,53
DATA 10,19,23,32,39,45,52,54
DATA 20,22,33,38,46,51,55,60
DATA 21,34,37,47,50,56,59,61
DATA 35,36,48,49,57,58,62,63
!---HT selection
DATA 0,2,2 !DC. Y Cb Cr
DATA 1,3,3 !AC. Y Cb Cr
!---QT selection
DATA 0,1,1 ! Y Cb Cr
!Page-2 へ続く
-
- [6]
Page-2
- 投稿者：SECOND
- 投稿日：2009年12月15日(火)04時55分18秒
-
-
!Page-2 の始め
!==========================================
! analizing frequency SV(nnnn,ssss) for DHT
SUB ZFRE0
MAT SV=ZER
LET B2(0)=0 ! Y.DC( start prediction)
LET B2(1)=0 ! Cb.DC
LET B2(2)=0 ! Cr.DC
!---
FOR V09=0 TO DV_-1 STEP 8*MV(0)
FOR U09=0 TO DU-1 STEP 8*MH(0)
!---MCU
FOR P=0 TO CMO ! ( 0=Y 1=Cb 2=Cr)
FOR V0=V09 TO V09+8*MV(P)-1 STEP 8
FOR U0=U09 TO U09+8*MH(P)-1 STEP 8
CALL F_BLK0
NEXT U0
NEXT V0
NEXT P
!---
NEXT U09
NEXT V09
!---
CALL MAKE_DHT !DH( ,J) DV( ,J) <--S_( ,J)
END SUB
! haffman Transform. main
SUB ZZRL0
! ---pass-1 analize frequency SV(nnnn,ssss) -->DV( ,J) DH( ,J)
IF 0< MHT OR DH(0,0)=0 THEN CALL ZFRE0
!---pass-2
CALL MAKE_H2 ! huffman code=B(V,J) len.=L(V,J) <-- DH( ,J) DV( ,J)
!---
LET byt=0 !!!
CALL WOPEN
CALL W_BIN31
!---
LET Hw=0 !bits stream buffer
LET BC=0 !bits in Hw
!---
LET B2(0)=0 ! Y.DC( start prediction)
LET B2(1)=0 ! Cb.DC
LET B2(2)=0 ! Cr.DC
!---
FOR V09=0 TO DV_-1 STEP 8*MV(0)
FOR U09=0 TO DU-1 STEP 8*MH(0)
!---MCU
FOR P=0 TO CMO !( 0=Y 1=Cb 2=Cr)
FOR V0=V09 TO V09+8*MV(P)-1 STEP 8
FOR U0=U09 TO U09+8*MH(P)-1 STEP 8
CALL W_BLK0
NEXT U0
NEXT V0
NEXT P
!---
NEXT U09
NEXT V09
CALL W_FLUSH
!---EOI
CALL WRT_H("FFD9")
CLOSE #1
PRINT "byte size:";byt
END SUB
!------
SUB F_BLK0
LET J=2*SGN(P) ! ( 0=Y 1=Cb 2=Cr)
!---D.C.part
LET W=D2(U0+U(0),V0+V(0),P)
LET SY=W-B2(P)
LET B2(P)=W ! previous of D.C.difference
IF SY<>0 THEN LET SS=LEN(BSTR$(ABS( SY),2)) ELSE LET SS=0 ! bit_length
LET SV(SS,J)=SV(SS,J)+1
!---A.C.parts
FOR AE=63 TO 0 STEP -1
IF 0<>D2(U0+U(AE),V0+V(AE),P) THEN EXIT FOR
NEXT AE
!---
LET Z=0 !zero run counter
FOR A_=1 TO AE
LET SY=D2(U0+U(A_),V0+V(A_),P)
IF SY=0 AND Z< 15 THEN
LET Z=Z+1
ELSE
IF SY<>0 THEN LET SS=LEN(BSTR$(ABS(SY),2)) ELSE LET SS=0 ! bit_length
LET W=Z*16+SS
LET Z=0
LET SV(W,J+1)=SV(W,J+1)+1
END IF
NEXT A_
IF A_< 64 THEN LET SV(0,J+1)=SV(0,J+1)+1 !End Of Block
END SUB
SUB W_BLK0
LET J=2*SGN(P) ! ( 0=Y 1=Cb 2=Cr)
!---D.C.part absolute SY bits length
LET W=D2(U0+U(0),V0+V(0),P)
LET SY=W-B2(P)
LET B2(P)=W ! previous of D.C.difference
IF SY<>0 THEN LET SS=LEN(BSTR$(ABS( SY),2)) ELSE LET SS=0 ! bit_length
LET L_=L(SS,J)
LET Ww=B(SS,J)
CALL W_HUFF
!---D.C.extent
IF SS<>0 THEN
IF SY< 0 THEN LET SY=SY+2^SS-1 !add maxim. in same bit_length.
LET L_=SS
LET Ww=SY*2^(16-L_)
CALL W_HUFF
END IF
!---A.C.parts
FOR AE=63 TO 0 STEP -1
IF 0<>D2(U0+U(AE),V0+V(AE),P) THEN EXIT FOR
NEXT AE
!---
LET Z=0 !zero run counter
FOR A_=1 TO AE
LET SY=D2(U0+U(A_),V0+V(A_),P)
IF SY=0 AND Z< 15 THEN
LET Z=Z+1
ELSE
IF SY<>0 THEN LET SS=LEN(BSTR$(ABS(SY),2)) ELSE LET SS=0 ! bit_length
LET W=Z*16+SS
LET Z=0
LET L_=L(W,J+1)
LET Ww=B(W,J+1)
CALL W_HUFF
!---A.C.extent
IF SS<>0 THEN
IF SY< 0 THEN LET SY=SY+2^SS-1 !add maxim. in same bit_length.
LET L_=SS
LET Ww=SY*2^(16-L_)
CALL W_HUFF
END IF
END IF
NEXT A_
IF A_< 64 THEN
LET L_=L(0,J+1)
LET Ww=B(0,J+1)
CALL W_HUFF !End Of Block
END IF
END SUB
!-----
!Ww：b15 ~0 左詰め 入力 bit_stream. L_：bit長
SUB W_HUFF
LET Hw=Hw+Ww*2^(-BC-8)
LET BC=BC+L_
DO WHILE 8<=BC
CALL WRT_D( IP(Hw))
IF IP(Hw)=255 THEN CALL WRT_D( 0)
LET Hw=FP(Hw)*256
LET BC=BC-8
LOOP
END SUB
!flush bit buffer with byte_bound( fill"1"in blank)
SUB W_FLUSH
IF BC<>0 THEN
LET w=Hw +2^(8-BC)-1
CALL WRT_D( w)
IF w=255 THEN CALL WRT_D( 0)
END IF
END SUB
!=====================
!pre hafman Transform.
!sort frequency S_() of nnnnssss( zero_run_length data)
!make DH() DV()
SUB MAKE_DHT
MAT DH=ZER
!---debug monitor
PRINT "-----------------------------------------"
FOR J=0 TO CMO+1 ! CMO =0=mono =2=color
PRINT "Zero_Run_Length 頻度表 (→)画素bit幅0~15、(↓)直前０数0~15"
PRINT W_$(J)
CALL msg00( SV, 0,255,J, 5) ! SV( 0~255, J)
PRINT "total=";Tx
NEXT J
PRINT
!---
FOR P=0 TO CMO+1 ! P(0~1=Y.DC~AC 2~3=C.DC~AC), CMO( 0=mono 2=color)
!--- make S_(,)DV(,)<-- SV(,)
LET SE=-1
FOR i=0 TO 255
IF SV(i,P)<>0 THEN
LET SE=SE+1
LET S_(SE,P)=SV(i,P)
LET DV(SE,P)=i
END IF
NEXT i
PRINT "===================================="
PRINT W_$(P)
PRINT "Zero Run Length 頻度表を詰めたもの"
CALL msg00( S_, 0,SE,P, 5) ! S_(0~SE, P)
PRINT "表座標(0~F=直前０数：0~F=画素bit長)"
CALL msg0x( DV, 0,SE,P, 5) ! DV(0~SE, P)
!---
CALL Qsort(0,SE)
CALL TREE3
!---
PRINT
PRINT " Encoder DHT table"
PRINT " (→)コード長1~16の、各個数"
CALL msg0x( DH, 1,16 ,P, 3) ! DH(1~16, P)
PRINT " 頻度順の、表座標(0~F=直前０数：0~F=画素bit長)"
CALL msg0x( DV, 0,Tx-1,P, 3) ! DV(0~Tx-1, P)
PRINT
NEXT P
END SUB
SUB msg00( M(,), S,E,J, w)
LET Tx=0
LET w$=""
FOR i=S TO E
LET Tx=Tx+M(i,J)
LET w$=w$& USING$( REPEAT$("#",w),M(i,J))
IF MOD(i-S,16)=15 THEN LET w$=w$& crlf$
NEXT i
IF MOD(i-S,16)=0 THEN PRINT w$; ELSE PRINT w$
END SUB
!Page-3 へ続く
-
- [7]
Page-3
- 投稿者：SECOND
- 投稿日：2009年12月15日(火)04時57分23秒
-
-
!Page-3 の始め
SUB msg0x( M(,), S,E,J, w)
LET Tx=0
LET w$=""
FOR i=S TO E
LET Tx=Tx+M(i,J)
LET w$=w$& REPEAT$(" ",w-2)& RIGHT$("0"& BSTR$(M(i,J),16),2)
IF MOD(i-S,16)=15 THEN LET w$=w$& crlf$
NEXT i
IF MOD(i-S,16)=0 THEN PRINT w$; ELSE PRINT w$
END SUB
!-------------------
! make huffman tree
!
!ハフマン符号木の、最後尾に１つ使わない枝を設けると、その分、符号木の
!枝ぶりが増え、全体のビット長も増え、圧縮コードとしては良くないですが、
!一部のビューアが、それを必要とする。（恐らくバグだと思います。）
!
!下の文中、" !←(+1) 符号木の最下に、空席を１つ作る。" …の行、２つ。
! SE+1 を、SE に戻すと空席は無くなり、本来の状態にも、戻せます。
!
SUB TREE3
MAT Tr=ZER
FOR i=0 TO SE
LET F_(i)=S_(i,P) !数値を壊すので、コピー F_(i)で実行
NEXT i
LET F_(SE+1)=0 !← 空席用
!---minimum pair
DO
LET w=1e8
FOR i=0 TO SE+1 !←(+1) 符号木の最下に、空席を１つ作る。
IF F_(i)< w THEN
LET w=F_(i)
LET Ad1=i ! minimum1 !頻度最小の分岐アドレスAd1
END IF
NEXT i
LET w=1e8
FOR i=0 TO SE+1 !←(+1) 符号木の最下に、空席を１つ作る。
IF F_(i)< w AND i<>Ad1 THEN
LET w=F_(i)
LET Ad2=i ! minimum2 !頻度最小の分岐アドレスAd2
END IF
NEXT i
IF w=1e8 THEN EXIT DO !分岐の組が無くなるまで
!---
LET F_(Ad1)=F_(Ad1)+F_(Ad2) !次の頻度最小の組探しは、２分岐合計を１つにし、
LET F_(Ad2)=1e9 !他方を外して行なう
!---
FOR Le1=16 TO 1 STEP -1 !アドレスAd1の最上 節点レベルLe1 を探す(最初のLe1=0)
IF Tr(Le1,Ad1,1)>0 OR Tr(Le1,Ad1,3)>0 THEN EXIT FOR
NEXT Le1
FOR Le2=16 TO 1 STEP -1 !アドレスAd2の最上 節点レベルLe2 を探す(最初のLe2=0)
IF Tr(Le2,Ad2,1)>0 OR Tr(Le2,Ad2,3)>0 THEN EXIT FOR
NEXT Le2
LET Le0=MAX( Le1,Le2 )+1 !両者何れよりも１つ上の節点レベル(Le0,Ad1)に、
!---
LET Tr(Le0,Ad1,0)=Le1 !分岐先( 節点レベル,アドレス)として２組記入
LET Tr(Le0,Ad1,1)=Ad1
LET Tr(Le0,Ad1,2)=Le2
LET Tr(Le0,Ad1,3)=Ad2
LOOP
!---make DH()
LET k=0
CALL bitl(Le0,Ad1) !全アドレスの nested 段数を求める。
FOR Ad=0 TO SE !nested 段数が同じ Ad の総数を、段数毎に集計
LET DH(Tr(0,Ad,1),P)=DH(Tr(0,Ad,1),P)+1
NEXT Ad
LET DH(0,P)=Ad
END SUB
SUB bitl(Le,Ad) !最上 節点(Le0,Ad1)より全分岐を、底まで辿る
IF 0< Le THEN
LET k=k+1
CALL bitl( Tr(Le,Ad,0), Tr(Le,Ad,1) ) !分岐 nested 1
CALL bitl( Tr(Le,Ad,2), Tr(Le,Ad,3) ) !分岐 nested 2
LET k=k-1
ELSE
LET Tr(0,Ad,1)=k !最上 節点から底までの nested 段数kを書く
END IF
END SUB
!-------------------
! Quick Sort S_()
SUB Qsort(L,R) ! 降順にセット。
local i,j
LET i=L
LET j=R
LET Tx=S_(IP((L+R)/2),P)
DO
DO WHILE S_(i,P) >Tx ! 降順＞、昇順＜
LET i=i+1
LOOP
DO WHILE Tx >S_(j,P) ! 降順＞、昇順＜
LET j=j-1
LOOP
IF j< i THEN EXIT DO ! 等号付 j<=i は、暴走。
SWAP S_(i,P),S_(j,P)
SWAP DV(i,P),DV(j,P)
LET i=i+1
LET j=j-1
LOOP UNTIL j< i ! 等号付 j<=i は、低速。
IF L< j THEN CALL Qsort(L,j)
IF i< R THEN CALL Qsort(i,R)
END SUB
!===================================
! make encorder table B()L()<-- DH()
SUB MAKE_H2
MAT L=ZER
FOR J=0 TO CMO+1
LET I=0 ! コード生成 順番(短い順)
LET Hx=0
LET Tx=BVAL("8000",16)
FOR L_=1 TO 16
FOR N=1 TO DH(L_,J)
LET V_=DV(I,J) ! 座標DV(頻度降順)
LET L(V_,J)=L_
LET B(V_,J)=Hx ! コード(座標V_)
LET I=I+1
LET Hx=Hx+Tx
NEXT N
LET Tx=Tx/2
NEXT L_
LET B(256,J)=1
NEXT J
END SUB
!==============================
!Fast Discrete Cosin Transform.( M=8x8, DCT-2 )
SUB DDCT8X8
FOR P=0 TO CMO ! (0=Y,1=Cb,2=Cr)
FOR V0=0 TO DV_-1 STEP 8*MV(0)/MV(P)
FOR U0=0 TO DU-1 STEP 8*MH(0)/MH(P)
FOR Y_=0 TO 7
LET w=Y_*MV(0) !sampling pt.Y
FOR X_=0 TO 7 !level shift, sampling CbCr from MCU
IF P=0 THEN LET X(X_)=D2(U0+X_,V0+Y_,P)-128 ELSE LET X(X_)=D2(U0+X_*MH(0),V0+w,P)
NEXT X_
CALL WANG
FOR U_=0 TO 7
LET T(U_,Y_)=X(U_)
NEXT U_
NEXT Y_
FOR U_=0 TO 7
FOR Y_=0 TO 7
LET X(Y_)=T(U_,Y_)
NEXT Y_
CALL WANG
FOR V_=0 TO 7
LET D2(U0+U_,V0+V_,P)=ROUND( X(V_)/DQ(U_,V_,QS(P)) ) ! Quantization
NEXT V_
NEXT U_
NEXT U0
NEXT V0
NEXT P
END SUB
!=============================
!Fast Discrete Cosin Transform
!Wang.( M=8, DCT-2 )
SUB WANG
LET XO(0)=X(0)+X(7)
LET XO(1)=X(1)+X(6)
LET XO(2)=X(2)+X(5)
LET XO(3)=X(3)+X(4)
LET XO(4)=X(3)-X(4)
LET XO(5)=X(2)-X(5)
LET XO(6)=X(1)-X(6)
LET XO(7)=X(0)-X(7)
!
LET X(0)=XO(0)+XO(3)
LET X(1)=XO(1)+XO(2)
LET X(2)=XO(1)-XO(2)
LET X(3)=XO(0)-XO(3)
LET X(4)=XO(7)*SQR(2)
LET X(5)=XO(6)-XO(5)
LET X(6)=XO(6)+XO(5)
LET X(7)=XO(4)*SQR(2)
!
LET XO(0)=(COS(PI/4)*X(0)+COS(PI/4)*X(1))
LET XO(1)=(COS(PI/4)*X(0)-COS(PI/4)*X(1)) !fin 0(0),1(4)
LET XO(2)=(COS(PI/8)*X(3)+SIN(PI/8)*X(2))
LET XO(3)=(COS(PI*3/8)*X(3)-SIN(PI*3/8)*X(2)) !fin 2(2),3(6)
LET XO(4)=X(4)
LET XO(5)=X(6)
LET XO(6)=X(5)
LET XO(7)=X(7)
!
LET X(4)=XO(4)+XO(5)
LET X(5)=XO(4)-XO(5)
LET X(6)=-XO(6)+XO(7)
LET X(7)=XO(6)+XO(7)
!
LET XO(4)=(COS(PI/16)*X(4)+SIN(PI/16)*X(7))
LET XO(5)=(COS(PI*5/16)*X(5)+SIN(PI*5/16)*X(6))
LET XO(6)=(SIN(PI*5/16)*X(5)-COS(PI*5/16)*X(6))
LET XO(7)=(SIN(PI/16)*X(4)-COS(PI/16)*X(7))
!
LET X(0)=SQR(2/8)*XO(0)
LET X(4)=SQR(2/8)*XO(1)
LET X(2)=SQR(2/8)*XO(2)
LET X(6)=SQR(2/8)*XO(3)
LET X(1)=SQR(1/8)*XO(4)
LET X(5)=SQR(1/8)*XO(5)
LET X(3)=SQR(1/8)*XO(6)
LET X(7)=SQR(1/8)*XO(7)
END SUB
!Page-4 へ続く
-
- [8]
Page-4
- 投稿者：SECOND
- 投稿日：2009年12月15日(火)04時59分7秒
-
-
!Page-4 の始め
!========================
SUB W_BIN31
!---SOI
CALL WRT_H("FFD8")
!---APP0
CALL WRT_H("FFE0")
CALL WRT_W( 2+14) !size
CALL WRT_M("JFIF"& CHR$(0))
CALL WRT_H("0102") !ver 1.2
CALL WRT_D(0) !0=NoUnit 1=dpi 2=dpcm
CALL WRT_W(100) !Xd 100 アスペクト比 IE6(no) imaging(ok)
CALL WRT_W(100) !Yd 100 アスペクト比
CALL WRT_D(0) !Xt 0
CALL WRT_D(0) !Yt 0
!---DQT.Y.C
CALL WRT_H("FFDB")
LET W=2 !start size
FOR J=0 TO CMO/2 !CMO=0(mono.) CMO=2(color)
LET W=W+1+64 !+ID +QT
NEXT J
CALL WRT_W( W) !size
FOR J=0 TO CMO/2 ! CMO=0(mono.) CMO=2(color)
CALL WRT_D(J) !J= (0)Y.DQ (1)C.DQ
FOR I=0 TO 63
CALL WRT_D( DQ(U(i),V(i),J) )
NEXT I
NEXT J
!---SOF0
CALL WRT_H("FFC0")
IF CMO=0 THEN CALL WRT_W( 11) ELSE CALL WRT_W( 17) !size
CALL WRT_D( 8) !8bit at RGB
CALL WRT_W( DY) !V.pixels
CALL WRT_W( DX) !H.pixels
IF CMO=0 THEN
CALL WRT_H("01") !1 items
CALL WRT_H("011100") !Y MCU(H:V) DQT
ELSE
CALL WRT_H("03") !3 items
CALL WRT_H("01"& STR$(MH(0))& STR$(MV(0))& "0"& STR$(QS(0))) !Y MCU(H:V) DQT
CALL WRT_H("02"& STR$(MH(1))& STR$(MV(1))& "0"& STR$(QS(1))) !Cb - -
CALL WRT_H("03"& STR$(MH(2))& STR$(MV(2))& "0"& STR$(QS(2))) !Cr - -
END IF
!---DHT
CALL WRT_H("FFC4")
LET W=2 !start size
FOR J=0 TO CMO+1 !CMO=0(mono.) CMO=2(color)
LET W=W+1+16+DH(0,J) !+ID +HT +VT
NEXT J
CALL WRT_W( W) !size
FOR J=0 TO CMO+1 !CMO=0(mono.) CMO=2(color)
CALL WRT_D( 16*MOD(J,2)+IP(J/2) ) ! 00h=YDC 10h=YAC 01h=CDC 11h=CAC
FOR I=1 TO 16 !(J) 0 =YDC 1 =YAC 2 =CDC 3 =CAC
CALL WRT_D( DH(I,J))
NEXT I
FOR I=0 TO DH(0,J)-1
CALL WRT_D( DV(I,J))
NEXT I
NEXT J
!---SOS
CALL WRT_H("FFDA")
IF CMO=0 THEN
CALL WRT_W(8) !size
CALL WRT_D(1) !1 items
CALL WRT_H("0100") !No. Y.HDC/HAC
ELSE
CALL WRT_W(12) !size
CALL WRT_D( 3) !3 items
CALL WRT_H("0100") !No. Y.HDC/HAC
CALL WRT_H("0211") !No. Cb.HDC/HAC
CALL WRT_H("0311") !No. Cr.HDC/HAC
END IF
CALL WRT_H("003F00") !band 0~63, AhAl=00
END SUB
!----------------
!open write binary
SUB WOPEN
OPEN #1:NAME FL$
ERASE #1
END SUB
!sequential write byte
SUB WRT_D( D)
PRINT #1 :CHR$(D);
LET byt=byt+1 !!!
END SUB
!sequential write word
SUB WRT_W( w)
PRINT #1 :CHR$(IP(w/256));CHR$(MOD(w,256));
LET byt=byt+2 !!!
END SUB
!sequential write binary massage W$
SUB WRT_M( w$)
FOR w=1 TO LEN(w$)
PRINT #1 :MID$(w$,w,1);
NEXT w
LET byt=byt+w-1 !!!
END SUB
!sequential write binary hex.massage w$
SUB WRT_H( w$)
FOR w=1 TO LEN(w$) STEP 2
PRINT #1 :CHR$(BVAL(MID$(w$,w,2),16));
NEXT w
LET byt=byt+LEN(w$)/2 !!!
END SUB
!=====================
! print huffman_table
SUB list_HT( m$,J)
PRINT m$;" 頻度 座標 bit コード(";
IF MOD(B(256,J),2)=1 THEN PRINT "座標順)" ELSE PRINT "生成順、頻度降順)"
LET sum=0
FOR i=0 TO 255
IF L(i,J)<>0 THEN
IF MOD(B(256,J),2)=1 THEN ! 1=Sort_value. Encorder
LET V_=i
LET w$=" "& RIGHT$(" "& STR$(SV(i,J)),4)& " " ! times
LET sum=sum+L(i,J)*SV(i,J)
ELSE ! 0=Sort_length. Decorder
LET V_=DV(i,J)
LET w$=" "& RIGHT$(" "& STR$(S_(i,J)),4)& " " ! times
LET sum=sum+L(i,J)*S_(i,J)
END IF
LET w$=w$& RIGHT$("0"& BSTR$(V_,16),2)& " " ! value
LET w$=w$& RIGHT$(" "& STR$(L(i,J)),2)& " " ! length
!--- huffman code
LET H_=B(i,J) ! code
LET L_=L(i,J) ! length
!--- bit_pattern
LET w$=w$& left$(right$("0000000"& BSTR$(H_,2),16),L_)
PRINT w$
END IF
NEXT i
PRINT " 合計( 頻度 * bit)=";sum
END SUB
END
EXTERNAL SUB sample(DX$,DY$)
! マンデルブロー（Complex\mandelbm.bas の着色改変）
OPTION ARITHMETIC COMPLEX
SET COLOR MODE "REGULAR"
SET POINT STYLE 1
FOR n=1 TO 51
SET COLOR MIX( n) 0 ,0 ,n/51 !BLACK < < BLUE
SET COLOR MIX( 51+n) 0 ,n/51 ,1 !BLUE < < CYAN
SET COLOR MIX(102+n) 0 ,1 ,1-n/51 !CYAN < < GREEN
SET COLOR MIX(153+n) n/51,1 ,0 !GREEN < < YELLOW
SET COLOR MIX(204+n) 1 ,1-n/51,n/51 !YELLOW< < MAGENTA
NEXT n
LET XL=-2
LET XR=.8
LET w1=XR-XL
LET w2=w1/2
SET WINDOW XL, XR,-w2,w2
ASK PIXEL SIZE(XL,-w2; XR,w2) px,py
!
FOR x=XL TO XR+1e-6 STEP w1/(px-1)
FOR y=-w2-.49*w1/(py-1) TO w2 STEP w1/(py-1) !故意に(x,0)を描点の間に挟む。
LET z=0
FOR n=1 TO 255
LET z=z^2+COMPLEX(x,y)
IF 2< ABS(z) THEN
IF n< 64 THEN SET POINT COLOR n*4 ELSE SET POINT COLOR 253
PLOT POINTS :x,y !上下の対象プロットをしない。
EXIT FOR
END IF
NEXT n
NEXT y
NEXT x
LET DX$=STR$(px)
LET DY$=STR$(py)
END SUB
-
- [9]
簡易 JPGファイル解析ツール
- 投稿者：SECOND
- 投稿日：2009年12月17日(木)03時40分41秒
-
-
!----------------------------
! 簡易 JPGファイル解析ツール
!
!表示される記号で、( ? | ? ) のように、真中を縦線で区切った( ) は、
!全体が１バイトで、左右 各々 4bit の数です。
!
!例、ハフマンテーブルの先頭１バイトの注釈は、(DC=0 AC=1| HT_0~3) 、
! 上位４bit が、DC/AC の区別、下位４bit は、標識番号の意味です。
!
!オフセットアドレスは、マーカー先頭にのみ、offset: 0x00… で、表示されます。
!左端の数値は、認識結果で区切られたデーター内容で、オフセットでありません。
!
!すき間の無い４桁数値は、16bit 数値ですが、バイト順は、逆になっておらず、
!ファイルデータ自体が、Big Endian です。
DEBUG ON
!----------------------------
!テキスト・ウィンドウの、左上位置(x0,y0)と、幅(xw,yw)。
CALL SetWindowPos( WinHandle("TEXT" ),0, 15,172,500,520, 0)
SUB SetWindowPos( handle,C2, x0,y0,xw,yw, nFLG) !nFLG: 0=x0y0xwyw 1=x0y0 2=xwyw
ASSIGN "user32.dll","SetWindowPos"
END SUB
!----------------------------------------------------------
OPTION ARITHMETIC NATIVE
OPTION CHARACTER byte
OPTION BASE 0
SET TEXT background "OPAQUE"
SET ECHO "OFF"
LET crlf$=CHR$(13)& CHR$(10)
LET q$=CHR$(34)
DIM d(500), m$( BVAL("C0",16) TO BVAL("FF",16) ), Ybr$(2), AhAl(20,0 TO 3)
DIM QS(2),CoID(255) !qT.<--color selection, color<--ID selection
MAT READ m$
MAT READ Ybr$
CALL zigzag
!
FILE GETNAME file$, "画像(*.jpg)|*.jpg"
IF file$="" THEN
PRINT "入力ファイル名が、ありません。"
STOP
END IF
PRINT "入力ファイル："& file$
CALL marker(file$)
CALL lstAhAl
IF M<>BVAL("D9",16) THEN
PRINT "Errors in File!,"
beep
END IF
PRINT "file end: 0x";BSTR$(byt-1,16)
SUB lstAhAl
IF 1< SOSn THEN
PRINT " (Ah|Al)"
PRINT "Ss Se Y Cb Cr プログレッシブ・シーケンスの表示"
FOR j=1 TO SOSn
LET w$=""
LET w$=w$& right$("0"& BSTR$(IP(AhAl(j,0)/256),16),2)& " "& right$("0"& BSTR$(MOD(AhAl(j,0),256),16),2)& " "
FOR i=1 TO 3
IF AhAl(j,i)< 0 THEN LET w$=w$& "-- " ELSE LET w$=w$& right$("0"& BSTR$(AhAl(j,i),16),2)& " "
NEXT i
PRINT w$
NEXT j
PRINT
END IF
END SUB
!------------
SUB marker(f$)
PRINT "----------------------------------------"
PRINT f$
OPEN #1 :NAME f$ ,ACCESS INPUT
LET byt=0 !offset counter
DO
DO
LET M=256 ! end of file
CHARACTER INPUT #1,IF MISSING THEN EXIT DO :D$
LET M=ORD(D$)
LET byt=byt+1 !offset counter
LOOP UNTIL M=255 ! 1st.mark
IF M<>255 THEN EXIT DO
CHARACTER INPUT #1 :D$
LET byt=byt+1 !offset counter
!---
LET M=ORD(D$)
IF M=BVAL("00",16) THEN ! 0xff in Picture data
LET ff0=ff0+1
ELSE
IF 0< SOSn AND 0< ff0 THEN PRINT right$(" "& STR$(ff0),3);"* 0xff00 in data"
!--- marker
PRINT
IF BVAL("C0",16)<=M THEN PRINT m$(M)(1:5);q$;m$(M)(6:10);q$; ELSE PRINT "FF";right$("0"& BSTR$(M,16),2);
PRINT TAB(14);"offset: 0x";BSTR$(byt-2,16)
!---
IF BVAL("D0",16)<=M AND M<=BVAL("D7",16) THEN ! RSI_0~7 in Picture data
ELSEIF M=BVAL("D8",16) THEN ! SOI
LET SOSn=0
MAT AhAl=(-1)*CON
ELSEIF M=BVAL("D9",16) THEN ! EOI
EXIT DO ! close#1, exit sub
ELSE
!---segment size
CHARACTER INPUT #1 :D$
LET w=ORD(D$)
CHARACTER INPUT #1 :D$
LET N=w*256+ORD(D$)-2 ! N=remain size
PRINT right$("000"& BSTR$(N+2,16),4);" size"
LET byt=byt+2 !offset counter
!---
IF M=BVAL("E0",16) THEN ! APP0
CALL FFE0_lst
ELSEIF M=BVAL("DB",16) THEN ! DQT
CALL FFDB_lst
ELSEIF m$(M)(6:8)="SOF" THEN ! SOF0~3 SOF5~7 SOF9~11 SOF13~15
CALL FFCx_lst
ELSEIF M=BVAL("C4",16) THEN ! DHT
CALL FFC4_lst
ELSEIF M=BVAL("DD",16) THEN ! DRI
CALL readD(N)
PRINT d(1)*256+d(2);"(decimal)"
ELSEIF M=BVAL("DA",16) THEN ! SOS
CALL FFDA_lst
!---
PRINT
PRINT "Image data, offset: 0x";BSTR$(byt,16)
LET ff0=0
ELSE
CALL skip_(N)
END IF
END IF
END IF
LOOP
CLOSE #1
PRINT
END SUB
!-----
SUB FFE0_lst ! APP0
CALL readW$(5)
PRINT q$;w$(1:4);q$ !"JFIF" ,00h
CALL readD(N-5)
IF w$(1:4)="JFIF" THEN
PRINT right$("0"& BSTR$(d(1),16),2);" ";right$("0"& BSTR$(d(2),16),2);
PRINT TAB(13);"Ver ";STR$(d(1));".";STR$(d(2))
PRINT right$("0"& BSTR$(d(3),16),2);" (0=No_unit 1=dot/inch 2=dot/cm)"
LET i=d(4)*256+d(5)
PRINT right$("000"& BSTR$(i,16),4);" ";i;"Xd アスペクト比" ! IE6(no_efect), imaging(active)
LET i=d(6)*256+d(7)
PRINT right$("000"& BSTR$(i,16),4);" ";i;"Yd アスペクト比"
PRINT right$("0"& BSTR$(d(8),16),2);" ";d(8);"Xt"
PRINT right$("0"& BSTR$(d(9),16),2);" ";d(9);"Yt"
FOR i=1 TO d(8)*d(9) ! (RGB)pixel*Xt*Yt
PRINT BSTR$(d(9+i)*65536+d(10+i)*256+d(11+i),16);" "
NEXT i
ELSEIF w$(1:4)="JFXX" THEN
PRINT BSTR$(d(1),16);" extension_code"
ELSE
PRINT "unknown file"
beep
STOP
END IF
END SUB
SUB FFDB_lst ! DQT
LET i=1
DO
CHARACTER INPUT #1 :D$
LET byt=byt+1 !offset counter
LET d(i)=ORD(D$)
LET w=IP(d(i)/16) !0=8bit 1=16bit
LET i=i+1
FOR s=0 TO 63
CHARACTER INPUT #1 :D$
LET byt=byt+1 !offset counter
LET j=ORD(D$)
IF w=1 THEN
CHARACTER INPUT #1 :D$
LET byt=byt+1 !offset counter
LET j=j*256+ORD(D$)
END IF
LET d( i+V(s)*8+U(s) )=j
NEXT s
LET i=i+s*(1+w)
LOOP UNTIL N< i
!---
LET i=1
DO
PRINT right$("0"& BSTR$(d(i),16),2);" (byte=0 word=1| QT_0~3)";
LET w$=""
LET j=i+1
FOR i=i+1 TO i+64
IF MOD(i-j,8)=0 THEN LET w$=w$& crlf$& " "
LET w$=w$& right$(REPEAT$(" ",3+w*2)& STR$(d(i)),4+w*2)
NEXT i
LET i=i+64*w
PRINT w$
LOOP UNTIL N< i
END SUB
SUB FFCx_lst ! SOF0~3 SOF5~7 SOF9~11 SOF13~15
CALL readD(N)
PRINT right$("0"& BSTR$(d(1),16),2);" precision_8~12bits in RGB source"
PRINT right$("0"& BSTR$(d(2),16),2);" ";right$("0"& BSTR$(d(3),16),2);" V.pix=";d(2)*256+d(3)
PRINT right$("0"& BSTR$(d(4),16),2);" ";right$("0"& BSTR$(d(5),16),2);" H.pix=";d(4)*256+d(5)
PRINT right$("0"& BSTR$(d(6),16),2);" number_in_flame_1~255"
FOR i=7 TO 7+3*(d(6)-1) STEP 3
LET co=(i-7)/3 !co_0~2= Ybr_0~2
LET CoID(d(i))=co !CoID( ID_0~255 )= co_0~2
PRINT right$("0"& BSTR$(d(i),16),2);" ";BSTR$(d(i+1),16);" ";BSTR$(d(i+2),16);" ID_0~255 (H|V)MCU QT …"& Ybr$(co)
LET QS(co)=d(i+2) !QT.numb= QS(co)
NEXT i
END SUB
SUB FFC4_lst ! DHT
CALL readD(N)
LET i=1
DO
PRINT right$("0"& BSTR$(d(i),16),2);" (DC=0 AC=1| HT_0~3)";
LET i=i+1
LET w=16
CALL dumph(i,w,16) ! d(i)~w bytes, return with(i=i+w, w=sum_d(i) )
CALL dumph(i,w,16) ! d(i)~w bytes, return with(i=i+w, w=sum_d(i) )
PRINT
LOOP UNTIL N< i
END SUB
!Page-2 へ続く
-
- [10]
Page-2
- 投稿者：SECOND
- 投稿日：2009年12月17日(木)03時43分35秒
-
-
!Page-2 の始め
SUB dumph(i,w,h) ! d(i)~w bytes, return with(i=i+w, w=sum_d(i) )
LET w$=""
LET j=i
FOR i=i TO i+w-1
IF MOD(i-j,h)=0 THEN LET w$=w$& crlf$& " " !right end, left end
LET w$=w$& " "& right$("0"& BSTR$(d(i),16),2) !Value.Table
IF i=j THEN LET w=d(i) ELSE LET w=w+d(i)
NEXT i
PRINT w$;
END SUB
SUB FFDA_lst !SOS
CALL readD(N)
LET SOSn=SOSn+1
PRINT right$("0"& BSTR$(d(1),16),2);" number_in_scan_1~4"
LET AhAl( SOSn,0)=d(N-2)*256+d(N-1)
FOR i=2 TO N-3 STEP 2
PRINT right$("0"& BSTR$(d(i),16),2);" ";right$("0"& BSTR$(d(i+1),16),2);" ID_0~255 (DC.HT_0~3|AC.HT_0~3) …"& Ybr$(CoID(d(i)))
LET AhAl( SOSn, CoID(d(i))+1 )=d(N) ! Ybr_0~2=CoID( ID_0~255 )
NEXT i
PRINT right$("0"& BSTR$(d(N-2),16),2);" ";right$("0"& BSTR$(d(N-1),16),2);" ";right$("0"& BSTR$(d(N),16),2);" Ss Se (Ah|Al)"
END SUB
SUB skip_(N)
FOR i=1 TO N
CHARACTER INPUT #1 :D$
NEXT i
LET byt=byt+N !offset counter
END SUB
SUB readD(N)
FOR i=1 TO N
CHARACTER INPUT #1 :D$
LET d(i)=ORD(D$)
NEXT i
LET byt=byt+N !offset counter
END SUB
SUB readW$(N)
LET w$=""
FOR i=1 TO N
CHARACTER INPUT #1 :D$
PRINT right$("0"& BSTR$(ORD(D$),16),2);" ";
LET w$=w$& D$
NEXT i
LET byt=byt+N !offset counter
END SUB
DATA "FFC0 SOF0" ,"FFC1 SOF1" ,"FFC2 SOF2" ,"FFC3 SOF3"
DATA "FFC4 DHT" ,"FFC5 SOF5" ,"FFC6 SOF6" ,"FFC7 SOF7"
DATA "FFC8 JPG" ,"FFC9 SOF9" ,"FFCA SOF10","FFCB SOF11"
DATA "FFCC DAC" ,"FFCD SOF13","FFCE SOF14","FFCF SOF15"
DATA "FFD0 RST0" ,"FFD1 RST1" ,"FFD2 RST2" ,"FFD3 RST3"
DATA "FFD4 RST4" ,"FFD5 RST5" ,"FFD6 RST6" ,"FFD7 RST7"
DATA "FFD8 SOI" ,"FFD9 EOI" ,"FFDA SOS" ,"FFDB DQT"
DATA "FFDC DNL" ,"FFDD DRI" ,"FFDE DHP" ,"FFDF EXP"
DATA "FFE0 APP0" ,"FFE1 APP1" ,"FFE2 APP2" ,"FFE3 APP3"
DATA "FFE4 APP4" ,"FFE5 APP5" ,"FFE6 APP6" ,"FFE7 APP7"
DATA "FFE8 APP8" ,"FFE1 APP9" ,"FFE2 APP10","FFE3 APP11"
DATA "FFEC APP12","FFE5 APP13","FFE6 APP14","FFE7 APP15"
DATA "FFF0 JPG0" ,"FFF1 JPG1" ,"FFF2 JPG2" ,"FFF3 JPG3"
DATA "FFF4 JPG4" ,"FFF5 JPG5" ,"FFF6 JPG6" ,"FFF7 JPG7"
DATA "FFF8 JPG8" ,"FFF9 JPG9" ,"FFFA JPG10","FFFB JPG11"
DATA "FFFC JPG12","FFFD JPG13","FFFE COM" ,"FFFF "
DATA "Y ","Cb","Cr"
SUB zigzag
DIM U(63),V(63) !zigzag
FOR V_=0 TO 7
FOR U_=0 TO 7
READ i
LET U(i)=U_
LET V(i)=V_
NEXT U_
NEXT V_
DATA 0, 1, 5, 6,14,15,27,28
DATA 2, 4, 7,13,16,26,29,42
DATA 3, 8,12,17,25,30,41,43
DATA 9,11,18,24,31,40,44,53
DATA 10,19,23,32,39,45,52,54
DATA 20,22,33,38,46,51,55,60
DATA 21,34,37,47,50,56,59,61
DATA 35,36,48,49,57,58,62,63
END SUB
END
-
- [11]
GIF ファイルの解析ツール Ver.22（画像付）
- 投稿者：SECOND
- 投稿日：2009年12月17日(木)09時51分44秒
-
-
! GIF ファイルの解析ツール Ver.22（画像付）
!-------
!先の投稿 http://6317.teacup.com/basic/bbs/bbs?M=RAL&&CID=464 では、
!大きな GIF ファイルの処理速度が、あまりに遅くなっていくので、
!一括処理を止め、細切れに処理するようにし、処理速度が同じ速さを保てる
!ようにしたものです。少し読み難い。
!
! Disposal Method 0,1,2,3 、インタレース画像、GIF アニメ動画などを、
! BASIC 文だけで、シンボリック・ダンプと、画像化を行います。
OPTION CHARACTER BYTE
LET bitfull=12 ! LZW max.bits (~~ 12)
DIM dic$(0 TO 2^bitfull) ! decoder 辞書
!
FILE GETNAME file$, "gif"
IF file$="" THEN
PRINT "入力ファイル名が、ありません。"
STOP
END IF
!
PRINT "入力ファイル："& file$
OPEN #1: NAME file$, ACCESS INPUT
PRINT "---------"
SET COLOR mode "native"
DIM ncp(0 TO 256) ! native color palette
LET ncp(256)=BVAL("eeeedd",16) ! 256= システム背景色 BGR
SET AREA COLOR ncp(256)
PLOT AREA: 0,0;1,0;1,1;0,1
CALL gif_head
LET i=MAX( Xsw,Ysw )
SET WINDOW -.1*i,1.1*i, 1.1*i,-.1*i
SET LINE STYLE 3
PLOT LINES:-1,-1;Xsw+.5,-1;Xsw+.5,Ysw+.5;-1,Ysw+.5;-1,-1
DIM m(0 TO Xsw-1,0 TO Ysw-1), m3(0 TO Xsw-1,0 TO Ysw-1)
MAT m=ncp(256)*CON
DO
CALL blocks_main
LOOP UNTIL b1$=CHR$(BVAL("3B",16))
PRINT "GIF 終端ブロック"
CALL dump(b1$,2,"block label")
PRINT "---------"
CLOSE #1
!----
SUB gif_head
LET h$=""
CALL readb( h$,13 )
IF h$(1:3)="GIF" THEN PRINT "GIF ヘッダー" ELSE CALL error
LET Xsw= ORD(h$( 8: 8))*256+ORD(h$(7:7))
LET Ysw= ORD(h$(10:10))*256+ORD(h$(9:9))
LET sflg=ORD(h$(11:11)) ! スクリーン情報のフラグ
LET BGco=ORD(h$(12:12))
LET aspect=ORD(h$(13:13))
LET b$=right$("0000000"& BSTR$(sflg,2) ,8)
LET colpix=2^(BVAL(b$(2:4),2)+1)
LET compal=2^(BVAL(b$(6:8),2)+1)
CALL dump(h$( 1: 6),6,"Asc") ! GIF識別文字
CALL dump(h$( 7: 8),2,"screen X_width= "& STR$(Xsw))
CALL dump(h$( 9:10),2,"screen Y_width= "& STR$(Ysw))
CALL dump(h$(11:11),2,"flags= "& b$ )
PRINT TAB(13);b$(1:1);":common_palette on=1/off=0"
PRINT TAB(11);b$(2:4);":colors/pixel 2^(";b$(2:4);"b+1)= ";STR$(colpix)
PRINT TAB(13);b$(5:5);":sort on=1/off=0 パレットの、重要度 色順ソート"
PRINT TAB(11);b$(6:8);":common_palette colors 2^(";b$(6:8);"b+1)= ";STR$(compal)
CALL dump(h$(12:12),2,"back_ground color= "& STR$(BGco))
IF aspect=0 THEN LET b$=".." ELSE LET b$=STR$(aspect)
CALL dump(h$(13:13),2,"アスペクト比 H:V=, 0 は 1:1 その他("& b$& "+15):64" )
CALL palette("common_パレット", c_p$, sflg)
END SUB
!----
SUB blocks_main
LET b1$=""
CALL readb( b1$, 1)
IF b1$=CHR$(BVAL("21",16)) THEN !追加データ・ブロック
CALL option_block
ELSEIF b1$=CHR$(BVAL("2C",16)) THEN !画像ブロック
CALL picture_block
ELSEIF b1$=CHR$(BVAL("3B",16)) THEN !GIF 終端ブロック
ELSE
CALL error
END IF
END SUB
!---追加データ・ブロック。
SUB option_block
PRINT "追加データ・ブロック"
CALL readb( b1$,1) ! w9$= readb_last_byte
IF w9$=CHR$(BVAL("F9",16)) THEN
LET im$=b1$
CALL dump(im$,2,"イメージコントロール・ブロック")
CALL blocks(im$,0,"",1) ! non print (data, byte/行, 注釈, 要求blocks)
LET iflg=ORD(im$(4:4)) ! フラグ
LET imtm=ORD(im$(6:6))*256+ORD(im$(5:5))
LET tcol=ORD(im$(7:7))
LET b$=right$("0000000"& BSTR$(iflg,2),8)
LET t_on=VAL(b$(8:8))
LET tmod=BVAL(b$(4:6),2)
CALL dump(im$(4:4),2,"flags= "& b$ )
PRINT TAB(11);b$(1:3);":blank"
PRINT TAB(11);b$(4:6);":"
PRINT TAB(15);"000= 描画後、そのまま、次の背景へ渡す"
PRINT TAB(15);"001= 描画後、そのまま、次の背景へ渡す"
PRINT TAB(15);"010= 描画後、back_ground color と取替、次の背景へ渡す"
PRINT TAB(15);"011= 描画後、描画前の画像を、次の背景へ渡す"
PRINT TAB(13);b$(7:7);":user click on=1/off=0"
PRINT TAB(13);b$(8:8);":透明色のスイッチ on=1/off=0"
CALL dump(im$(5:6),2,"フレーム表示時間= "& STR$(imtm)& " x10ms" )
CALL dump(im$(7:7),2,"透明色= "& STR$(tcol) )
CALL blocks(im$,16,"",999) ! (data, byte/行, 注釈, 要求blocks)
ELSEIF w9$=CHR$(BVAL("FE",16)) THEN
LET co$=b1$
CALL dump(co$,2,"コメント・ブロック")
CALL blocks(co$,8,"Asc",999) ! (data, byte/行, 注釈, 要求blocks)
ELSEIF w9$=CHR$(BVAL("FF",16)) THEN
LET ap$=b1$
CALL dump(ap$,2,"アプリケーション・ブロック")
CALL blocks(ap$,8,"Asc",1) ! (data, byte/行, 注釈, 要求blocks)
IF ap$(4:11)="NETSCAPE" THEN
CALL blocks(ap$,0,"",1) ! non print (data, byte/行, 注釈, 要求blocks)
CALL dump(ap$(16:16),2,"extension code= 01~07 next data type")
IF MOD(ORD(ap$(16:16)),8)=1 THEN
LET rept=ORD(ap$(18:18))*256+ORD(ap$(17:17))
CALL dump(ap$(17:18),2,"繰返し回数= "& STR$(rept)& " (0=endless)")
ELSE
CALL dump(ap$(17:LEN(ap$)),2,"02= 32bit buffering size, 03~07= ?")
END IF
END IF
CALL blocks(ap$,8,"Asc",999) ! (data, byte/行, 注釈, 要求blocks)
ELSEIF w9$=CHR$(BVAL("01",16)) THEN
LET tx$=b1$
CALL dump(tx$,2,"テキスト・イメージ・ブロック")
CALL blocks(d$,8,"Asc",999) ! (data, byte/行, 注釈, 要求blocks)
ELSE
CALL error
END IF
END SUB
!-------
SUB blocks(d$,m,t$,n) ! d$=data, m=byte/行, t$=注釈, n=要求blocks
FOR n=1 TO n
CALL readb(d$,1) ! w9$= readb_last_byte ! =block Size
CALL dump(w9$,1,"block size")
IF w9$=CHR$(0) THEN EXIT SUB ! block End
LET s=LEN(d$)
CALL readb(d$,ORD(w9$)) ! block data
IF 0< m THEN CALL dump(d$(s+1:LEN(d$)),m,t$)
NEXT n
END SUB
SUB readb(d$,cx) !cx=bytes size
FOR i=1 TO cx
CHARACTER INPUT #1,IF MISSING THEN EXIT FOR :w9$
LET d$=d$& w9$
NEXT i
IF i<=cx THEN CALL error
END SUB
SUB dump(d$,m,t$) ! t$="comment" →;comment t$="Asc" →;"ascii.dump"
FOR j=1 TO LEN(d$) STEP m
LET ww$=right$("0000"& BSTR$(adr,16),5)& " "
FOR i=j TO MIN(j+m-1, LEN(d$))
LET ww$=ww$& " "& right$("0"& BSTR$( ORD(d$(i:i)),16),2)
LET adr=adr+1
NEXT i
IF t$>"" THEN
LET ww$=ww$& REPEAT$(" ",6+3*m-LEN(ww$))
IF t$="Asc" THEN ! ascii.dump
LET ww$=ww$& " ;"""
FOR i=j TO MIN(j+m-1, LEN(d$))
IF " "<=d$(i:i) THEN LET ww$=ww$& d$(i:i) ELSE LET ww$=ww$& "."
NEXT i
LET ww$=ww$& """"
ELSE
IF m=3 THEN LET ww$=ww$& " ;"& STR$(IP(j/3)) ! パレット色番号
IF j=1 THEN LET ww$=ww$& " ;"& t$ ! comment
END IF
END IF
PRINT ww$ ! 行単位、テキスト画面のピカつき減少、高速。
NEXT j
END SUB
!-------
SUB picture_block
PRINT "画像ブロック"
LET pi$=b1$
CALL readb( pi$,9)
CALL dump(pi$(1:1),2,"block label")
LET Xp0=ORD(pi$(3:3))*256+ORD(pi$(2:2))
LET Yp0=ORD(pi$(5:5))*256+ORD(pi$(4:4))
LET Xpw=ORD(pi$(7:7))*256+ORD(pi$(6:6))
LET Ypw=ORD(pi$(9:9))*256+ORD(pi$(8:8))
CALL dump(pi$(2:3),2,"picture.X0_position left= "& STR$(Xp0))
CALL dump(pi$(4:5),2,"picture.Y0_position top= "& STR$(Yp0))
CALL dump(pi$(6:7),2,"picture.X_width= "& STR$(Xpw))
CALL dump(pi$(8:9),2,"picture.Y_width= "& STR$(Ypw))
LET pflg=ORD(pi$(10:10)) ! 画像情報のフラグ
LET b$=right$("0000000"& BSTR$(pflg,2),8)
CALL dump(pi$(10:10),2,"flags= "& b$ )
LET interl=VAL(b$(2:2))
LET pripal=2^(BVAL(b$(6:8),2)+1)
PRINT TAB(13);b$(1:1);":private_palette on=1/off=0"
PRINT TAB(13);b$(2:2);":interlace on=1/off=0, 1~step8 5~step8 3~step4 2~step2"
PRINT TAB(13);b$(3:3);":sort on=1/off=0 パレットの、重要度 色順ソート"
PRINT TAB(12);b$(4:5);":blank"
PRINT TAB(11);b$(6:8);":private_palette colors 2^(";b$(6:8);"b+1)= ";STR$(pripal)
CALL palette("private_パレット", p_p$, pflg)
!---
PRINT "画像データ"
LET LZW$=""
CALL readb( LZW$,1)
CALL dump(LZW$,1,"最小データ・ビット長")
CALL decomp_data
END SUB
!Page-2 へ続く
-
- [12]
Page-2
- 投稿者：SECOND
- 投稿日：2009年12月17日(木)09時53分33秒
-
-
!Page-2 の始め
SUB palette(n$, p$, pf)
PRINT n$;
LET p$=""
IF IP(pf/128)>0 THEN ! pf AND 0x80
PRINT
CALL readb(p$, 3*2^(MOD(pf,8)+1)) ! pf AND 0x07
CALL dump(p$,3,"R G B")
CALL palette10(n$, p$, pf)
ELSE
LET ww$=" 無し。"
IF n$(1:1)="p" AND c_p$>"" THEN
LET ww$=ww$& "→ common_パレット 使用。"
IF bk_p$(1:1)<>"c" THEN CALL palette10("common", c_p$, sflg)
END IF
PRINT ww$
END IF
END SUB
SUB palette10(n$, p$, pf)
FOR i=1 TO 3*2^(MOD(pf,8)+1) STEP 3
LET ncp(IP(i/3))= ORD(p$(i:i))+256*ORD(p$(i+1:i+1))+65536*ORD(p$(i+2:i+2))
NEXT i
LET bk_p$=n$
END SUB
SUB error
beep
PRINT "File Error Stop"
STOP
END SUB
!------------------- 画像データ LZW$ から、原画を見る。------------------
SUB decomp_data
LET obits0=ORD(LZW$(1:1))
LET N000=2^obits0
LET li=2 !start input pointer
LET blkend=0 !clear input block pointer
LET iacc$="" !clear input bit buffer
LET pdata$="" !clear output byte buffer
!---
IF tmod=3 THEN MAT m3=m
LET p_=1
IF interl=0 THEN
CALL intlace( 0, 1) ! start_raster, step !---NO interlace
ELSE
CALL intlace( 0, 8) ! start_raster, step !---0~7step8
WAIT DELAY .5
CALL intlace( 4, 8) ! start_raster, step !---4~7step8
WAIT DELAY .5
CALL intlace( 2, 4) ! start_raster, step !---2~3step4
WAIT DELAY .5
CALL intlace( 1, 2) ! start_raster, step !---1~1step2
END IF
IF tmod=2 THEN MAT m=ncp(256)*CON ! 256= システム背景色
IF tmod=3 THEN MAT m=m3
PRINT "******************* 復元終り"
END SUB
SUB intlace( ss, stp)
FOR j_=Yp0 TO Yp0+Ypw-1
IF MOD(j_-Yp0, stp) >=ss THEN
IF MOD(j_-Yp0, stp) >ss THEN LET p_=p_-Xpw
!---
IF LEN(pdata$)-Xpw< p_ THEN
LET pdata$=pdata$(p_:LEN(pdata$))
CALL LZW_decoder
LET p_=1
END IF
!---
FOR i_=Xp0 TO Xp0+Xpw-1
WHEN EXCEPTION IN
LET col=ORD(pdata$(p_:p_))
IF t_on<>1 OR tcol<>col THEN LET m(i_,j_)=ncp(col)
LET p_=p_+1
USE
END WHEN
NEXT i_
END IF
NEXT j_
MAT PLOT CELLS,IN 0,0; Xsw-1,Ysw-1 :m
END SUB
!-------------
SUB LZW_decoder
DO
LET dicnum=N000+2-1 !start dic.number-1
DO
LET iwidth=LEN( BSTR$(dicnum,2) ) !remake iwidth
IF bitfull< iwidth THEN LET iwidth=bitfull !to handle BAD file
LET dicnum=dicnum+1
CALL inpcode !data on bx
IF bx=-1 OR bx=N000+1 THEN
EXIT SUB
ELSEIF bx=N000 THEN
EXIT DO
ELSE
IF bx< N000 THEN
LET dic$(dicnum)=CHR$(bx)
ELSE
LET dic$(dicnum)=dic$(bx)
LET dic$(dicnum)=dic$(dicnum)& dic$(bx+1)(1:1)
END IF
LET pdata$=pdata$& dic$(dicnum)
END IF
LOOP
IF Xpw< LEN(pdata$) THEN EXIT SUB
LOOP
END SUB
SUB inpcode
LET bx=-1
DO WHILE LEN(iacc$)< iwidth
IF blkend<=li THEN
!--- LZW データ（size:data, size:data, … 0 ）
IF LEN(LZW$)< li THEN
IF 2< li THEN PRINT "******************* 復元休止"
LET LZW$=""
CALL blocks(LZW$,16,"",8) ! (data, byte/行, 注釈, 要求blocks)
LET li=1
PRINT "******************* 復元( LZW デコード )"
END IF
!---
LET blksize=ORD(LZW$(li:li))
CALL prthex(CHR$(blksize),"block size") !モニター
IF blksize=0 THEN EXIT SUB
LET li=li+1
LET blkend=li+blksize
END IF
LET iacc$=right$("0000000"& BSTR$(ORD(LZW$(li:li)),2),8)& iacc$
LET li=li+1
LOOP
LET bx=BVAL(right$(iacc$,iwidth),2)
LET iacc$=iacc$(1:LEN(iacc$)-iwidth)
END SUB
SUB prthex(d$,w$)
LET ww$=""
FOR i9=1 TO LEN(d$)
LET ww$=ww$& right$("0"& BSTR$(ORD(d$(i9:i9)),16),2)& " "
NEXT i9
IF w$>"" THEN LET ww$=ww$& "; "& w$
PRINT ww$
END SUB
END
!関連投稿（１ページ未満）
![メインスレッド458] GIF アニメ－ション を作る。
![メインスレッド460] LZW エンコーダーと、デコーダー
-
- [13]
ＦＦＴ プログラム
- 投稿者：SECOND
- 投稿日：2010年 3月25日(木)05時05分48秒
-
-
DEBUG ON
!******************************************************************
! ＦＦＴ プログラム
!******************************************************************
!------------------------------------------------------------------------------
! デジタルでのフーリエ変換。 逆変換。
!
! ◆DFT Discrete Fourier Transform ◆IDFT Inverse Discrete Fourier Transform
!
! N-1 N-1
! F(k)= ∑X(r)*exp(-j2π* k*r/N ) X(r)= ( 1/N )* ∑F(k)*exp( j2π* r*k/N )
! r=0 k=0
!
! ※k= 0～N-1 ※r= 0～N-1
!
!------------------------------------------------------------------------------
! DFT が Ｎ個の周期波を扱う場合は、その高速算法がある。
! Cooley-Tukey(クーリー・テューキー)型 Fast Fourier Transform
!
! 入力信号の番号を、奇数 偶数の 時間間引き ２分割にする方法と、
! 出力信号の番号を、奇数 偶数の 周波数間引き２分割にする方法の、２通り。
!
! 周期 N の DFT 計算量は、N*N 個の積和であるが、N/2 個の DFT ２つに
! 分ける事が出来れば、(N/2)*(N/2) 個の積和２つになって、ざっとみても
! １回で、半分に計算量が減る。これを、限界まで進めたもの。
! ( N=512 程度の実用サイズで、計算量が 約1/100 に。詳細は、以下)
!----------------------------------------
! FFT 部分の初期設定
!----------------------------------------
OPTION ARITHMETIC COMPLEX
OPTION BASE 0
!
LET N1=8
LET N1=16
!LET N1=32
!LET N1=64
!LET N1=128
!LET N1=256
!LET N1=512
!LET N1=1024
!LET N1=2048
!LET N1=4096
!LET N1=65536
!
DIM Z(N1*2) !FFT 信号入力、変換出力を、Z() 配列だけの双方向で行なう。
! !データは N1 個であるが、ビットリバースを使わない型だけが、
! !同じサイズの後部余白を、必要とする。Clear は不要。
!
DIM ZO(N1*2) !DFT だけが使用。
!
DIM W19(N1)
LET W2=2*PI/N1
FOR I=0 TO N1
LET w19(I)=EXP(COMPLEX(0,W2*I))
NEXT I
!------------------------------------------------------------------------------
! ◆ 時間間引き型 Decimation In Time
!
! 入力番号を、奇数 偶数の 時間間引き２分割、出力番号を、連続の周波数２分割
! とした ２つの DFT が、合計で、元の DFT と同じになる様に、書くと、
!
! ※W( a/b )=exp(-j2π *a/b )
! N/2-1 N/2-1
! F(k) = ∑X(2r)*W( k*2r/N ) + ∑X(2r+1)*W( k*(2r+1)/N ) ・・・①
! r=0 r=0
! ※k= 0～N/2-1
! N/2-1 N/2-1
! F(k+N/2)= ∑X(2r)*W( (k+N/2)*2r/N ) + ∑X(2r+1)*W( (k+N/2)*(2r+1)/N ) ・・・②
! r=0 r=0
!
!------------------------------------------------------------------------------
! 上の式を、それぞれ整理したもの。 ※W( 整数 )=1, W( 整数＋0.5 )=-1
!
! N/2-1 N/2-1
! F(k) = ∑X(2r)*W( k*r/(N/2) ) + W(k/N)*∑X(2r+1)*W( k*r/(N/2) ) ・・・①'
! r=0 r=0
! ※k= 0～N/2-1
! N/2-1 N/2-1
! F(k+N/2)= ∑X(2r)*W( k*r/(N/2) ) - W(k/N)*∑X(2r+1)*W( k*r/(N/2) ) ・・・②'
! r=0 r=0
!
!------------------------------------------------------------------------------
! 半分ずつの X(2r) と X(2r+1) それぞれを新しい入力
! とする 基数=N/2 のDFT 出力の合成になっている。
! この DFT も、さらに２分割し、入れ子構造にしていく。
!
! ２分割を、N(基数)=2 、入力側が、基数=1 の DFT 出力になるまで繰返すと、
! 下の様になる。 ※添え字は、N(基数)=2 の DFT を基準に付け直している。
!
! 基数=1 の出力とは、入力それ自身なので、
! 入力 X(0) X(1) を直接、基数=1 の DFT 出力とみなせる。結果、F(0)=X(0)+X(1)
! F(1)=X(0)-X(1)
! となって、これは、N(基数)=2 のまま、DFT 変換する状況とも同一になる。
! 結局、FFT は、殆どDFT 計算を、しないまま分割だけで、その結果を得ている。
!
! 0 0
! F(0)= ∑X(0)*W(0/1) + W(0/2)*∑X(1)*W(0/1) ・・・①'
! r=0 r=0
!
! 0 0
! F(1)= ∑X(2r)*W(0/1) - W(0/2)*∑X(1)*W(0/1) ・・・②'
! r=0 r=0
!
!------------------------------------------------------------------------------
! 実際のプログラム
!
! ◆ 時間間引き型 Decimation In Time
! N(基数)=8 の場合の例。 数式は右端から、(← 方向)へ２分割を繰返すが、
! 実計算は左端から、(→ 方向)へ
! ※ Wa/b はexp(-j2π*a/b)を乗ずる意。W1/4 は右へ1/4 回転。W1/8 は右へ1/8回転。
! ① は負号を付ける意。 ╋ は加算の意。
!
! 入力信号の順を、番号の上位下位ビットを逆に読む順 (Index Bit Reverse) に入替る。
! ↓
! X0 X0 F0 F0 F0
! ────┬─╋─────┬────╋─────┬──────────╋─
! → ＼／ ＼ ／ ＼ ／
! X1 X4 W0/2 ／＼ F4 ＼ ／ F2 ＼ ／ F1
! ────┴①╋─────┬────╋─────┬──────────╋─
! ＼／ ＼／ ＼ ＼ ／ ／
! X2 X2 F2 W0/4 ／＼ ／＼ F4 ＼ ＼ ／ ／ F2
! ────┬─╋─────┴───①╋─────┬──────────╋─
! ＼／ ／ ＼ ＼ ＼／ ＼／ ／
! X3 X6 W0/2 ／＼ F6 W1/4 ／ ＼ F6 ＼ ／＼ ／＼ ／ F3
! ────┴①╋─────┴───①╋─────┬──────────╋─
! ＼／ ＼／ ＼／ ＼／
! X4 X1 F1 F1 W0/8 ／＼ ／＼ ／＼ ／＼ F4
! ────┬─╋─────┬────╋─────┴─────────①╋─
! ＼／ ＼ ／ ／ ＼／ ＼／ ＼
! X5 X5 W0/2 ／＼ F5 ＼ ／ F3 W1/8 ／ ／＼ ／＼ ＼ F5
! ────┴①╋─────┬────╋─────┴─────────①╋─
! ＼／ ＼／ ／ ／ ＼ ＼
! X6 X3 F3 W0/4 ／＼ ／＼ F5 W2/8 ／ ／ ＼ ＼ F6
! ────┬─╋─────┴───①╋─────┴─────────①╋─
! ＼／ ／ ＼ ／ ＼
! X7 X7 W0/2 ／＼ F7 W1/4 ／ ＼ F7 W3/8 ／ ＼ F7
! ────┴①╋─────┴───①╋─────┴─────────①╋─
! ┠基数２のＤＦＴ─╂←─ ネスト(Nested)された周波数出力の２分割の重ね ─→┨
!
!
! 計算プログラムは、和と差の１対 の繰返しであり、同じ処理で実行するには、
! 各段の計算結果が、次段の計算の入力にもなるため、同じ場所に書き戻す必要があった。
!
! X0~,F0~ の添え字を、直接の物理アドレスとして、配列を組むと、書き込む前に有った
! データーが使用される前に、上書きされてしまう不都合が発生する。その解決に、
!
! X0~,F0~ の添え字を、論理アドレスとして無視、図の上から下への配置を物理アドレス
! として、配列を用いる計算方法をとる。その為には、
!
! 左の入力信号内容 X0~ を、X0,X4,X2,, の順に入換えておかなければならない。
! この操作を、Index Bit Reverse と呼ぶ。
!
! X0~,F0~ の添え字を、直接の物理アドレスとしても、上書されないように出来れば、
! Index Bit Reverse は、不要なので、その例を最後に紹介する。
!
!---------INDEX BIT REVERSE----
SUB IBrev ! Z() を、その添え字 番号の上位下位ビット順が、
LET J=0 ! 逆になった番号 の位置に並べ替える。
LET N2=N1/2
FOR I=0 TO N1-2
IF I< J THEN swap Z(J),Z(I)
!----- J=000.. 100.. 010.. 110..ビット逆順のＵＰカウンター
LET K=N2 !K の初期値、100... 最上位ビットだけの数
DO WHILE K<=J !J の K ビット位置が 1 である場合。
LET J=J-K !J の桁上りの為、そのビットを消す
LET K=K/2 !K を、次のビットへ。
LOOP !J の K ビット位置が 0 で loop を出る。
LET J=J+K !J の K ビット位置に 1 を置く
!-----
NEXT I
END SUB
!=========================================
! 高速フーリエ変換 Fast Fourier Transform
! 時間間引き型 DIT (Decimation In Time)
!
SUB FFT_DIT
CALL IBrev !Index Bit Reverse
!-------------
!CALL ref_DIT(0,N1) !再帰型
CALL normal_DIT !通常型（高速)
END SUB
!再帰型（通常型より約１.５倍に遅くなるが、簡潔。)
SUB ref_DIT(J,N)
IF 1< N THEN
CALL ref_DIT(J ,N/2)
CALL ref_DIT(J+N/2,N/2)
LET H=N1 !FFT: N1 IFFT: 0
FOR I=J TO J+N/2-1
LET IK=I+N/2
LET ZT=w19(H)*Z(IK)
LET Z(IK)=Z(I)-ZT !F(k)-Wk/N*F(k+N/2)
LET Z(I)= Z(I)+ZT !F(k)+Wk/N*F(k+N/2)
LET H=H-N1/N !FFT: -N1/N IFFT: +N1/N
NEXT I
END IF
END SUB
!----
!通常型（高速)
SUB normal_DIT
LET K=1
DO WHILE K< N1
LET H=N1 !FFT: N1 IFFT: 0
LET K2=K*2
LET D=N1/K2
FOR J=0 TO K-1
FOR I=J TO N1-1 STEP K2
LET IK=I+K
LET ZT=w19(H)*Z(IK)
LET Z(IK)=Z(I)-ZT !F(k)-Wk/N*F(k+N/2)
LET Z(I)= Z(I)+ZT !F(k)+Wk/N*F(k+N/2)
NEXT I
LET H=H-D !FFT: -D IFFT: +D
NEXT J
LET K=K2
LOOP
END SUB
!Page-2 へ続く
-
- [14]
Page-2
- 投稿者：SECOND
- 投稿日：2010年 3月25日(木)05時08分52秒
-
-
!Page-2 の始め
!=========================================
! 高速フーリエ 逆変換 Inverse Fast Fourier Transform
! 時間間引き型 DIT (Decimation In Time) と同文で、
! 回転子W だけを 逆回しにしたもの。内容は周波数間引きになる。
!
SUB IFFT_DIT
CALL IBrev !Index Bit Reverse
!-------------
!CALL ref_IDIT(0,N1) !再帰型
CALL normal_IDIT !通常型（高速)
!-------------
MAT Z=(1/N1)*Z
END SUB
!再帰型（通常型より約１.５倍に遅くなるが、簡潔。)
SUB ref_IDIT(J,N)
IF 1< N THEN
CALL ref_IDIT(J ,N/2)
CALL ref_IDIT(J+N/2,N/2)
LET H=0 !FFT: N1 IFFT: 0
FOR I=J TO J+N/2-1
LET IK=I+N/2
LET ZT=w19(H)*Z(IK)
LET Z(IK)=Z(I)-ZT !F(k)-Wk/N*F(k+N/2)
LET Z(I)= Z(I)+ZT !F(k)+Wk/N*F(k+N/2)
LET H=H+N1/N !FFT: -N1/N IFFT: +N1/N
NEXT I
END IF
END SUB
!----
!通常型（高速)
SUB normal_IDIT
LET K=1
DO WHILE K< N1
LET H=0 !FFT: N1 IFFT: 0
LET K2=K*2
LET D=N1/K2
FOR J=0 TO K-1
FOR I=J TO N1-1 STEP K2
LET IK=I+K
LET ZT=w19(H)*Z(IK)
LET Z(IK)=Z(I)-ZT !F(k)-Wk/N*F(k+N/2)
LET Z(I)= Z(I)+ZT !F(k)+Wk/N*F(k+N/2)
NEXT I
LET H=H+D !FFT: -D IFFT: +D
NEXT J
LET K=K2
LOOP
END SUB
!------------------------------------------------------------------------------
! ◆ 周波数間引き型 Decimation In Frequency
!
! 入力番号を、連続の時間２分割、出力番号を、奇数 偶数の 周波数間引き２分割
! とした ２つのＤＦＴが、合計で、元のＤＦＴと同じになる様に、書くと、
!
! ※W( a/b )=exp(-j2π *a/b )
! N/2-1 N/2-1
! F(2k) = ∑X(r)*W( 2k*r/N ) + ∑X(r+N/2)*W( 2k*(r+N/2)/N ) ・・・①
! r=0 r=0
! ※k= 0～N/2-1
! N/2-1 N/2-1
! F(2k+1)= ∑X(r)*W( (2k+1)*r/N ) + ∑X(r+N/2)*W( (2k+1)*(r+N/2)/N ) ・・・②
! r=0 r=0
!
!------------------------------------------------------------------------------
! 上の式を、それぞれ整理したもの。 ※W( 整数 )=1, W( 整数＋0.5 )=-1
!
! N/2-1
! F(2k) = ∑( X(r)+X(r+N/2) ) *W( k*r/(N/2) ) ・・・①'
! r=0
! ※k= 0～N/2-1
! N/2-1
! F(2k+1)= ∑( X(r)-X(r+N/2) )*W(r/N) *W( k*r/(N/2) ) ・・・②'
! r=0
!
!------------------------------------------------------------------------------
! 半分ずつの X(r)+X(r+N/2) と ( X(r)-X(r+N/2) )*W(r/N) それぞれを新しい入力
! とする基数=N/2 のDFT 出力 ２つになっている。
! この DFT も、さらに２分割し、入れ子構造にしていく。
!
! ２分割を、N(基数)=2 、和入力、差入力が、各々１つになるまで繰返すと、
! 下の様になる。 ※添え字は、N(基数)=2 の DFT を基準に付け直している。
!
! これは、基数=1 の DFT で、その入力自身が、その出力になり、
! 入力 X(0)+X(1), (X(0)-X(1))*W(0/2) を直接 DFT 出力とみなせる。
! 結果、F(0)=X(0)+X(1)
! F(1)=X(0)-X(1)
! となって、これは、N(基数)=2 のまま、DFT 変換する状況とも同一になる。
! 結局、FFT は、殆どDFT 計算を、しないまま、分割だけで、終っています。
!
! 0
! F(0)= ∑( X(0)+X(1) ) *W(0/1) ・・・①'
! r=0
!
! 0
! F(1)= ∑( X(0)-X(1) )*W(0/2) *W(0/1) ・・・②'
! r=0
!
!------------------------------------------------------------------------------
! 実際のプログラム
!
! ◆ 周波数間引き型 Decimation In Frequency
! N(基数)=8 の場合の例。 数式は左端から、(→ 方向)へ２分割を繰返す。
! 実計算も左端から、(→ 方向)へ
! ※ Wa/b はexp(-j2π*a/b)を乗ずる意。W1/4 は右へ1/4 回転。W1/8 は右へ1/8回転。
! ① は負号を付ける意。 ╋ は加算の意。
!
! 番号の上位下位ビットを逆に読む順 (Index Bit Reverse) に、出力信号の順を入替る。
! ↓
! X0 X0 X0 F0 F0
! ─┬──────────╋─────┬────╋─────┬─╋────
! ＼ ／ ＼ ／ ＼／ →
! X1 ＼ ／ X2 ＼ ／ X4 ／＼ W0/2 F4 F1
! ─┬──────────╋─────┬────╋─────┴①╋────
! ＼ ＼ ／ ／ ＼／ ＼／
! X2 ＼ ＼ ／ ／ X4 ／＼ ／＼ W0/4 X2 F2 F2
! ─┬──────────╋─────┴───①╋─────┬─╋────
! ＼ ＼／ ＼／ ／ ／ ＼ ＼／
! X3 ＼ ／＼ ／＼ ／ X6 ／ ＼ W1/4 X6 ／＼ W0/2 F6 F3
! ─┬──────────╋─────┴───①╋─────┴①╋────
! ＼／ ＼／ ＼／ ＼／
! X4 ／＼ ／＼ ／＼ ／＼ W0/8 X1 X1 F1 F4
! ─┴─────────①╋─────┬────╋─────┬─╋────
! ／ ＼／ ＼／ ＼ ＼ ／ ＼／
! X5 ／ ／＼ ／＼ ＼ W1/8 X3 ＼ ／ X5 ／＼ W0/2 F5 F5
! ─┴─────────①╋─────┬────╋─────┴①╋────
! ／ ／ ＼ ＼ ＼／ ＼／
! X6 ／ ／ ＼ ＼ W2/8 X5 ／＼ ／＼ W0/4 X3 F3 F6
! ─┴─────────①╋─────┴───①╋─────┬─╋────
! ／ ＼ ／ ＼ ＼／
! X7 ／ ＼ W3/8 X7 ／ ＼ W1/4 X7 ／＼ W0/2 F7 F7
! ─┴─────────①╋─────┴───①╋─────┴①╋────
!┠←─ ネスト(Nested)された時間域入力の２分割の重ね ─→╂基数２のＤＦＴ─┨
!
!
! 計算プログラムは、和と差の１対 の繰返しであり、同じ処理で実行するには、
! 各段の計算結果が、次段の計算の入力にもなるため、同じ場所に書き戻す必要があった。
!
! X0~,F0~ の添え字を、直接の物理アドレスとして、配列を組むと、書き込む前に有った
! データーが使用される前に、上書きされてしまう不都合が発生する。その解決に、
!
! X0~,F0~ の添え字を、論理アドレスとして無視、図の上から下への配置を物理アドレス
! として、配列を用いる計算方法をとる。その為には、
!
! 右の出力信号内容 F0,F4,F2,, を、最後に F0,F1,F2,, へ戻さなければならない。
! この操作を、Index Bit Reverse と呼ぶ。
!
! X0~,F0~ の添え字を、直接の物理アドレスとしても、上書されないように出来れば、
! Index Bit Reverse は、不要になるので、その例を最後に紹介する。
!
!=========================================
! 高速フーリエ変換 Fast Fourier Transform
! 周波数間引き型 DIF (Decimation In Frequency)
!
SUB FFT_DIF
!CALL ref_DIF(0,N1) !再帰型
CALL normal_DIF !通常型（高速)
!-------------
CALL IBrev !Index Bit Reverse
END SUB
!再帰型（通常型より約１.５倍に遅くなるが、簡潔。)
SUB ref_DIF(J,N)
IF 1< N THEN
LET H=N1 !FFT: N1 IFFT: 0
FOR I=J TO J+N/2-1
LET IK=I+N/2
LET ZT=Z(I)-Z(IK)
LET Z(I)=Z(I)+Z(IK) ! X(r)+X(r+N/2)
LET Z(IK)=ZT*w19(H) !( X(r)-X(r+N/2) )*Wr/N
LET H=H-N1/N !FFT: -N1/N IFFT: +N1/N
NEXT I
CALL ref_DIF(J ,N/2)
CALL ref_DIF(J+N/2,N/2)
END IF
END SUB
!----
!通常型（高速)
SUB normal_DIF
LET N2=N1
DO WHILE 1< N2
LET N11=N2
LET N2=N2/2
LET D=N1/N11
LET H=N1 !FFT: N1 IFFT: 0
FOR J=0 TO N2-1
FOR I=J TO N1-1 STEP N11
LET IK=I+N2
LET ZT=Z(I)-Z(IK)
LET Z(I)=Z(I)+Z(IK) ! X(r)+X(r+N/2)
LET Z(IK)=ZT*w19(H) !( X(r)-X(r+N/2) )*Wr/N
NEXT I
LET H=H-D !FFT: -D IFFT: +D
NEXT J
LOOP
END SUB
!Page-3 へ続く
-
- [15]
Page-3
- 投稿者：SECOND
- 投稿日：2010年 3月25日(木)05時11分12秒
-
-
!Page-3 の始め
!=========================================
! 高速フーリエ 逆変換 Inverse Fast Fourier Transform
! 周波数間引き型 DIF (Decimation In Frequency) と同文で、
! 回転子W だけを 逆回しにしたもの。内容は時間間引きになる。
!
SUB IFFT_DIF
!CALL ref_IDIF(0,N1) !再帰型
CALL normal_IDIF !通常型（高速)
!-------------
CALL IBrev !Index Bit Reverse
!-------------
MAT Z=(1/N1)*Z
END SUB
!再帰型（通常型より約１.５倍に遅くなるが、簡潔。)
SUB ref_IDIF(J,N)
IF 1< N THEN
LET H=0 !FFT: N1 IFFT: 0
FOR I=J TO J+N/2-1
LET IK=I+N/2
LET ZT=Z(I)-Z(IK)
LET Z(I)=Z(I)+Z(IK) ! X(r)+X(r+N/2)
LET Z(IK)=ZT*w19(H) !( X(r)-X(r+N/2) )*Wr/N
LET H=H+N1/N !FFT: -N1/N IFFT: +N1/N
NEXT I
CALL ref_IDIF(J ,N/2)
CALL ref_IDIF(J+N/2,N/2)
END IF
END SUB
!----
!通常型（高速)
SUB normal_IDIF
LET N2=N1
DO WHILE 1< N2
LET N11=N2
LET N2=N2/2
LET D=N1/N11
LET H=0 !FFT: N1 IFFT :0
FOR J=0 TO N2-1
FOR I=J TO N1-1 STEP N11
LET IK=I+N2
LET ZT=Z(I)-Z(IK)
LET Z(I)=Z(I)+Z(IK) ! X(r)+X(r+N/2)
LET Z(IK)=ZT*w19(H) !( X(r)-X(r+N/2) )*Wr/N
NEXT I
LET H=H+D !FFT: -D IFFT: +D
NEXT J
LOOP
END SUB
!=========================================
! INDEX BIT REVERSE を使用しない方法
!
!=========================================
! 高速フーリエ変換 Fast Fourier Transform
! 時間間引き型 DIT (Decimation In Time) without IndexBitReverse
!
SUB FFT_DIT_NB
LET I=0
LET O=MOD(LEN(BSTR$(N1,2)),2)*N1 !bitand(N1,0x5555)
LET N2=N1/2
LET K=1
LET D=N2
DO WHILE K<=N2 ! FFT: FOR H=N1 TO N1-D*(K-1) STEP -D
FOR H=N1 TO N1-D*(K-1) STEP -D !IFFT: FOR H=0 TO D*(K-1) STEP D
FOR I=I TO I+D-1
LET ZT=w19(H)*Z(I+D)
LET Z(O+N2)=Z(I)-ZT
LET Z(O) =Z(I)+ZT
LET O=O+1
NEXT I
LET I=I+D
NEXT H
LET K=K+K
LET D=D/2
IF O<=N1 THEN LET I=0 ELSE LET I=N1
LET O=N1-I
LOOP
END SUB
!=========================================
! 高速フーリエ 逆変換 Inverse Fast Fourier Transform
! 時間間引き型 DIT (Decimation In Time) without IndexBitReverse と同文で、
! 回転子W を 逆回しにしたもの。内容は周波数間引きになる。
!
SUB IFFT_DIT_NB
LET I=0
LET O=MOD(LEN(BSTR$(N1,2)),2)*N1 !bitand(N1,0x5555)
LET N2=N1/2
LET K=1
LET D=N2
DO WHILE K<=N2 ! FFT: FOR H=N1 TO N1-D*(K-1) STEP -D
FOR H=0 TO D*(K-1) STEP D !IFFT: FOR H=0 TO D*(K-1) STEP D
FOR I=I TO I+D-1
LET ZT=w19(H)*Z(I+D)
LET Z(O+N2)=Z(I)-ZT
LET Z(O) =Z(I)+ZT
LET O=O+1
NEXT I
LET I=I+D
NEXT H
LET K=K+K
LET D=D/2
IF O<=N1 THEN LET I=0 ELSE LET I=N1
LET O=N1-I
LOOP
!-------------
MAT Z=(1/N1)*Z
END SUB
!=========================================
! 高速フーリエ変換 Fast Fourier Transform
! 周波数間引き型 DIF (Decimation In Frequency) without IndexBitReverse
!
SUB FFT_DIF_NB
LET I=0
LET O=N1
LET N2=N1/2
LET K=1
LET D=N2
DO WHILE K<=N2 ! FFT: FOR H=N1 TO N1-K*(D-1) STEP -K
FOR H=N1 TO N1-K*(D-1) STEP -K !IFFT: FOR H=0 TO K*(D-1) STEP K
FOR I=I TO I+K-1
LET ZT=Z(I+N2)
LET Z(O+K)=(Z(I)-ZT)*w19(H)
LET Z(O )= Z(I)+ZT
LET O=O+1
NEXT I
LET O=O+K
NEXT H
LET K=K+K
LET D=D/2
IF I>N1 THEN LET I=0 ELSE LET I=N1
IF K=N2 THEN LET O=0 ELSE LET O=N1-I
LOOP
END SUB
!=========================================
! 高速フーリエ 逆変換 Inverse Fast Fourier Transform
! 周波数間引き型 DIF (Decimation In Frequency) without IndexBitReverse と同文で、
! 回転子W を 逆回しにしたもの。内容は時間間引きになる。
!
SUB IFFT_DIF_NB
LET I=0
LET O=N1
LET N2=N1/2
LET K=1
LET D=N2
DO WHILE K<=N2 ! FFT: FOR H=N1 TO N1-K*(D-1) STEP -K
FOR H=0 TO K*(D-1) STEP K !IFFT: FOR H=0 TO K*(D-1) STEP K
FOR I=I TO I+K-1
LET ZT=Z(I+N2)
LET Z(O+K)=(Z(I)-ZT)*w19(H)
LET Z(O )= Z(I)+ZT
LET O=O+1
NEXT I
LET O=O+K
NEXT H
LET K=K+K
LET D=D/2
IF I>N1 THEN LET I=0 ELSE LET I=N1
IF K=N2 THEN LET O=0 ELSE LET O=N1-I
LOOP
!-------------
MAT Z=(1/N1)*Z
END SUB
!Page-4 へ続く
-
- [16]
Page-4
- 投稿者：SECOND
- 投稿日：2010年 3月25日(木)05時12分54秒
-
-
!Page-4 の始め
!******************************************************************
! 上記の評価用プログラム
!******************************************************************
LET ym=10
DIM XC(4),YC(ym),XG(4),YG(ym)
!
LET XCd=7 !dots.chr.Xw
LET YCd=13+1 !dots.chr.Yw+1
LET YFL=1.52 !Data max.amplitude
!
LET XC(0)=5 !left margin
LET XG(0)=XC(0)*XCd
FOR i=1 TO 4
LET XC(i)=XC(i-1)+18
LET XG(i)=XC(i)*XCd
NEXT i
LET YC(0)=YFL+1
LET YG(0)=YC(0)*YCd
FOR i=1 TO ym
LET YC(i)=YC(i-1)+YC(0)+YFL
LET YG(i)=YC(i)*YCd
NEXT i
!
LET i=(XC(4)+13)*XCd
LET j=(YC(ym)+5)*YCd
SET bitmap SIZE i,j
SET WINDOW 0,i-1,j-1,0
!
FOR GS=1 TO 10000
IF FP(N1/GS)=0 THEN
LET XD=(XG(2)-XG(0))/N1*GS !Data pitch
IF 2<=XD THEN EXIT FOR
END IF
NEXT GS
LET YD=YFL*YCd !Data max.amplitude
!
FOR test=0 TO 2
MAT Z=ZER
FOR i=0 TO test
LET Z(i)=1 !試験用 "SIGNAL" 作成。
NEXT i
!----------グリッド-----
CLEAR
SET LINE COLOR "gray"
FOR i=0 TO ym
PLOT LINES: XG(0),YG(i);XG(4),YG(i) !,,,0x1111
NEXT i
FOR i=0 TO 4
PLOT LINES: XG(i),YG(0);XG(i),YG(ym) !,,,0x1111
NEXT i
LET i=XG(0)+XCd*1.3
LET j=YG(ym)+YCd*2.5
PLOT TEXT,AT i+XCd*2,j+YCd*.5 :"実数部 虚数部 絶対値"
SET AREA COLOR 4 ! red
DRAW disk WITH SCALE(10,1)*SHIFT(i,j)
SET AREA COLOR 10 ! dark green
DRAW disk WITH SCALE(10,1)*SHIFT(i+XCd*15,j)
SET AREA COLOR 2 ! blue
DRAW disk WITH SCALE(10,1)*SHIFT(i+XCd*30,j)
!---------------------[0]
PLOT TEXT,AT XG(4)+XCd,YG(0) :"SIGNAL"
CALL SCALE
CALL DRAW_re_im(YG(0))
!CALL DRAW_abs(YG(0))
!---------------------[1] DFT は、著しく遅いため 512 までで休止
IF N1<=512 THEN !! させているが、必要なら 数値を大きく。
CALL FFT(" DFT",YG(1))
CALL SCALE
CALL DRAW_re_im(YG(1))
CALL DRAW_abs(YG(1))
!---------------------[2]
CALL FFT("IDFT",YG(2))
CALL SCALE
CALL DRAW_re_im(YG(2))
!CALL DRAW_abs(YG(2))
ELSE
PLOT TEXT,AT XG(4)+XCd,YG(1)+YCd :" DFT 休止中"
PLOT TEXT,AT XG(4)+XCd,YG(2)+YCd :"IDFT 休止中"
END IF
!---------------------[3]
CALL FFT(" FFT_DIT",YG(3))
CALL SCALE
CALL DRAW_re_im(YG(3))
CALL DRAW_abs(YG(3))
!---------------------[4]
CALL FFT("IFFT_DIT",YG(4))
CALL SCALE
CALL DRAW_re_im(YG(4))
!CALL DRAW_abs(YG(4))
!---------------------[5]
CALL FFT(" FFT_DIF",YG(5))
CALL SCALE
CALL DRAW_re_im(YG(5))
CALL DRAW_abs(YG(5))
!---------------------[6]
CALL FFT("IFFT_DIF",YG(6))
CALL SCALE
CALL DRAW_re_im(YG(6))
!CALL DRAW_abs(YG(6))
!---------------------[7]
CALL FFT(" FFT_DIT_NB",YG(7))
CALL SCALE
CALL DRAW_re_im(YG(7))
CALL DRAW_abs(YG(7))
!---------------------[8]
CALL FFT("IFFT_DIT_NB",YG(8))
CALL SCALE
CALL DRAW_re_im(YG(8))
!CALL DRAW_abs(YG(8))
!---------------------[9]
CALL FFT(" FFT_DIF_NB",YG(9))
CALL SCALE
CALL DRAW_re_im(YG(9))
CALL DRAW_abs(YG(9))
!---------------------[10]
CALL FFT("IFFT_DIF_NB",YG(10))
CALL SCALE
CALL DRAW_re_im(YG(10))
!CALL DRAW_abs(YG(10))
!---------------------
WAIT DELAY 2
NEXT test
beep
!------- full scale of X() Y()-------
SUB SCALE
LET S9=0
FOR I=0 TO N1-1
IF S9< ABS(Z(I)) THEN LET S9=ABS(Z(I))
NEXT I
IF S9=0 THEN LET S9=1E-4
END SUB
!------- draw re(),im()-----------
SUB DRAW_re_im(YG)
SET BEAM MODE "IMMORTAL" !pen-on/off, only plot lines (not_JIS)
PLOT TEXT,AT XG(2)+XCd,YG-YD ,USING"####.###":S9
!----
SET LINE COLOR 4 !red
LET X9=XG(0)
FOR I=0 TO N1+N1 STEP GS
LET Y9=YG-YD*re(Z(MOD(I,N1)))/S9
PLOT LINES: X9,Y9;
IF 5< XD THEN DRAW circle WITH SCALE(1)*SHIFT(X9,Y9)
LET X9=X9+XD
NEXT I
PLOT LINES
!----
SET LINE COLOR 10 !dark green
LET X9=XG(0)
FOR I=0 TO N1+N1 STEP GS
LET Y9=YG-YD*im(Z(MOD(I,N1)))/S9
PLOT LINES: X9,Y9;
IF 5< XD THEN DRAW circle WITH SCALE(1)*SHIFT(X9,Y9)
LET X9=X9+XD
NEXT I
PLOT LINES
END SUB
!------- draw abs()-----------
!v.line&circle
SUB DRAW_abs(YG)
SET LINE COLOR 2 !blue
LET R=INT(XD/5)
IF 3< R THEN LET R=3
LET X9=XG(0)
FOR I=0 TO N1+N1 STEP GS
LET Y9=YG-YD*ABS(Z(MOD(I,N1)))/S9
PLOT LINES: X9,YG;X9,Y9
DRAW circle WITH SCALE(R)*SHIFT(X9,Y9)
LET X9=X9+XD
NEXT I
END SUB
!--------------------------------------
SUB FFT(a$,YG)
PLOT TEXT,AT XG(4)+XCd,YG :a$
LET tm0=TIME
SELECT CASE a$
!
!---- ＤＦＴ(Discrete Fourier Transform)
CASE " DFT"
CALL DFT !離散フーリエ変換（比較用基準）
CASE "IDFT"
CALL IDFT !離散フーリエ 逆変換（比較用基準）
!
!---- ＦＦＴ(Fast Fourier Transform)
CASE " FFT_DIT"
CALL FFT_DIT !入力が、時間間引き( Decimation In Time)
CASE "IFFT_DIT"
CALL IFFT_DIT !時間間引き型での、逆変換（周波数間引き入力）
CASE " FFT_DIF"
CALL FFT_DIF !出力が、周波数間引き( Decimation In Frequency)
CASE "IFFT_DIF"
CALL IFFT_DIF !周波数間引き型での、逆変換（時間間引き出力）
!
!---- 以下は、index bit reverse を、不要にした同型。
CASE " FFT_DIT_NB"
CALL FFT_DIT_NB
CASE "IFFT_DIT_NB"
CALL IFFT_DIT_NB
CASE " FFT_DIF_NB"
CALL FFT_DIF_NB
CASE "IFFT_DIF_NB"
CALL IFFT_DIF_NB
CASE ELSE
END SELECT
PLOT TEXT,AT XG(4)+XCd,YG+YCd,USING"#####ms":(TIME-tm0)*1000
END SUB
!=========================================
! Discrete Fourier Transform
! 離散フーリエ変換（比較用基準）
!
SUB DFT
MAT ZO=ZER
FOR K=0 TO N1-1
LET H=N1
FOR I=0 TO N1-1
LET ZO(K)=ZO(K)+w19(H)*Z(I)
LET H=MOD( H-K,N1)
NEXT I
NEXT K
MAT Z=ZO
END SUB
!=========================================
! Inverse Discrete Fourier Transform (standard)
! 離散フーリエ逆変換（比較用基準）
!
SUB IDFT
MAT ZO=ZER
FOR K=0 TO N1-1
LET H=0
FOR I=0 TO N1-1
LET ZO(I)=ZO(I)+w19(H)*Z(K)
LET H=MOD( H+K,N1)
NEXT I
NEXT K
MAT Z=(1/N1)*ZO
END SUB
END
戻る