• [0]
  • 複数ページ長編プログラム(新規投稿)

  • 投稿者:---
  • 投稿日:2009年12月14日(月)20時32分6秒
 
◆「スレッド管理メニュー」を開ける際の、
 スレッド管理画面パスワードは「 BASIC 」。スレッド作成パスワード と同じです。

◆投稿内容の初期フォントが、MS ゴシックでなく、MS Pゴシック になっていて、
 桁位置が崩れます。冒頭にタグ <font face="MS ゴシック"> を置くとよいです。

[1] 改訂 十進 BASIC による プログレッシブ JPG の展開と画像化。
[5] 改訂 十進 BASIC で、JPG ファイルを作る。
[9] 簡易 JPGファイル解析ツール
[11] GIF ファイルの解析ツール Ver.22(画像付)
[13] FFT プログラム

  • [1]
  • 改訂 十進 BASIC による プログレッシブ JPG の展開と画像化。

  • 投稿者:SECOND
  • 投稿日:2009年12月14日(月)20時36分38秒
  • 編集済
  • 返信
 

!改訂 十進 BASIC による プログレッシブ JPG の展開と画像化。
!
!先投稿では、プログレッシブ JPG 再生過程の画像を、何枚も表示していたが、
!途中画像は、最初のDC成分1枚だけ と、最終完成画、全2枚とした。
!Baseline JPG は、全1枚なので、画数で両者を区別できる。
!(必要なら、再生過程 全ての画像も、表示できるよう、SUB IZZRL0 に注釈行がある)
!
!色差成分 Cb Cr の間引き走査復元の塗潰しを、SUB IDDCT8X8 に統合し、簡素化した。
!その他、各所の手入れ多数。
!
!具体的、可視的なプログラムで、実行し画像化しますので、詳細事項の追跡と御参考に。
!再生できるファイルは、1000x1000 までの JPG だけで、
! baseline , spectral selection , successive approximation の3種類( web 上の、ほぼ全種)
! 参考資料:http://www.w3.org/Graphics/JPEG/itu-t81.pdf
!
!1)successive approximation AC.subsequence(Y,Cb,Cr 別々、1bit づつの処理)
!
!      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)
!
!2)successive approximation DC.subsequence(Y,Cb,Cr 別々、1bit づつの処理)
!
!  ハフマン・コード RRRRssss 部は、存在せず、
!  頭からの bit_stream.で、1bit づつ、全てのblock の DC係数 に加える。
!
!  AC と同様、0 → 無変化。1 → ±符号は上位桁に合せて加算。(絶対値が+1)
!
!※上記、successive approximation AC, DC とも、加える1は、
!  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 ファイルを作る。
!
!ハフマン符号木に、1つ空席を残さないと、開けないビューアの存在や、
!Linux 対策などが、先の投稿で、あちこち散乱していたのを、1つにし、
!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数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数: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数: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つ使わない枝を設けると、その分、符号木の
!枝ぶりが増え、全体のビット長も増え、圧縮コードとしては良くないですが、
!一部のビューアが、それを必要とする。(恐らくバグだと思います。)
!
!下の文中、" !←(+1) 符号木の最下に、空席を1つ作る。"   …の行、2つ。
!        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) 符号木の最下に、空席を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) 符号木の最下に、空席を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)  !次の頻度最小の組探しは、2分岐合計を1つにし、
      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     !両者何れよりも1つ上の節点レベル(Le0,Ad1)に、
      !---
      LET Tr(Le0,Ad1,0)=Le1        !分岐先( 節点レベル,アドレス)として2組記入
      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ファイル解析ツール
!
!表示される記号で、( ? | ? ) のように、真中を縦線で区切った( ) は、
!全体が1バイトで、左右 各々 4bit の数です。
!
!例、ハフマンテーブルの先頭1バイトの注釈は、(DC=0 AC=1| HT_0~3) 、
!    上位4bit が、DC/AC の区別、下位4bit は、標識番号の意味です。
!
!オフセットアドレスは、マーカー先頭にのみ、offset: 0x00… で、表示されます。
!左端の数値は、認識結果で区切られたデーター内容で、オフセットでありません。
!
!すき間の無い4桁数値は、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

!関連投稿(1ページ未満)
![メインスレッド458] GIF アニメ-ション を作る。
![メインスレッド460] LZW エンコーダーと、デコーダー

  • [13]
  •  FFT プログラム

  • 投稿者:SECOND
  • 投稿日:2010年 3月25日(木)05時05分48秒
  • 編集済
  • 返信
 
DEBUG ON
!******************************************************************
!  FFT プログラム
!******************************************************************
!------------------------------------------------------------------------------
!    デジタルでのフーリエ変換。          逆変換。
!
! ◆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 が N個の周期波を扱う場合は、その高速算法がある。
!  Cooley-Tukey(クーリー・テューキー)型  Fast Fourier Transform
!
!  入力信号の番号を、奇数 偶数の 時間間引き  2分割にする方法と、
!  出力信号の番号を、奇数 偶数の 周波数間引き2分割にする方法の、2通り。
!
!  周期 N の DFT 計算量は、N*N 個の積和であるが、N/2 個の DFT 2つに
!  分ける事が出来れば、(N/2)*(N/2) 個の積和2つになって、ざっとみても
!  1回で、半分に計算量が減る。これを、限界まで進めたもの。
!  ( 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
!
! 入力番号を、奇数 偶数の 時間間引き2分割、出力番号を、連続の周波数2分割
! とした 2つの 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 も、さらに2分割し、入れ子構造にしていく。
!
! 2分割を、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 の場合の例。 数式は右端から、(← 方向)へ2分割を繰返すが、
!                            実計算は左端から、(→ 方向)へ
!  ※ 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
!         ────┴①╋─────┴───①╋─────┴─────────①╋─
!      ┠基数2のDFT─╂←─ ネスト(Nested)された周波数出力の2分割の重ね ─→┨
!
!
!  計算プログラムは、和と差の1対 の繰返しであり、同じ処理で実行するには、
!  各段の計算結果が、次段の計算の入力にもなるため、同じ場所に書き戻す必要があった。
!
!  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..ビット逆順のUPカウンター
      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

!再帰型(通常型より約1.5倍に遅くなるが、簡潔。)
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

!再帰型(通常型より約1.5倍に遅くなるが、簡潔。)
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
!
! 入力番号を、連続の時間2分割、出力番号を、奇数 偶数の 周波数間引き2分割
! とした 2つのDFTが、合計で、元のDFTと同じになる様に、書くと、
!
!                                                   ※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 出力 2つになっている。
! この DFT も、さらに2分割し、入れ子構造にしていく。
!
! 2分割を、N(基数)=2 、和入力、差入力が、各々1つになるまで繰返すと、
! 下の様になる。 ※添え字は、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 の場合の例。 数式は左端から、(→ 方向)へ2分割を繰返す。
!                            実計算も左端から、(→ 方向)へ
!  ※ 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)された時間域入力の2分割の重ね ─→╂基数2のDFT─┨
!
!
!  計算プログラムは、和と差の1対 の繰返しであり、同じ処理で実行するには、
!  各段の計算結果が、次段の計算の入力にもなるため、同じ場所に書き戻す必要があった。
!
!  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

!再帰型(通常型より約1.5倍に遅くなるが、簡潔。)
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

!再帰型(通常型より約1.5倍に遅くなるが、簡潔。)
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$
   !
   !---- DFT(Discrete Fourier Transform)
   CASE " DFT"
      CALL DFT  !離散フーリエ変換(比較用基準)
   CASE "IDFT"
      CALL IDFT !離散フーリエ 逆変換(比較用基準)
      !
      !---- FFT(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



 <思いやりのあるコミュニティ宣言>
 teacup.掲示板は、皆様の権利を守りながら、思いやり、温かみのあるコミュニティづくりを応援します。
 いつもご協力いただきありがとうございます。

投稿者
題名
*内容 入力補助画像・ファイル<IMG>タグが利用可能です。(詳細)
URL
sage
レンタル掲示板