Re: 多次元配列 - nagram

2025/01/16 (Thu) 22:13:01

一般化してみました.

既存の配列 A を, bound(m,2) で下限,上限を設定し m次元の仮想配列に再定義する.
仮想配列の次元は何次元でもよい. Aの次元より大きくても小さくてもよい. (次元が同じ場合は MAT REDIM と同等)
仮想配列の全要素数は, Aの全要素数と等しくなくてもよい. (超えた場合はエラー)

◎ 外部関数 remat* … subscript の添字で指定した要素が仮想配列の何番目にあたるかを算出し,
同順位の A の要素を戻り値とする. (remat* の * は 1～4 の数字)
・ remat1(A,bound,subscript)
・ remat2(A,bound,subscript)
・ remat3(A,bound,subscript)
・ remat4(A,bound,subscript)
A … 元の配列. remat に続く数字の次元を持つ.
bound … 2次元配列 bound(m,2). m次元の仮想配列の各次元の下限,上限.
subscript … 1次元配列 subscript(m). 値を求めたい仮想配列の 1～m次の添字を指定.
配列要素の順位は次元の下位から変化する. DIM A(2,2,3) と定義された配列の順位は,
A(1,1,1),A(1,1,2),A(1,1,3),A(1,2,1),A(1,2,2),A(1,2,3),A(2,1,1),A(2,1,2),A(2,1,3),A(2,2,1),A(2,2,2),A(2,2,3)

◎ 外部副プログラム remat_print* … 仮想配列を一括出力する. (remat_print* の * は 1～4 の数字)
・ remat_print1(A,bound,width,intvl)
・ remat_print2(A,bound,width,intvl)
・ remat_print3(A,bound,width,intvl)
・ remat_print4(A,bound,width,intvl)
A … 元の配列. remat_print に続く数字の次元を持つ.
bound … 2次元配列 bound(m,2). m次元の仮想配列の各次元の下限,上限.
width … 仮想配列の出力欄の幅の仕様. width=1 は MAT PRINT V; に, width=2 は MAT PRINT V に相当.
intvl … 3次元以上(m>=3)の仮想配列の出力で, 間隔をあける仕様.
intvl=1 … 間隔を1行に固定 , intvl=2 … 第 j次元内の区切りで (m-j-1)行の間隔をあける
([参考] 十進BASIC の MAT PRINT で4次元配列を出力すると, 第1次元内・第2次元内の間隔は1行固定)


DECLARE EXTERNAL FUNCTION remat1, remat2, remat3, remat4
DECLARE EXTERNAL SUB remat_print1, remat_print2, remat_print3, remat_print4, mat_init
DECLARE NUMERIC m,j

DIM X(1 TO 3, 5 TO 6, 11 TO 15, 31 TO 39) ! X は 4次元配列
CALL mat_init(X) ! 元の配列の値を決定
PRINT SIZE(X,1);SIZE(X,2);SIZE(X,3);SIZE(X,4) ;" 元の配列 X の各次元のサイズ"
MAT PRINT X;
PRINT

LET m=5 ! 仮想配列の次元
DIM virtual(m,2), index(m)
DATA 1,2, 5,7, 11,12, -8,-5, 21,25 ! 仮想配列の下限,上限 (サイズは 2,3,2,4,5)
MAT READ virtual ! DIM V(1 TO 2, 5 TO 7, 11 TO 12, -8 TO -5, 21 TO 25)
PRINT "仮想配列の下限,上限"
MAT PRINT virtual;

DATA 2,7,11,-6,24 ! 値を求めたい仮想配列の指標
MAT READ index
PRINT "仮想配列の添え字を指定"
MAT PRINT index;
PRINT "仮想配列の下限を 1とした位置"
FOR j=1 TO m
PRINT index(j)-virtual(j,1)+1;
NEXT j
PRINT
PRINT
PRINT remat4(X,virtual,index) ;" 関数 remat4 の戻り値"

PRINT
PRINT
FOR j=1 TO m
PRINT virtual(j,2)-virtual(j,1)+1;
NEXT j
PRINT " 仮想配列の各次元のサイズ"
CALL remat_print4(X,virtual,1,2) ! 仮想配列の一括出力
END


EXTERNAL FUNCTION remat1(A(),bound(,),subscript()) ! 1次元配列 A を bound で仮想配列に再定義し,
DECLARE EXTERNAL SUB element_order ! subscript で指定した値
DECLARE NUMERIC n
CALL element_order(SIZE(A),bound,subscript,n)
!PRINT " A(";STR$(n);")" ! 元の配列での位置(下限を1とする)
LET remat1=A(LBOUND(A)+n-1)
END FUNCTION

EXTERNAL FUNCTION remat2(A(,),bound(,),subscript()) ! 2次元配列 A を bound で仮想配列に再定義し,
DECLARE EXTERNAL SUB element_order ! subscript で指定した値
DECLARE NUMERIC n,sz2,a2
CALL element_order(SIZE(A),bound,subscript,n)
LET sz2=SIZE(A,2)
LET a2=MOD(n,sz2)
IF a2=0 THEN LET a2=sz2
!PRINT n;" A(";CEIL(n/sz2);",";a2;")" ! 元の配列での位置(下限を1とする)
LET remat2=A(LBOUND(A,1)+CEIL(n/sz2)-1,LBOUND(A,2)+a2-1)
END FUNCTION

EXTERNAL FUNCTION remat3(A(,,),bound(,),subscript()) ! 3次元配列 A を bound で仮想配列に再定義し,
DECLARE EXTERNAL SUB element_order ! subscript で指定した値
DECLARE NUMERIC n,sz2,sz3,a2,a3
CALL element_order(SIZE(A),bound,subscript,n)
LET sz3=SIZE(A,3)
LET a3=MOD(n,sz3)
LET sz2=SIZE(A,2)
LET a2=MOD(CEIL(n/sz3),sz2)
IF a3=0 THEN LET a3=sz3
IF a2=0 THEN LET a2=sz2
!PRINT n;" A(";CEIL(n/(sz2*sz3));",";a2;",";a3;")" ! 元の配列での位置(下限を1とする)
LET remat3=A(LBOUND(A,1)+CEIL(n/(sz2*sz3))-1,LBOUND(A,2)+a2-1,LBOUND(A,3)+a3-1)
END FUNCTION

EXTERNAL FUNCTION remat4(A(,,,),bound(,),subscript()) ! 4次元配列 A を bound で仮想配列に再定義し,
DECLARE EXTERNAL SUB element_order ! subscript で指定した値
DECLARE NUMERIC n,sz2,sz3,sz4,a2,a3,a4
CALL element_order(SIZE(A),bound,subscript,n)
LET sz4=SIZE(A,4)
LET a4=MOD(n,sz4)
LET sz3=SIZE(A,3)
LET a3=MOD(CEIL(n/sz4),sz3)
LET sz2=SIZE(A,2)
LET a2=MOD(CEIL(n/(sz3*sz4)),sz2)
IF a4=0 THEN LET a4=sz4
IF a3=0 THEN LET a3=sz3
IF a2=0 THEN LET a2=sz2
!PRINT n;" A(";CEIL(n/(sz2*sz3*sz4));",";a2;",";a3;",";a4;")" ! 元の配列での位置(下限を1とする)
LET remat4=A(LBOUND(A,1)+CEIL(n/(sz2*sz3*sz4))-1,LBOUND(A,2)+a2-1,LBOUND(A,3)+a3-1,LBOUND(A,4)+a4-1)
END FUNCTION

EXTERNAL SUB element_order(sza,bd(,),ss(),n) ! 仮想配列中のssで指定した順位 → n
DECLARE NUMERIC m,sr,sz,j
LET m=SIZE(bd,1) ! 仮想配列の次元数
IF m<>SIZE(ss) OR SIZE(bd,2)<>2 THEN CAUSE EXCEPTION 6001 ! 配列の寸法の誤り
FOR j=1 TO m
LET bd(j,1)=ROUND(bd(j,1))
LET bd(j,2)=ROUND(bd(j,2))
LET ss(j)=ROUND(ss(j))
NEXT j
LET n=ss(m)-bd(m,1)+1
LET sr=bd(m,2)-bd(m,1)+1
IF ss(m)<bd(m,1) OR ss(m)>bd(m,2) OR sr<=0 THEN CAUSE EXCEPTION 2001 ! 配列の添字が範囲外
FOR j=m-1 TO 1 STEP -1
LET sz=bd(j,2)-bd(j,1)+1
IF ss(j)<bd(j,1) OR ss(j)>bd(j,2) OR sz<=0 THEN CAUSE EXCEPTION 2001 ! 配列の添字が範囲外
LET n=n+sr*(ss(j)-bd(j,1))
LET sr=sr*sz ! sr=仮想配列の下位からj次までの要素数
NEXT j
IF sr>sza THEN CAUSE EXCEPTION 5001 ! 仮想配列の全要素数が, 元の配列Aの全要素数より多い
END SUB


EXTERNAL SUB remat_print1(A(),bound(,),width,intvl) ! 元の配列が1次元である仮想配列を一括出力
DECLARE EXTERNAL SUB elements
DECLARE NUMERIC m, m1,k, i1, j
LET m=SIZE(bound,1) ! 仮想配列の次元
DIM r(m)
CALL elements(SIZE(A),bound,m,r)
IF intvl=2 THEN LET m1=1 ELSE LET m1=m-1
LET k=0
FOR i1=LBOUND(A) TO UBOUND(A)
LET k=k+1
IF width=1 THEN PRINT A(i1); ELSE PRINT A(i1),
FOR j=m TO m1 STEP -1
IF MOD(k,r(j))=0 THEN
PRINT
IF k=r(1) THEN EXIT SUB
ELSE
EXIT FOR
END IF
NEXT j
NEXT i1
END SUB

EXTERNAL SUB remat_print2(A(,),bound(,),width,intvl) ! 元の配列が2次元である仮想配列を一括出力
DECLARE EXTERNAL SUB elements
DECLARE NUMERIC m, m1,k, i1,i2, j
LET m=SIZE(bound,1) ! 仮想配列の次元
DIM r(m)
CALL elements(SIZE(A),bound,m,r)
IF intvl=2 THEN LET m1=1 ELSE LET m1=m-1
LET k=0
FOR i1=LBOUND(A,1) TO UBOUND(A,1)
FOR i2=LBOUND(A,2) TO UBOUND(A,2)
LET k=k+1
IF width=1 THEN PRINT A(i1,i2); ELSE PRINT A(i1,i2),
FOR j=m TO m1 STEP -1
IF MOD(k,r(j))=0 THEN
PRINT
IF k=r(1) THEN EXIT SUB
ELSE
EXIT FOR
END IF
NEXT j
NEXT i2
NEXT i1
END SUB

EXTERNAL SUB remat_print3(A(,,),bound(,),width,intvl) ! 元の配列が3次元である仮想配列を一括出力
DECLARE EXTERNAL SUB elements
DECLARE NUMERIC m, m1,k, i1,i2,i3, j
LET m=SIZE(bound,1) ! 仮想配列の次元
DIM r(m)
CALL elements(SIZE(A),bound,m,r)
IF intvl=2 THEN LET m1=1 ELSE LET m1=m-1
LET k=0
FOR i1=LBOUND(A,1) TO UBOUND(A,1)
FOR i2=LBOUND(A,2) TO UBOUND(A,2)
FOR i3=LBOUND(A,3) TO UBOUND(A,3)
LET k=k+1
IF width=1 THEN PRINT A(i1,i2,i3); ELSE PRINT A(i1,i2,i3),
FOR j=m TO m1 STEP -1
IF MOD(k,r(j))=0 THEN
PRINT
IF k=r(1) THEN EXIT SUB
ELSE
EXIT FOR
END IF
NEXT j
NEXT i3
NEXT i2
NEXT i1
END SUB

EXTERNAL SUB remat_print4(A(,,,),bound(,),width,intvl) ! 元の配列が4次元である仮想配列を一括出力
DECLARE EXTERNAL SUB elements
DECLARE NUMERIC m, m1,k, i1,i2,i3,i4, j
LET m=SIZE(bound,1) ! 仮想配列の次元
DIM r(m)
CALL elements(SIZE(A),bound,m,r)
IF intvl=2 THEN LET m1=1 ELSE LET m1=m-1
LET k=0
FOR i1=LBOUND(A,1) TO UBOUND(A,1)
FOR i2=LBOUND(A,2) TO UBOUND(A,2)
FOR i3=LBOUND(A,3) TO UBOUND(A,3)
FOR i4=LBOUND(A,4) TO UBOUND(A,4)
LET k=k+1
IF width=1 THEN PRINT A(i1,i2,i3,i4); ELSE PRINT A(i1,i2,i3,i4),
FOR j=m TO m1 STEP -1
IF MOD(k,r(j))=0 THEN
PRINT
IF k=r(1) THEN EXIT SUB
ELSE
EXIT FOR
END IF
NEXT j
NEXT i4
NEXT i3
NEXT i2
NEXT i1
END SUB

EXTERNAL SUB elements(sza,bb(,),m,r())
DECLARE NUMERIC j
IF SIZE(bb,2)<>2 THEN CAUSE EXCEPTION 6001 ! 配列の寸法の誤り
FOR j=1 TO m
LET bb(j,1)=ROUND(bb(j,1))
LET bb(j,2)=ROUND(bb(j,2))
NEXT j
LET r(m)=bb(m,2)-bb(m,1)+1
IF r(m)<=0 THEN CAUSE EXCEPTION 2001 ! 仮想配列の添字が範囲外
IF m>=2 THEN
FOR j=m-1 TO 1 STEP -1
LET r(j)=r(j+1)*(bb(j,2)-bb(j,1)+1)
IF r(j)<=0 THEN CAUSE EXCEPTION 2001 ! 仮想配列の添字が範囲外
NEXT j
END IF
IF r(1)>sza THEN CAUSE EXCEPTION 5001 ! 仮想配列の全要素数が, 元の配列の全要素数より多い
END SUB


EXTERNAL SUB mat_init(A(,,,)) ! 元の配列の値を決定
DECLARE NUMERIC L(4), i, i1,i2,i3,i4
FOR i=1 TO 4
LET L(i)=LBOUND(A,i)
NEXT i
FOR i1=1 TO SIZE(A,1)
FOR i2=1 TO SIZE(A,2)
FOR i3=1 TO SIZE(A,3)
FOR i4=1 TO SIZE(A,4)
LET A(L(1)+i1-1,L(2)+i2-1,L(3)+i3-1,L(4)+i4-1)=1000*i1+100*i2+10*i3+i4
! LET A(L(1)+i1-1,L(2)+i2-1,L(3)+i3-1,L(4)+i4-1)=INT(10000*RND)
NEXT i4
NEXT i3
NEXT i2
NEXT i1
END SUB

Re: 数独 - SHIRAISHI Kazuo

2025/01/08 (Wed) 08:01:25

とりあえず気になったのが，変数trueの値を変える文が見当たらないことです。十進BASICでは変数の初期値は0に設定されています。

Re: 数独 - 永野護

2025/01/08 (Wed) 14:46:19

アドバイスに感謝します。

Re: 数独 - しばっち

2025/01/12 (Sun) 07:59:52

次のプログラムはバックトラック法(後戻り法)と呼ばれる手法を
使用しています。

このバックトラック法は今回のようなパズル問題を解くのに有効な
方法ですがプログラミング初心者には少々難しいかもしれません。

この手法は覚えておいて損はないと思います。

https://ja.wikipedia.org/wiki/バックトラッキング
https://www.cc.kyoto-su.ac.jp/~yamada/ap/backtrack.html
http://www.tommylab.ynu.ac.jp/lecture/Algorithm/10/10.pdf
https://sevendays-study.com/algorithm/ex-day4.html


内部関数、内部副プログラムのままではうまくいかず
外部関数、外部副プログラムとするとうまくいったので
パックトラック法を実装する時は外部関数や外部副プログラム
として定義した方がいいかと思います。


OPTION BASE 0
PUBLIC NUMERIC Board(5, 5),True,False,COUNT

!' サンプルの6x6数独問題 (0は空きマス)
DATA 1, 0, 3, 0, 0, 0
DATA 0, 0, 0, 0, 0, 4
DATA 0, 0, 0, 0, 2, 0
DATA 0, 5, 0, 0, 0, 0
DATA 3, 0, 0, 0, 0, 0
DATA 0, 0, 0, 4, 0, 0
LET True=1
LET False=0

!' 盤面を初期化
FOR i = 0 TO 5
FOR j = 0 TO 5
READ Board(i, j)
NEXT j
NEXT i

PRINT "元の盤面:"
FOR i = 0 TO 5
FOR j = 0 TO 5
IF Board(i, j) = 0 THEN
PRINT ".";
ELSE
PRINT STR$(Board(i, j));
END IF
NEXT j
PRINT
NEXT i
CALL SolveSudoku(0)
END

EXTERNAL FUNCTION IsValid (Row, Col, Num)
!' 行と列をチェック
FOR i = 0 TO 5
IF Board(Row, i) = Num THEN
LET IsValid = False
EXIT FUNCTION
END IF
IF Board(i, Col) = Num THEN
LET IsValid = False
EXIT FUNCTION
END IF
NEXT i
LET IsValid = True
END FUNCTION

EXTERNAL SUB SolveSudoku(N)
IF N>35 THEN
LET COUNT=COUNT+1
PRINT "解けた盤面:";COUNT
FOR i = 0 TO 5
FOR j = 0 TO 5
PRINT Board(i, j);
NEXT j
PRINT
NEXT i
IF COUNT=100 THEN STOP
ELSE
LET Col=MOD(N,6)
LET Row=INT(N/6)
IF Board(Row, Col) = 0 THEN
FOR Num = 1 TO 6
IF IsValid(Row, Col, Num) = True THEN
LET Board(Row, Col) = Num
CALL SolveSudoku(N+1) ! 再帰呼び出し
LET Board(Row, Col) = 0
END IF
NEXT Num
ELSE
CALL SolveSudoku(N+1) ! 再帰呼び出し
END IF
END IF
END SUB


再帰呼び出しを使用することで多重ループ(ここではFORループ29個)を代替させることができます。
※下記プログラム参照
-----------------------------------------------------------------------------------
上記のプログラムを総当り法(※厳密には少し違う)で書き換えてみました。

6^29(36845653286788892983296)回という途方もない回数のループですが
if文によってループが制御されています。

プログラム的にやっていることは上記のプログラムと同じです。(プログラムの書き方が違うだけ!?)


OPTION BASE 0
DIM Board(5, 5)

!' サンプルの6x6数独問題 (0は空きマス)
DATA 1, 0, 3, 0, 0, 0
DATA 0, 0, 0, 0, 0, 4
DATA 0, 0, 0, 0, 2, 0
DATA 0, 5, 0, 0, 0, 0
DATA 3, 0, 0, 0, 0, 0
DATA 0, 0, 0, 4, 0, 0
LET True=1
LET False=0

!' 盤面を初期化
FOR i = 0 TO 5
FOR j = 0 TO 5
READ Board(i, j)
NEXT j
NEXT i

PRINT "元の盤面:"
FOR i = 0 TO 5
FOR j = 0 TO 5
IF Board(i, j) = 0 THEN
PRINT ".";
ELSE
PRINT STR$(Board(i, j));
END IF
NEXT j
PRINT
NEXT i

FOR A1=1 TO 6
IF IsValid (0,1,A1)=True THEN
LET Board(0,1)=A1
FOR A2=1 TO 6
IF IsValid (0, 3, A2)=True THEN
LET Board(0,3)=A2
FOR A3=1 TO 6
IF IsValid (0, 4, A3)=True THEN
LET Board(0,4)=A3
FOR A4=1 TO 6
IF IsValid (0, 5, A4)=True THEN
LET Board(0,5)=A4
FOR B1=1 TO 6
IF IsValid (1, 0, B1)=True THEN
LET Board(1,0)=B1
FOR B2=1 TO 6
IF IsValid (1, 1, B2)=True THEN
LET Board(1,1)=B2
FOR B3=1 TO 6
IF IsValid (1, 2, B3)=True THEN
LET Board(1,2)=B3
FOR B4=1 TO 6
IF IsValid (1, 3, B4)=True THEN
LET Board(1,3)=B4
FOR B5=1 TO 6
IF IsValid (1, 4, B5)=True THEN
LET Board(1,4)=B5
FOR C1=1 TO 6
IF IsValid (2, 0, C1)=True THEN
LET Board(2,0)=C1
FOR C2=1 TO 6
IF IsValid (2, 1, C2)=True THEN
LET Board(2,1)=C2
FOR C3=1 TO 6
IF IsValid (2, 2, C3)=True THEN
LET Board(2,2)=C3
FOR C4=1 TO 6
IF IsValid (2, 3, C4)=True THEN
LET Board(2,3)=C4
FOR C5=1 TO 6
IF IsValid (2, 5, C5)=True THEN
LET Board(2,5)=C5
FOR D1=1 TO 6
IF IsValid (3, 0, D1)=True THEN
LET Board(3,0)=D1
FOR D2=1 TO 6
IF IsValid (3, 2, D2)=True THEN
LET Board(3,2)=D2
FOR D3=1 TO 6
IF IsValid (3, 3, D3)=True THEN
LET Board(3,3)=D3
FOR D4=1 TO 6
IF IsValid (3, 4, D4)=True THEN
LET Board(3,4)=D4
FOR D5=1 TO 6
IF IsValid (3, 5, D5)=True THEN
LET Board(3,5)=D5
FOR E1=1 TO 6
IF IsValid (4, 1, E1)=True THEN
LET Board(4,1)=E1
FOR E2=1 TO 6
IF IsValid (4, 2, E2)=True THEN
LET Board(4,2)=E2
FOR E3=1 TO 6
IF IsValid (4, 3, E3)=True THEN
LET Board(4,3)=E3
FOR E4=1 TO 6
IF IsValid (4, 4, E4)=True THEN
LET Board(4,4)=E4
FOR E5=1 TO 6
IF IsValid (4, 5, E5)=True THEN
LET Board(4,5)=E5
FOR F1=1 TO 6
IF IsValid (5, 0, F1)=True THEN
LET Board(5,0)=F1
FOR F2=1 TO 6
IF IsValid (5, 1, F2)=True THEN
LET Board(5,1)=F2
FOR F3=1 TO 6
IF IsValid (5, 2, F3)=True THEN
LET Board(5,2)=F3
FOR F4=1 TO 6
IF IsValid (5, 4, F4)=True THEN
LET Board(5,4)=F4
FOR F5=1 TO 6
IF IsValid (5, 5, F5)=True THEN
LET Board(5,5)=F5
LET COUNT=COUNT+1
PRINT "解けた盤面:";COUNT
FOR i = 0 TO 5
FOR j = 0 TO 5
PRINT Board(i, j);
NEXT j
PRINT
NEXT i
IF COUNT=100 THEN STOP
LET Board(5,5)=0
END IF
NEXT F5
LET Board(5,4)=0
END IF
NEXT F4
LET Board(5,2)=0
END IF
NEXT F3
LET Board(5,1)=0
END IF
NEXT F2
LET Board(5,0)=0
END IF
NEXT F1
LET Board(4,5)=0
END IF
NEXT E5
LET Board(4,4)=0
END IF
NEXT E4
LET Board(4,3)=0
END IF
NEXT E3
LET Board(4,2)=0
END IF
NEXT E2
LET Board(4,1)=0
END IF
NEXT E1
LET Board(3,5)=0
END IF
NEXT D5
LET Board(3,4)=0
END IF
NEXT D4
LET Board(3,3)=0
END IF
NEXT D3
LET Board(3,2)=0
END IF
NEXT D2
LET Board(3,0)=0
END IF
NEXT D1
LET Board(2,5)=0
END IF
NEXT C5
LET Board(2,3)=0
END IF
NEXT C4
LET Board(2,2)=0
END IF
NEXT C3
LET Board(2,1)=0
END IF
NEXT C2
LET Board(2,0)=0
END IF
NEXT C1
LET Board(1,4)=0
END IF
NEXT B5
LET Board(1,3)=0
END IF
NEXT B4
LET Board(1,2)=0
END IF
NEXT B3
LET Board(1,1)=0
END IF
NEXT B2
LET Board(1,0)=0
END IF
NEXT B1
LET Board(0,5)=0
END IF
NEXT A4
LET Board(0,4)=0
END IF
NEXT A3
LET Board(0,3)=0
END IF
NEXT A2
LET Board(0,1)=0
END IF
NEXT A1

FUNCTION IsValid (Row,COL,Num)
!' 行と列をチェック
FOR i = 0 TO 5
IF Board(Row, i) = Num THEN
LET IsValid = False
EXIT FUNCTION
END IF
IF Board(i, Col) = Num THEN
LET IsValid = False
EXIT FUNCTION
END IF
NEXT i
LET IsValid = True
END FUNCTION
END

Re: SET LINE WIDTHのバグ - SHIRAISHI Kazuo

2024/12/15 (Sun) 09:36:35

JIS Full BASICのPLOT LINES文にはやっかいな規定があります。
描点（BEAM）がオフになる条件がいくつか定められていて，PLOT LINES以外の図形出力文の実行時もその一つです。
独自拡張命令の　SET BEAM MODE "IMMORTAL"　を追加すると，この規定を無視します。

SET BITMAP SIZE 400,400
SET WINDOW 0,399,399,0
SET TEXT HEIGHT 50
SET TEXT COLOR 5
SET LINE WIDTH 20
SET BEAM MODE "IMMORTAL"
FOR X=0 TO 399
PLOT POINTS: 0,0
SET LINE WIDTH 20
PLOT TEXT ,AT 100,100:"ABCDE"
PLOT LINES:X,X;
NEXT X
END

Re: Collatz problem 6*n + 4.と奇数列の関係 - Tarosa

2024/11/05 (Tue) 09:34:52

計算結果出力
4 1 4
10 3 10
16 5 16
22 7 22
28 9 28
34 11 34
40 13 40
46 15 46
52 17 52
58 19 58
64 21 64
70 23 70
76 25 76
82 27 82
88 29 88
94 31 94
100 33 100

2n の　IF MOD((s-1),3)=0 　は　6*n + 4

Re: Collatz problem 6*n + 4.と奇数列の関係 - Tarosa

2024/11/05 (Tue) 17:50:10

!Collatz problem 3
OPTION ARITHMETIC NATIVE
!OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
LET cc=1
LET c1=0
LET odd=0
LET eve=0
FOR n= 0 TO 1E8
LET a1=3*3*n+2
IF MOD(a1,2)=0 THEN
LET b1=a1/2
ELSE
LET b1=a1*3+1
END IF
!PRINT a1; b1 !1

IF MOD(b1,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF

LET a2=3*3*n+4
IF MOD(a2,2)=0 THEN
LET b2=a2/2
ELSE
LET b2=a2*3+1
END IF

!PRINT a2;b2 !2

IF MOD(b2,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF

IF MOD(n,2)=0 THEN
LET a31=3*c1+1
IF MOD(a31,2)=0 THEN
LET b31=a31/2
ELSE
LET b31=a31*3+1
END IF
!PRINT a31;b31 !3-1
LET c1=c1+1
IF MOD(b31,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF

END IF

IF MOD(n,2)=1 THEN
LET a32=18*cc-2
IF MOD(a32,2)=0 THEN
LET b32=a32/2
ELSE
LET b32=a32*3+1
END IF

!PRINT a32;b32 !3-2
LET cc=cc+1
IF MOD(b32,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF

END IF

LET a4=3*3*6*n+52
IF MOD(a4,2)=0 THEN
LET b4=a4/2
ELSE
LET b4=a4*3+1
END IF
!PRINT a4;b4 !4
IF MOD(b4,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF

NEXT n

PRINT
PRINT "eve";eve
PRINT "odd";odd

LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"

END

6n+4は全て偶数

even number/2 or odd number*3+1


(6n+4)/2 ～　3n+2 更に　3n+1 と　18n+16 その後の数列は

計算結果出力

FOR n= 0 TO 100
eve 266
odd 138

FOR n= 0 TO 1000
eve 2628
odd 1376

FOR n= 0 TO 1E4
eve 26253
odd 13751

FOR n= 0 TO 1E5
eve 262503
odd 137501

FOR n= 0 TO 1E6
eve 2625003
odd 1375001

FOR n= 0 TO 1E7
eve 26250003
odd 13750001
5.24秒
eve 26250003
odd 13750001
154.43秒　　　　1000桁モード

FOR n= 0 TO 1E8
eve 262500003
odd 137500001
51.54秒

FOR n= 0 TO 1E9
eve 2625000003
odd 1375000001

eve 2625000003
odd 1375000001
39.87秒

FOR n= 0 TO 1E10
eve 26250000003
odd 13750000001
399.59秒

規則正しく偶数が多い。事は照明可能か?

Re: Collatz problem 6*n + 4.と奇数列の関係 - Tarosa

2024/11/10 (Sun) 21:23:53

!Collatz problem 3
OPTION ARITHMETIC NATIVE
LET t0=TIME
LET cc=1
LET c1=0
LET eve=0
LET odd=0
FOR n= 0 TO 1E5
LET a1=3*3*n+2
IF MOD (a1,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF
!PRINT !1
LET a2=3*3*n+4
IF MOD (a2,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF
!PRINT !2
IF MOD(n,2)=0 THEN
LET a31=3*c1+1
IF MOD (a31,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF
!PRINT !31
LET c1=c1+1
END IF
IF MOD(n,2)=1 THEN
LET a32=18*cc-2
IF MOD (a32,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF
!PRINT !32
LET cc=cc+1
END IF
LET a4=3*3*6*n+52
IF MOD (a4,2)=0 THEN
LET eve=eve+1
ELSE
LET odd=odd+1
END IF
!PRINT !4
NEXT n

PRINT "eve=";eve
PRINT "odd=";odd
LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"

END

!Collatz problem


FOR n= 0 TO 100
eve= 278
odd= 126

FOR n= 0 TO 1000
eve= 2753
odd= 1251

FOR n= 0 TO 10000
eve= 27503
odd= 12501

FOR n= 0 TO 1E5
eve= 275003
odd= 125001

FOR n= 0 TO 1E6
eve= 2750003
odd= 1250001

FOR n= 0 TO 1E7
eve= 27500003
odd= 12500001

FOR n= 0 TO 1E8
eve= 275000003
odd= 125000001

FOR n= 0 TO 1E9
eve= 2750000003
odd= 1250000001

FOR n= 0 TO 1E10
eve= 27500000003
odd= 12500000001
348.82秒

FOR n= 0 TO 1E11
eve= 275000000003
odd= 125000000001
3458.11秒




(((6n+4)/2)even/2 or odd*3+1)even/2 or odd*3+1
2,4,1,52,11,13,16,106,20,22,4,160,29,31,34,214,38,


偶数が多い。2の倍数の次に3の倍数が多いのは自明
不思議なのは・・・

Re: 6n-1,6n+1と五角数の関係 - Tarosa

2024/10/29 (Tue) 22:03:58

!6n-1,6n+1
OPTION ARITHMETIC NATIVE
FOR n=1 TO 10000/6
LET p=SQR((n*(3*n-1)/2)*24+1)
!LET p=SQR(12*n*(3*n-1)+ 1)
LET k=6*n-1
LET p1=SQR((n*(3*n+1)/2)*24+1)
!LET p1=SQR(12*n*(3*n+1)+ 1)
LET k1=6*n+1
PRINT p;k
PRINT p1;k1
NEXT n
END


計算結果出力
5 5
7 7
11 11
13 13
17 17
19 19
23 23
25 25
29 29
31 31
35 35
37 37
41 41
43 43
47 47
49 49
53 53
55 55
59 59

確率的素数と5角数の関係

Re: 6n-1,6n+1と五角数の関係 - Tarosa

2024/10/29 (Tue) 22:08:50

!https://oeis.org/A001318
!A001318 Generalized pentagonal numbers: m*(3*m - 1)/2, m = 0, +-1, +-2, +-3, ....
OPTION ARITHMETIC NATIVE
LET t0=TIME
LET k6=1299709 !31607 !10億,31607 1億,9973
LET k2=100000 !3401 !10億,3401 1億,1229
!エラトステネスの篩

DIM P(k6)
DIM A(k2) !素数
MAT P=ZER
MAT A=ZER
LET A(1)=2
LET H1=1
FOR I=3 TO SQR(k6) STEP 2
IF P(I)=0 THEN
FOR J=I*I TO k6 STEP I
LET P(J)=1
NEXT J
END IF
NEXT I
FOR I=3 TO k6 STEP 2
IF P(I)=0 THEN
LET H1=H1+1
LET A(H1)=I
END IF
NEXT I


LET Z=10000
DIM Pn(z)
DIM Pm(z)
MAT Pn = ZER


LET cc=1
FOR n=1 TO z/2
LET pp=n*(3*n-1)/2
LET Pn(cc)=pp
LET cc=cc+1
LET p1=n*(3*n+1)/2
LET Pn(cc)=p1
LET cc=cc+1
NEXT n

LET S=0
FOR n=1 TO z
LET Pm(n)=SQR(Pn(n)*24+1)
!PRINT
next n

LET c1=1
LET cc=3
FOR n=1 TO z
LET DD=Pm(n)
FOR nn=c1 TO z
IF DD=A(nn) THEN
PRINT cc;A(nn)
LET cc=cc+1
LET c1=c1+1
EXIT FOR
END IF
NEXT nn
NEXT n

LET TM=TIME-t0
PRINT USING"####." & REPEAT$("#",2):TM;
PRINT "秒"

END

計算結果出力
3 5
4 7
5 11
6 13
7 17
8 19
9 23
10 29
11 31
12 37
13 41
14 43
15 47
16 53
17 59
18 61
19 67
20 71
21 73
22 79
23 83
24 89
25 97
26 101

素数の個数と素数　　

Re: FIRフィルタ - しばっち

2024/10/13 (Sun) 08:16:05

下記のようにすれば移動平均(LPF)から別のフィルタを生成できる。

https://vision.kuee.kyoto-u.ac.jp/lecture/dsp/pdf/2013/dsp2013_12.pdf



LET TAP=5 ! 次数
LET SAMPLINGFREQ=2000 ! サンプリング周波数
LET N=SAMPLINGFREQ
DIM IN(N),OUT(N),COEF(TAP)
FOR I=1 TO N
LET IN(I)=SIN(30*I/SAMPLINGFREQ*2*PI)*.8+SIN(800*I/SAMPLINGFREQ*2*PI)*.2
NEXT I
FOR I=1 TO TAP ! 移動平均
LET COEF(I)=1/TAP
NEXT I
LET W0=SAMPLINGFREQ/4
LET MODE=0
SELECT CASE MODE
CASE 0 ! LPF
CASE 1 ! HPF
FOR I=1 TO TAP
LET COEF(I)=(-1)^I*COEF(I)
NEXT I
CASE 2 ! BPF
FOR I=1 TO TAP
LET COEF(I)=2*COS(TAP*W0/SAMPLINGFREQ)*COEF(I)
NEXT I
CASE 3 ! BSF
LET COEF(1)=1-2*COEF(1)
FOR I=2 TO TAP
LET COEF(I)=2*COS(TAP*W0/SAMPLINGFREQ)*COEF(I)
NEXT I
END SELECT
SET WINDOW 1,N/10,-1.5,1.5
DRAW GRID(N/20,.3)
FOR I=1 TO N/10
PLOT LINES:I,IN(I);
NEXT I
PLOT LINES
CALL FIR(IN,N,COEF,TAP,OUT)
SET LINE COLOR "RED"
FOR I=1 TO N/10
PLOT LINES:I,OUT(I);
NEXT I
END

以下略

Re: FIRフィルタ - しばっち

2024/10/13 (Sun) 08:17:07

DCT(離散COS変換)によるフィルタ処理

CALL DCT(X,N,Y)を実行するとN個の配列Xを配列YにDCT変換します。
Y(0)はDC(直流成分)でY(1),Y(2),Y(3)と添え字が大きくなるとAC(交流成分)が高くなり
Y(N-1)が一番高い高周波成分になります。

つまりN=8ではY(7)が一番高い高周波成分になります。
Y(0)=0とすると直流成分遮断になります。
Y(7)=0,Y(6)=0,Y(5)=0のようにすると高周波成分遮断となります。

また、Y(0)=Y(0)*1.2のようにすると増幅となり、Y(0)=Y(0)*.8のように1より小さくすると減衰させることができます。

※NUMはNで割り切れる必要があります。MOD(NUM,N)=0


DECLARE EXTERNAL FUNCTION C
PUBLIC NUMERIC N
LET N=8 ! DCTフィルタ次数
OPTION BASE 0
LET NUM=200
DIM X(N),Y(N),OUT(NUM)
FOR I=1 TO NUM
LET OUT(I)=SIN(400/48000*I*2*PI)*.6+SIN(15000/48000*I*2*PI)*.4
NEXT I
SET WINDOW 0,NUM,-1,1
DRAW GRID(NUM/5,.2)
LET MODE=0
FOR I=1 TO NUM
PLOT LINES:I,OUT(I);
NEXT I
PLOT LINES
FOR I=1 TO NUM STEP N
FOR J=0 TO N-1
LET X(J)=OUT(I+J)
NEXT J
CALL DCT2(X,N,Y) !'離散COSINE変換
SELECT CASE MODE ! フィルタ処理
CASE 0
LET Y(0)=0
CASE 1
LET Y(N-3)=0
LET Y(N-2)=0
LET Y(N-1)=0
CASE 2
LET Y(0)=Y(0)*.7
CASE 3
LET Y(N-3)=Y(N-3)*.5
LET Y(N-2)=Y(N-2)*.5
LET Y(N-1)=Y(N-1)*.5
END SELECT
CALL DCT3(Y,N,X) !'離散COSINE逆変換
FOR J=0 TO N-1
LET OUT(I+J)=X(J)
NEXT J
NEXT I
SET LINE COLOR "RED"
FOR I=1 TO NUM
PLOT LINES:I,OUT(I);
NEXT I
END

EXTERNAL FUNCTION C(X)
IF X=0 OR X=N THEN LET C=SQR(.5) ELSE LET C=1
END FUNCTION

EXTERNAL SUB DCT2(A(),N,B())
FOR I=0 TO N-1
LET S=0
FOR K=0 TO N-1
LET S=S+A(K)*COS((2*K+1)*I*PI/2/N)
NEXT K
LET B(I)=S*SQR(2/N)*C(I)
NEXT I
END SUB

EXTERNAL SUB DCT3(A(),N,B())
FOR I=0 TO N-1
LET S=0
FOR K=0 TO N-1
LET S=S+C(K)*A(K)*COS((2*I+1)*K*PI/2/N)
NEXT K
LET B(I)=S*SQR(2/N)
NEXT I
END SUB
-----------------------------------------------------------------------------------
周波数特性調べてみた。


RANDOMIZE
DECLARE EXTERNAL FUNCTION C
LET M=12
OPTION BASE 0
LET SAMPLINGFREQ=44100
LET NUM=2^M
DIM RR(NUM),II(NUM),X(64),Y(64)
LET NN=8
DO
CLEAR
MAT RR=ZER
MAT II=ZER
FOR I=1 TO NUM
LET RR(I)=2*RND-1
NEXT I
LET MODE=1
FOR I=1 TO NUM STEP NN
FOR J=0 TO NN-1
LET X(J)=RR(I+J)
NEXT J
CALL DCT2(X,NN,Y) !'離散COSINE変換
SELECT CASE MODE
CASE 0
LET Y(0)=0
LET Y(1)=0
LET Y(2)=0
CASE 1
FOR J=1 TO NN-1
LET Y(J)=0
NEXT J
END SELECT
CALL DCT3(Y,NN,X) !'離散COSINE逆変換
FOR J=0 TO NN-1
LET RR(I+J)=X(J)
NEXT J
NEXT I
CALL FFT(M,RR,II,1)
SET WINDOW -15,SAMPLINGFREQ/2,-30,40
DRAW GRID(SAMPLINGFREQ/10,10)
LET LL=SAMPLINGFREQ/NUM
FOR I=1 TO NUM/2
IF RR(I)=0 AND II(I)=0 THEN
PLOT LINES:I*LL,-30;
ELSE
PLOT LINES:I*LL, 10*LOG10(RR(I)^2+II(I)^2);
END IF
NEXT I
WAIT DELAY 1
LET NN=NN*2
LOOP UNTIL NN>64
END

以下略

Re: FIRフィルタ - しばっち

2024/10/13 (Sun) 08:18:06

下記プログラムはDCTフィルタ次数を8として逆DCTの次数を増減させている。
これは画像処理でいう拡大・縮小処理にあたります。

つまり拡大にあたる処理では周波数は減衰し縮小では周波数が増大します。
これは全体のサイズが変化しているためです。

黒線が処理前のパワースペクトルで200Hzが高くなっています。
赤線がDCT処理後のパワースペクトルです。


DECLARE EXTERNAL FUNCTION C
PUBLIC NUMERIC NN
LET NN=8 ! DCT次数
LET M=8 ! FFT次数
OPTION BASE 0
LET NUM=2^M
LET SAMPLINGFREQ=1000 ! サンプリング周波数
DIM X(NN*2),Y(NN*2),RR(NUM*2),II(NUM*2)
DIM TR(NUM*2)
CALL SETDATA
CALL DISPLAY

CALL SETDATA
FOR I=1 TO NUM STEP NN
FOR J=0 TO NN-1
LET X(J)=RR(I+J)
NEXT J
LET MM=NN/2 ! 縮小
!LET MM=NN*2 ! 拡大
CALL DCT2(X,NN,Y) !'離散COSINE変換
CALL DCT3(Y,MM,X) !'離散COSINE逆変換
FOR J=0 TO MM-1
LET TR(JJ+J)=X(J)
NEXT J
LET JJ=JJ+MM
NEXT I
MAT RR=TR
SET LINE COLOR "RED"
CALL DISPLAY

SUB SETDATA
MAT II=ZER
FOR I=1 TO NUM
LET RR(I)=SIN(200/SAMPLINGFREQ*I*2*PI)
NEXT I
END SUB

SUB DISPLAY
MAT II=ZER
CALL FFT(M,RR,II,1)
SET WINDOW -15,SAMPLINGFREQ/2,-10,40
DRAW GRID(SAMPLINGFREQ/10,10)
LET LL=SAMPLINGFREQ/NUM
FOR I=1 TO 2^M/2
IF RR(I)=0 AND II(I)=0 THEN
PLOT LINES:I*LL,-10;
ELSE
PLOT LINES:I*LL, 10*LOG10(RR(I)^2+II(I)^2);
END IF
NEXT I
PLOT LINES
END SUB
END

以下略

Re: FIRフィルタ - しばっち

2024/10/13 (Sun) 08:19:23

ウェーブレット変換の周波数特性を調べてみた。


RANDOMIZE
OPTION BASE 0
LET M=16
LET NUM=2^M
FOR LEV=1 TO 4
LET SIZE=INT((NUM+2^LEV-1)/2^LEV)*2^LEV
DIM DAT1(SIZE),DAT2(SIZE),II(SIZE)
FOR I=0 TO NUM-1
LET DAT1(I)=2*RND-1
LET DAT2(I+1)=DAT1(I)
NEXT I
MAT II=ZER
CALL FFT(M,DAT2,II,1)
CLEAR
SET WINDOW 1,NUM/2,-50,70
DRAW GRID(NUM/10,10)
FOR I=1 TO NUM/2
IF DAT2(I)<>0 OR II(I)<>0 THEN
PLOT LINES:I,10*LOG10(DAT2(I)^2+II(I)^2);
ELSE
PLOT LINES:I,-50;
END IF
NEXT I
PLOT LINES
MAT DAT2=ZER
WAIT DELAY 1
CLEAR
CALL WAVELET1D(LEV,1,SIZE,DAT1,DAT2)
LET MODE=1
SELECT CASE MODE
CASE 0
FOR I=0 TO SIZE/2^LEV-1
LET DAT2(I)=0
NEXT I
CASE 1
FOR I=SIZE/2^LEV TO SIZE-1
LET DAT2(I)=0 ! フィルタ処理
NEXT I
END SELECT
MAT DAT1=ZER
CALL IWAVELET1D(LEV,1,SIZE,DAT2,DAT1)
FOR I=0 TO SIZE-1
LET DAT2(I+1)=DAT1(I)
NEXT I
MAT II=ZER
CALL FFT(M,DAT2,II,1)
SET WINDOW 1,NUM/2,-50,70
DRAW GRID(NUM/10,10)
FOR I=1 TO NUM/2
IF DAT2(I)<>0 OR II(I)<>0 THEN
PLOT LINES:I,10*LOG10(DAT2(I)^2+II(I)^2);
ELSE
PLOT LINES:I,-50;
END IF
NEXT I
NEXT LEV
END

EXTERNAL SUB WAVELET1D(LEV,XS,SIZE,DAT1(),DAT2())
OPTION BASE 0
DIM X(SIZE),Y(SIZE)
IF LEV>0 THEN
FOR I=0 TO SIZE-1
LET X(I)=DAT1(I+XS)
NEXT I
CALL WAVELET(SIZE,X,Y)
FOR I=0 TO SIZE-1
LET DAT2(I+XS)=Y(I)
NEXT I
FOR I=0 TO SIZE-1
LET DAT1(I+XS)=DAT2(I+XS)
NEXT I
CALL WAVELET1D(LEV-1,XS,SIZE/2,DAT1,DAT2)
CALL WAVELET1D(LEV-1,XS+SIZE/2,SIZE/2,DAT1,DAT2)
END IF
END SUB

EXTERNAL SUB IWAVELET1D(LEV,XS,SIZE,DAT1(),DAT2())
OPTION BASE 0
DIM X(SIZE),Y(SIZE)
IF LEV>0 THEN
CALL IWAVELET1D(LEV-1,XS,SIZE/2,DAT1,DAT2)
CALL IWAVELET1D(LEV-1,XS+SIZE/2,SIZE/2,DAT1,DAT2)
FOR I=0 TO SIZE-1
LET X(I)=DAT1(I+XS)
NEXT I
CALL IWAVELET(SIZE,X,Y)
FOR I=0 TO SIZE-1
LET DAT2(I+XS)=Y(I)
NEXT I
FOR I=0 TO SIZE-1
LET DAT1(I+XS)=DAT2(I+XS)
NEXT I
END IF
END SUB

EXTERNAL SUB WAVELET(SIZE,DAT1(),DAT2())
FOR I=0 TO SIZE/2-1
LET DAT2(I)=.5*DAT1(I*2)+.5*DAT1(I*2+1)
LET DAT2(I+SIZE/2)=.5*DAT1(I*2)-.5*DAT1(I*2+1)
NEXT I
END SUB

EXTERNAL SUB IWAVELET(SIZE,DAT2(),DAT1())
FOR I=0 TO SIZE/2-1
LET DAT1(I*2)=DAT2(I)+DAT2(I+SIZE/2)
LET DAT1(I*2+1)=DAT2(I)-DAT2(I+SIZE/2)
NEXT I
END SUB

以下略

Re: FIRフィルタ - しばっち

2024/10/13 (Sun) 08:20:36

アダマール変換の周波数特性調べてみた。



LET M=12
LET SAMPLINGFREQ=44100
OPTION BASE 0
RANDOMIZE
LET NUM=2^M
DIM RR(NUM),II(NUM),X(2^8),Y(2^8)
LET LL=SAMPLINGFREQ/NUM
FOR BIT=3 TO 8
CLEAR
MAT RR=ZER
MAT II=ZER
LET N=2^BIT
FOR I=1 TO NUM
LET RR(I)=2*RND-1
NEXT I
FOR I=1 TO NUM STEP N
FOR J=0 TO N-1
LET X(J)=RR(I+J)
NEXT J
CALL FASTHADAMARD(X,N,Y)
FOR J=1 TO N-1
LET K=SEQUENCY(BIT,J)
LET Y(K)=0
NEXT J
CALL FASTHADAMARD(Y,N,X)
FOR J=0 TO N-1
LET RR(I+J)=X(J)/N
NEXT J
NEXT I
CALL FFT(M,RR,II,1)
SET WINDOW 1,SAMPLINGFREQ/2,-60,40
DRAW GRID(SAMPLINGFREQ/10,10)
FOR I=1 TO NUM/2
IF RR(I)=0 AND II(I)=0 THEN
PLOT LINES:I*LL,-60;
ELSE
PLOT LINES:I*LL,10*LOG10(RR(I)^2+II(I)^2);
END IF
NEXT I
WAIT DELAY 1
NEXT BIT
END

EXTERNAL SUB FASTHADAMARD(X(),N,A()) ! 高速アダマール変換
LET BIT=LOG2(N)
FOR J=1 TO BIT
LET K=0
LET L=0
FOR I=0 TO N-1
IF I<N/2 THEN
LET A(I)=X(2*K)+X(2*K+1)
LET K=K+1
ELSE
LET A(I)=X(2*L)-X(2*L+1)
LET L=L+1
END IF
NEXT I
MAT X=A
NEXT J
END SUB

EXTERNAL FUNCTION SEQUENCY(BIT,N)
FOR I=1 TO BIT
LET G=BITAND(N,1)
LET N=INT(N/2)
LET L=BITXOR(G,B)
LET A=A*2+L
LET B=L
NEXT I
LET SEQUENCY=A
END FUNCTION

Re: FIRフィルタ - しばっち

2024/10/13 (Sun) 08:21:53

位相シフト

位相を任意のシフト幅でずらします。


LET SAMPLINGFREQ=50 ! サンプリング周波数
SET WINDOW 1,SAMPLINGFREQ,-1,1
DRAW GRID(SAMPLINGFREQ/10,.2)
LET TH=PI/6 ! シフト幅
LET FREQ=4 ! 周波数
LET D=SAMPLINGFREQ/(4*FREQ)
DIM A(SAMPLINGFREQ+D)
FOR I=1 TO SAMPLINGFREQ+D
LET A(I)=SIN(FREQ/SAMPLINGFREQ*I*2*PI)
PLOT LINES:I,A(I); ! 元波形
NEXT I
PLOT LINES
SET LINE COLOR "BLUE"
FOR I=1 TO SAMPLINGFREQ
PLOT LINES:I,SIN(FREQ/SAMPLINGFREQ*I*2*PI+TH); !位相シフトした波形
NEXT I
PLOT LINES
SET LINE COLOR "RED"
FOR I=1 TO SAMPLINGFREQ
! PLOT LINES:I,SIN(FREQ/SAMPLINGFREQ*I*2*PI)*COS(TH)+COS(FREQ/SAMPLINGFREQ*I*2*PI)*SIN(TH); ! 加法定理 SIN(A+B)=SIN(A)*COS(B)+COS(A)*SIN(B)
PLOT LINES:I,A(I)*COS(TH)+A(I+D)*SIN(TH);
NEXT I
END

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:02:05

ナイキスト周波数

https://ja.wikipedia.org/wiki/ナイキスト周波数

ナイキスト周波数(サンプリング周波数の半分)を超えると下記のようになります。



LET SAMPLINGFREQ=8000 ! サンプリング周波数
LET FREQ1=1000
LET FREQ2=SAMPLINGFREQ-FREQ1 ! SAMPLINGFREQ/2<FREQ2
LET NUM=30
SET WINDOW 1,NUM,-1,1
DRAW GRID(NUM/5,.2)
FOR I=1 TO NUM
PLOT LINES:I,SIN(FREQ1/SAMPLINGFREQ*I*2*PI);
NEXT I
PLOT LINES
SET LINE COLOR "RED"
FOR I=1 TO NUM
PLOT LINES:I,SIN(FREQ2/SAMPLINGFREQ*I*2*PI);
NEXT I
END

位相がずれていますが、周期は同じになります。

位相のずれは画像では表示位置のずれとして認知できますが
音の位相ずれは聞き分けられません。


LET SAMPLINGFREQ=8000 ! サンプリング周波数
LET FREQ1=1000
LET FREQ2=SAMPLINGFREQ-FREQ1 ! SAMPLINGFREQ/2<FREQ2
LET NUM=SAMPLINGFREQ*3
LET VOL=.5
DIM LOUT(NUM)
FOR I=1 TO NUM
LET LOUT(I)=SIN(FREQ1/SAMPLINGFREQ*I*2*PI)*VOL
NEXT I
CALL WRITEWAV("原音.WAV",NUM,SAMPLINGFREQ,16,LOUT)
FOR I=1 TO NUM
LET LOUT(I)=SIN(FREQ2/SAMPLINGFREQ*I*2*PI)*VOL
NEXT I
CALL WRITEWAV("位相ずれ.WAV",NUM,SAMPLINGFREQ,16,LOUT)
END

以下略

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:03:18

周波数解析

FFTを使用しない簡易的な方法ですが周期性を探し出しています。


LET SAMPLINGFREQ=48000
!LET SAMPLINGFREQ=44100
LET N=20
LET SMIN=1E+10
LET FREQ=400
LET START=1
LET NUM=SAMPLINGFREQ
DIM IN(NUM)
FOR I=1 TO NUM
LET IN(I)=SIN(FREQ/SAMPLINGFREQ*I*2*PI)
NEXT I
FOR K=3 TO SAMPLINGFREQ/2-1 !ナイキスト周波数まで探索
LET S=0
FOR J=0 TO N-1
LET S=S+(IN(START+J)-IN(START+J+K))^2 ! S=0なら周期一致
NEXT J
IF S<SMIN THEN
IF SMIN<1E+10 THEN PRINT SAMPLINGFREQ/K;"Hz"
LET SMIN=S
END IF
NEXT K
END
-----------------------------------------------------------------------------
相関係数で周期性を探し出します。


LET SAMPLINGFREQ=48000
LET NUM=SAMPLINGFREQ
LET FREQ=400
LET N=50
DIM IN(NUM),A(N),B(N)
FOR I=1 TO NUM
LET IN(I)=SIN(FREQ/SAMPLINGFREQ*I*2*PI)
NEXT I
FOR K=3 TO 1000
FOR J=1 TO N
LET A(J)=IN(J)
LET B(J)=IN(J+K)
NEXT J
!LET S1=ABS(CORREL(N,A,B))
LET S2=ABS(RSQ(N,A,B))
LET S3=ABS(PEARSON(N,A,B))
LET S4=ABS(CORR(N,A,B))
IF S2=1 OR S3=1 OR S4=1 THEN PRINT SAMPLINGFREQ/K;"Hz"
NEXT K
END

EXTERNAL FUNCTION CORREL(N,X(),Y())
!'相関係数 -1<=R<=1
LET COV=COVARIANCE(N,X,Y)
LET PX=STDEV(N,X)
LET PY=STDEV(N,Y)
LET CORREL=COV/(PX*PY)
END FUNCTION

EXTERNAL FUNCTION RSQ(N,A(),B())
!'相関係数 0<=R<=1
LET RSQ=PEARSON(N,A,B)^2
END FUNCTION

EXTERNAL FUNCTION PEARSON(N,A(),B())
!'相関係数 -1<=R<=1
FOR I=1 TO N
LET X=X+A(I)
LET XX=XX+A(I)*A(I)
LET Y=Y+B(I)
LET YY=YY+B(I)*B(I)
LET XY=XY+A(I)*B(I)
NEXT I
LET PEARSON=(N*XY-X*Y)/SQR((N*XX-X*X)*(N*YY-Y*Y))
END FUNCTION

EXTERNAL FUNCTION CORR(N,A(),B())
!'相関係数 -1<=R<=1
FOR I=1 TO N
LET XX=XX+A(I)*A(I)
LET YY=YY+B(I)*B(I)
LET XY=XY+A(I)*B(I)
NEXT I
LET CORR=XY/SQR(XX)/SQR(YY)
END FUNCTION

EXTERNAL FUNCTION COVARIANCE(N,X(),Y())
!'共分散
LET NX=MEAN(N,X)
LET NY=MEAN(N,Y)
FOR K=1 TO N
LET S=S+(X(K)-NX)*(Y(K)-NY)
NEXT K
LET COVARIANCE=S/N
END FUNCTION

EXTERNAL FUNCTION STDEV(N,A())
!'標準偏差
FOR I=1 TO N
LET X=X+A(I)
LET XX=XX+A(I)*A(I)
NEXT I
LET STDEV=SQR((N*XX-X*X)/N/(N-1))
END FUNCTION

EXTERNAL FUNCTION MEAN(N,A())
!'相加平均
FOR I=1 TO N
LET X=X+A(I)
NEXT I
LET MEAN=X/N
END FUNCTION

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:04:26

通常FFTを行う場合は実数値なので虚部を0として処理しますが
偶数、奇数に分けることでFFTを半分のサンプル数で行えます。
メモリーの節約にもなります。


LET M=8
LET N=2^M
DIM RR(N),II(N)
FOR I=1 TO N
LET RR(I)=I ! 通常虚部は0とする
LET II(I)=0
NEXT I
CALL FFT(M,RR,II,1)
CALL FFT(M,RR,II,-1)
FOR I=1 TO N
PRINT RR(I);II(I)
NEXT I
MAT RR=ZER
MAT II=ZER
LET M=M-1
LET N=2^M
FOR I=1 TO N
LET RR(I)=2*I-1 ! 偶数、奇数に分ける
LET II(I)=2*I
NEXT I
CALL FFT(M,RR,II,1)
CALL FFT(M,RR,II,-1)
FOR I=1 TO N
PRINT RR(I);II(I)
NEXT I
END

以下略
---------------------------------------------------------------------
LET SAMPLINGFREQ=44100
LET M=10
LET N=2^M
LET LL=SAMPLINGFREQ/N
DIM RR(N),II(N)
FOR I=1 TO N
LET RR(I)=SIN(400/SAMPLINGFREQ*I*2*PI)*.5+SIN(5000/SAMPLINGFREQ*I*2*PI)*.5
LET II(I)=0 ! 通常虚部は0とする
NEXT I
CALL FFT(M,RR,II,1)
SET WINDOW 1,SAMPLINGFREQ/2,-40,70
DRAW GRID(SAMPLINGFREQ/10,10)
FOR I=1 TO N/2
PLOT LINES:I*LL,10*LOG10(RR(I)^2+II(I)^2);
NEXT I
PLOT LINES
MAT RR=ZER
MAT II=ZER
LET N=2^(M-1)
LET LL=SAMPLINGFREQ/N
FOR I=1 TO N
LET RR(I)=SIN(400/SAMPLINGFREQ*(I-1)*2*PI)*.5+SIN(5000/SAMPLINGFREQ*(I-1)*2*PI)*.5
LET II(I)=SIN(400/SAMPLINGFREQ*I*2*PI)*.5+SIN(5000/SAMPLINGFREQ*I*2*PI)*.5 ! 偶数、奇数に分けることで半分のサンプル数でできる
NEXT I
CALL FFT(M-1,RR,II,1)
SET LINE COLOR "RED"
FOR I=1 TO N/2
PLOT LINES:I*LL,10*LOG10(RR(I)^2+II(I)^2);
NEXT I
END

以下略

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:05:53

窓関数

https://ja.wikipedia.org/wiki/窓関数

窓関数はサンプル画像のような山型をしています。
正規化するとグラフの赤線のようになるようです。

https://watlab-blog.com/2019/04/20/window-correction/


LET N=128
FOR X=0 TO N
LET S=S+W(X,N)
NEXT X
LET ACF=1/(S/N)
SET WINDOW -N/20,N,-.1,ACF
DRAW GRID(N/10,.5)
FOR X=0 TO N
PLOT LINES:X,W(X,N);
NEXT X
PLOT LINES
SET LINE COLOR "RED"
FOR X=0 TO N
PLOT LINES:X,W(X,N)*ACF;
NEXT X
PRINT ACF
END

EXTERNAL FUNCTION W(X,N)
LET W=HANNING(X,N)
!LET W=TRIANGULAR(X,N)
!LET W=HAMMING(X,N)
!LET W=BARTLETT(X,N)
!LET W=BARTLETTHANN(X,N)
!LET W=NUTTALL(X,N)
!LET W=BLACKMANHARRIS(X,N)
!LET W=BLACKMANNUTTALL(X,N)
!LET W=FLATTOP(X,N)
!LET W=PARZAN(X)
!LET W=AKAIKE(X,N)
!LET W=WELCH(X,N)
!LET W=SINE(X,N)
!LET W=VORBIS(X,N)
END FUNCTION

EXTERNAL FUNCTION TRIANGULAR(X,N)
LET TRIANGULAR=1-2*ABS((2*X-N)/(2*N))
END FUNCTION

EXTERNAL FUNCTION HANNING(X,N)
LET XX=X/N
LET HANNING=0.5-0.5*COS(2*PI*XX)
END FUNCTION

EXTERNAL FUNCTION HAMMING(X,N)
LET XX=X/N
LET HAMMING=0.54-0.46*COS(2*PI*XX)
END FUNCTION

EXTERNAL FUNCTION BARTLETT(X,N)
LET XX=X/N
LET BARTLETT=1-2*ABS(XX-0.5)
END FUNCTION

EXTERNAL FUNCTION BARTLETTHANN(X,N)
LET XX=X/N
LET BARTLETTHANN=0.62-0.48*ABS(XX-0.5)-0.38*COS(2*PI*XX)
END FUNCTION

EXTERNAL FUNCTION NUTTALL(X,N)
LET XX=X/N
LET NUTTALL=0.355768-0.487396*COS(2*PI*XX)+0.144232*COS(4*PI*XX)-0.012604*COS(6*PI*XX)
END FUNCTION

EXTERNAL FUNCTION BLACKMANHARRIS(X,N)
LET XX=X/N
LET BLACKMANHARRIS=0.35875-0.48829*COS(2*PI*XX)+0.144128*COS(4*PI*XX)-0.01168*COS(6*PI*XX)
END FUNCTION

EXTERNAL FUNCTION BLACKMANNUTTALL(X,N)
LET XX=X/N
LET BLACKMANNUTTALL=0.3635819-0.4891775*COS(2*PI*XX)+0.1365995*COS(4*PI*XX)-0.0106411*COS(6*PI*XX)
END FUNCTION

EXTERNAL FUNCTION FLATTOP(X,N)
LET XX=X/N
LET FLATTOP=1-1.93*COS(2*PI*XX)+1.29*COS(4*PI*XX)+0.388*COS(6*PI*XX)+0.032*COS(8*PI*XX)
END FUNCTION

EXTERNAL FUNCTION PARZAN(X)
IF ABS(X)<=1 THEN
LET PARZAN=1-1.5*X^2+0.75*ABS(X^3)
ELSEIF ABS(X)>1 AND ABS(X)<=2 THEN
LET PARZAN=0.25*(2-ABS(X))^3
END IF
END FUNCTION

EXTERNAL FUNCTION AKAIKE(X,N)
LET XX=X/N
LET AKAIKE=0.625-0.5*COS(2*PI*XX)-0.125*COS(4*PI*XX)
END FUNCTION

EXTERNAL FUNCTION WELCH(X,N)
LET XX=X/N
LET WELCH=4*XX*(1-XX)
END FUNCTION

EXTERNAL FUNCTION SINE(X,N)
LET XX=X/N
LET SINE=SIN(PI*XX)
END FUNCTION

EXTERNAL FUNCTION VORBIS(X,N)
LET XX=X/N
LET VORBIS=SIN(PI/2*SIN(PI*XX)^2)
END FUNCTION
---------------------------------------------------------------------
LET N=1
SET WINDOW 0,N,0,1
FOR B=1 TO 32
FOR X=0 TO N STEP 1/256
LET Y=KAISER(B,X,N)
PLOT LINES:X,Y;
NEXT X
PLOT LINES
NEXT B
END

EXTERNAL FUNCTION BESSEL0(X)
!'I0(X)=Σ(X^M/2^M/M!)^2
LET A=1
LET S=1
FOR M=1 TO 1000
LET A=A*X/M/2
LET S=S+A*A
IF A*A<1E-8 THEN EXIT FOR
NEXT M
LET BESSEL0=S
END FUNCTION

EXTERNAL FUNCTION KAISER(B,X,N) !'カイザー窓
LET KAISER=BESSEL0(B*PI*SQR(1-(2*X/N-1)^2))/BESSEL0(B*PI)
END FUNCTION
---------------------------------------------------------------------
FOR A=.7 TO .1 STEP -1/16
FOR X=0 TO 1 STEP 1/64
LET Y=GAUSS(A,X,1)
PLOT LINES:X,Y;
NEXT X
PLOT LINES
NEXT A
END

EXTERNAL FUNCTION GAUSS(A,Y,N) ! 指数窓
LET X=2*Y-1
LET GAUSS=EXP(-(X/N*X/N)/A/A)
END FUNCTION

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:07:44

wavファイルを読み込んでフィルタ処理を行います。

バッファーサイズを最大2^24(32ビット版十進BASIC)として
サンプリング周波数44100Hzとすると2^24(16777216)/44100=約380秒迄の
PCMデータを読み込みできます。

また高速なFFTライブラリーの移植版を用意しました。
https://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html


サンプル数2^24個のFFT演算が500秒あまりでした。


※曲の始めは無音状態から始まり徐々に音量が増え(フェードイン)
曲の終わりは徐々に音量が減り(フェードアウト)無音状態で終わるものと想定して
窓関数は使用していません。この想定に反する場合は窓関数を使用してください。


2進モードで実行してください。


LET SIZE=2^24
DIM LIN(SIZE),RIN(SIZE)
DIM LI(SIZE),RI(SIZE)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET NN=2^INT(LOG2(NUM)+1) !'FFTポイント数
INPUT PROMPT "ローパス(1) ハイパス(2) バンドパス(3) バンドストップ(4) ":MODE
SELECT CASE MODE
CASE 1
INPUT PROMPT "通過域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
CASE 2
INPUT PROMPT "遮断域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
CASE 3
INPUT PROMPT "遮断域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
INPUT PROMPT "通過域周波数"&STR$(FREQ1)&"～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ2
CASE 4
INPUT PROMPT "通過域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
INPUT PROMPT "遮断域周波数 "&STR$(FREQ1)&"～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ2
END SELECT
CALL FFTFILTER(MODE,NN,SAMPLINGFREQ,FREQ1,FREQ2,LIN,LI)
IF CHANNEL=2 THEN CALL FFTFILTER(MODE,NN,SAMPLINGFREQ,FREQ1,FREQ2,RIN,RI)
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN)
END

EXTERNAL SUB FFTFILTER(MODE,NN,SAMPLINGFREQ,FREQ1,FREQ2,FR(),FI())
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
SELECT CASE MODE
CASE 1
LET F=FREQ1/SAMPLINGFREQ*NN
FOR I=F TO NN-F
LET FR(I)=0
LET FI(I)=0
NEXT I
CASE 2
LET F=FREQ1/SAMPLINGFREQ*NN
FOR I=1 TO F
LET FR(I)=0
LET FI(I)=0
NEXT I
FOR I=NN-F TO NN
LET FR(I)=0
LET FI(I)=0
NEXT I
CASE 3
LET FS=FREQ1/SAMPLINGFREQ*NN
LET FE=FREQ2/SAMPLINGFREQ*NN
FOR I=1 TO FS
LET FR(I)=0
LET FI(I)=0
NEXT I
FOR I=FE TO NN-FE
LET FR(I)=0
LET FI(I)=0
NEXT I
FOR I=NN-FS TO NN
LET FR(I)=0
LET FI(I)=0
NEXT I
CASE 4
LET FS=FREQ1/SAMPLINGFREQ*NN
LET FE=FREQ2/SAMPLINGFREQ*NN
FOR I=FS TO FE
LET FR(I)=0
LET FI(I)=0
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=0
LET FI(I)=0
NEXT I
END SELECT
CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),FR,FI)
END SUB

EXTERNAL SUB CDFT(N,WR,WI,AR(),AI())
OPTION BASE 0
DIM A(N)
FOR I=0 TO N/2-1
LET A(2*I)=AR(I+1)
LET A(2*I+1)=AI(I+1)
NEXT I
LET WMR=WR
LET WMI=WI
LET M=N
DO WHILE M>4
LET L=M/2
LET WKR=1
LET WKI=0
LET WDR=1-2*WMI*WMI
LET WDI=2*WMI*WMR
LET SS=2*WDI
LET WMR=WDR
LET WMI=WDI
FOR J=0 TO N-M STEP M
LET I=J+L
LET XR=A(J)-A(I)
LET XI=A(J+1)-A(I+1)
LET A(J)=A(J)+A(I)
LET A(J+1)=A(J+1)+A(I+1)
LET A(I)=XR
LET A(I+1)=XI
LET XR=A(J+2)-A(I+2)
LET XI=A(J+3)-A(I+3)
LET A(J+2)=A(J+2)+A(I+2)
LET A(J+3)=A(J+3)+A(I+3)
LET A(I+2)=WDR*XR-WDI*XI
LET A(I+3)=WDR*XI+WDI*XR
NEXT J
FOR K=4 TO L-4 STEP 4
LET WKR=WKR-SS*WDI
LET WKI=WKI+SS*WDR
LET WDR=WDR-SS*WKI
LET WDI=WDI+SS*WKR
FOR J=K TO N-M+K STEP M
LET I=J+L
LET XR=A(J)-A(I)
LET XI=A(J+1)-A(I+1)
LET A(J)=A(J)+A(I)
LET A(J+1)=A(J+1)+A(I+1)
LET A(I)=WKR*XR-WKI*XI
LET A(I+1)=WKR*XI+WKI*XR
LET XR=A(J+2)-A(I+2)
LET XI=A(J+3)-A(I+3)
LET A(J+2)=A(J+2)+A(I+2)
LET A(J+3)=A(J+3)+A(I+3)
LET A(I+2)=WDR*XR-WDI*XI
LET A(I+3)=WDR*XI+WDI*XR
NEXT J
NEXT K
LET M=L
LOOP
IF M>2 THEN
FOR J=0 TO N-4 STEP 4
LET XR=A(J)-A(J+2)
LET XI=A(J+1)-A(J+3)
LET A(J)=A(J)+A(J+2)
LET A(J+1)=A(J+1)+A(J+3)
LET A(J+2)=XR
LET A(J+3)=XI
NEXT J
END IF
IF N>4 THEN CALL BITRV2(N,A)
FOR I=0 TO N/2-1
LET AR(I+1)=A(2*I)/SQR(N)
LET AI(I+1)=A(2*I+1)/SQR(N)
NEXT I
END SUB

EXTERNAL SUB BITRV2(N,A())
LET M=N/4
LET M2=2*M
LET N2=N-2
LET K=0
FOR J=0 TO M2-4 STEP 4
IF J<K THEN
LET XR=A(J)
LET XI=A(J+1)
LET A(J)=A(K)
LET A(J+1)=A(K+1)
LET A(K)=XR
LET A(K+1)=XI
ELSEIF J>K THEN
LET J1=N2-J
LET K1=N2-K
LET XR=A(J1)
LET XI=A(J1+1)
LET A(J1)=A(K1)
LET A(J1+1)=A(K1+1)
LET A(K1)=XR
LET A(K1+1)=XI
END IF
LET K1=M2+K
LET XR=A(J+2)
LET XI=A(J+3)
LET A(J+2)=A(K1)
LET A(J+3)=A(K1+1)
LET A(K1)=XR
LET A(K1+1)=XI
LET L=M
DO WHILE K>=L
LET K=K-L
LET L=L/2
LOOP
LET K=K+L
NEXT J
END SUB

以下略
---------------------------------------------------------------------------
上記FFTフィルタ処理では全体に対してFFTしていましたが
下記プログラムでは短時間FFT(STFT)によるフィルタ処理しています。
但し、実行結果は上記プログラムとは異なります。

FFTポイント数とFFT実行回数の違いに注目してください。

https://ja.wikipedia.org/wiki/短時間フーリエ変換



LET SIZE=44100*200
DIM LIN(SIZE),RIN(SIZE),LOUT(SIZE),ROUT(SIZE)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET NN=2^10 !'FFTポイント数 短時間FFT
LET SFT=NN/8
DIM LR(NN),LI(NN),RR(NN),RI(NN)
INPUT PROMPT "ローパス(1) ハイパス(2) バンドパス(3) バンドストップ(4) ":MODE
SELECT CASE MODE
CASE 1
INPUT PROMPT "通過域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
CASE 2
INPUT PROMPT "遮断域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
CASE 3
INPUT PROMPT "遮断域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
INPUT PROMPT "通過域周波数"&STR$(FREQ1)&"～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ2
CASE 4
INPUT PROMPT "通過域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
INPUT PROMPT "遮断域周波数 "&STR$(FREQ1)&"～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ2
END SELECT
FOR K=1 TO NUM-NN STEP SFT ! オーバーラップさせる
MAT LI=ZER
MAT RI=ZER
MAT LR=ZER
MAT RR=ZER
FOR J=0 TO NN-1
IF K+J>=1 AND K+J<=NUM THEN
LET LR(J+1)=LIN(K+J)*W(J,NN) ! 窓関数をかける
LET RR(J+1)=RIN(K+J)*W(J,NN)
END IF
NEXT J
CALL FFTFILTER(MODE,NN,SAMPLINGFREQ,FREQ1,FREQ2,LR,LI)
IF CHANNEL=2 THEN CALL FFTFILTER(MODE,NN,SAMPLINGFREQ,FREQ1,FREQ2,RR,RI)
FOR J=0 TO NN-1
LET LOUT(K+J)=LOUT(K+J)+LR(J+1)/4
LET ROUT(K+J)=ROUT(K+J)+RR(J+1)/4
NEXT J
NEXT K
FOR K=1 TO NUM
LET LMAX=MAX(LMAX,ABS(LOUT(K)))
LET RMAX=MAX(RMAX,ABS(ROUT(K)))
NEXT K
MAT LOUT=(1/LMAX)*LOUT
IF CHANNEL=2 THEN MAT ROUT=(1/RMAX)*ROUT
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LOUT,ROUT)
END

以下略
---------------------------------------------------------------------------
こちらも短時間FFTによるフィルタ処理です。



LET SIZE=44100*200
DIM LIN(SIZE),RIN(SIZE),LOUT(SIZE),ROUT(SIZE)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET NN=256 !'FFTポイント数 短時間FFT
DIM LR(NN),LI(NN),RR(NN),RI(NN)
INPUT PROMPT "ローパス(1) ハイパス(2) バンドパス(3) バンドストップ(4) ":MODE
SELECT CASE MODE
CASE 1
INPUT PROMPT "通過域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
CASE 2
INPUT PROMPT "遮断域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
CASE 3
INPUT PROMPT "遮断域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
INPUT PROMPT "通過域周波数"&STR$(FREQ1)&"～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ2
CASE 4
INPUT PROMPT "通過域周波数 0 ～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ1
INPUT PROMPT "遮断域周波数 "&STR$(FREQ1)&"～"&STR$(SAMPLINGFREQ/2)&"Hz まで":FREQ2
END SELECT
FOR K=1 TO NUM
MAT LI=ZER
MAT RI=ZER
FOR J=-NN/2 TO NN/2-1
IF K+J>=1 AND K+J<=NUM THEN
LET FL=1
LET LR(NN/2+J+1)=LIN(K+J)*W(J+NN/2,NN) ! 窓関数をかける
IF CHANNEL=2 THEN LET RR(NN/2+J+1)=RIN(K+J)*W(J+NN/2,NN)
ELSE
LET FL=0
LET LOUT(K)=LIN(K)
LET ROUT(K)=RIN(K)
EXIT FOR
END IF
NEXT J
IF FL=1 THEN
CALL FFTFILTER(MODE,NN,SAMPLINGFREQ,FREQ1,FREQ2,LR,LI)
LET LOUT(K)=LR(NN/2)
IF CHANNEL=2 THEN
CALL FFTFILTER(MODE,NN,SAMPLINGFREQ,FREQ1,FREQ2,RR,RI)
LET ROUT(K)=RR(NN/2)
END IF
END IF
NEXT K
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LOUT,ROUT)
END

以下略
---------------------------------------------------------------------------
アニメーション振幅表示
但し、リアルタイム表示ではありません。


LET SIZE=44100*200
DIM LIN(SIZE),RIN(SIZE)
FILE GETOPENNAME F$,"読み込み WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET ST=SAMPLINGFREQ/30 ! 1/30秒毎
LET NN=2^INT(LOG2(ST)+1) !'FFTポイント数
LET JJ=SAMPLINGFREQ/NN ! 分解能
DIM LR(NN),RR(NN),LI(NN),RI(NN)
FOR I=1 TO NUM STEP ST
MAT LI=ZER
MAT RI=ZER
MAT LR=ZER
MAT RR=ZER
FOR J=0 TO ST-1
LET LR(J+1)=LIN(I+J)*W(J,ST) ! 窓関数を掛ける
LET RR(J+1)=RIN(I+J)*W(J,ST)
NEXT J
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),LR,LI)
IF CHANNEL=2 THEN CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),RR,RI)
SET VIEWPORT 0,1,0,.49
SET WINDOW 1,SAMPLINGFREQ/2,-50,30
DRAW GRID0(SAMPLINGFREQ/10,10)
FOR J=1 TO NN/2
IF LR(J)<>0 OR LI(J)<>0 THEN
PLOT LINES:J*JJ,-50;J*JJ,10*LOG10(LR(J)^2+LI(J)^2)
END IF
NEXT J
IF CHANNEL=2 THEN
SET VIEWPORT 0,1,.51,1
SET WINDOW 1,SAMPLINGFREQ/2,-50,30
DRAW GRID0(SAMPLINGFREQ/10,10)
FOR J=1 TO NN/2
IF RR(J)<>0 OR RI(J)<>0 THEN
PLOT LINES:J*JJ,-50;J*JJ,10*LOG10(RR(J)^2+RI(J)^2)
END IF
NEXT J
END IF
SET DRAW MODE EXPLICIT
WAIT DELAY .1
SET DRAW MODE HIDDEN
CLEAR
NEXT I
END

以下略
---------------------------------------------------------------------------
スペクトログラム表示

https://ja.wikipedia.org/wiki/スペクトログラム

カラーマップを定義してスペクトログラムを表示します。

サンプル画像参照


LET SIZE=44100*200
DIM LIN(SIZE),RIN(SIZE)
ASK COLOR MIX(15) RED,GREEN,BLUE
CALL JETCOLORMAP
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET ST=NUM/200
LET NN=2^INT(LOG2(ST)+1)
DIM LR(NN),RR(NN),LI(NN),RI(NN)
LET JJ=SAMPLINGFREQ/NN ! 分解能
LET VMIN=-70 ! db範囲
LET VMAX=0
FOR I=1 TO NUM STEP ST
MAT LI=ZER
MAT RI=ZER
MAT LR=ZER
MAT RR=ZER
FOR J=0 TO ST-1
LET LR(J+1)=LIN(I+J)*W(J,ST) ! 窓関数を掛ける
LET RR(J+1)=RIN(I+J)*W(J,ST)
NEXT J
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),LR,LI)
IF CHANNEL=2 THEN CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),RR,RI)
FOR J=1 TO NN/2-1
SET VIEWPORT 0,1,0,.49
SET WINDOW -NUM/20,NUM,-1000,SAMPLINGFREQ/2
IF LR(J)<>0 OR LI(J)<>0 THEN
LET DB=10*LOG10(LR(J)^2+LI(J)^2)
ELSE
LET DB=VMIN
END IF
LET C=INT(INTERPOLATE(DB,VMIN,VMAX,0,255))
SET AREA COLOR C
PLOT AREA:I,J*JJ;I+ST,J*JJ;I+ST,(J+1)*JJ;I,(J+1)*JJ
IF CHANNEL=2 THEN
SET VIEWPORT 0,1,.51,1
SET WINDOW -NUM/20,NUM,-1000,SAMPLINGFREQ/2
IF RR(J)<>0 OR RI(J)<>0 THEN
LET DB=10*LOG10(RR(J)^2+RI(J)^2)
ELSE
LET DB=VMIN
END IF
LET C=INT(INTERPOLATE(DB,VMIN,VMAX,0,255))
SET AREA COLOR C
PLOT AREA:I,J*JJ;I+ST,J*JJ;I+ST,(J+1)*JJ;I,(J+1)*JJ
END IF
NEXT J
NEXT I
SET COLOR MIX(15) RED,GREEN,BLUE
SET VIEWPORT 0,1,0,.49
SET WINDOW -NUM/20/SAMPLINGFREQ,NUM/SAMPLINGFREQ,-1000,SAMPLINGFREQ/2
DRAW GRID(NUM/10/SAMPLINGFREQ,5000)
SET VIEWPORT 0,1,.51,1
SET WINDOW -NUM/20/SAMPLINGFREQ,NUM/SAMPLINGFREQ,-1000,SAMPLINGFREQ/2
DRAW GRID(NUM/10/SAMPLINGFREQ,5000)
END

EXTERNAL FUNCTION INTERPOLATE(X,XMIN,XMAX,YMIN,YMAX) ! XMIN<=X<=XMAX → YMIN(0)～YMAX(255)
LET X=MIN(XMAX,MAX(X,XMIN))
LET INTERPOLATE=(X-XMIN)*(YMAX-YMIN)/(XMAX-XMIN) + YMIN
END FUNCTION

以下略
---------------------------------------------------------------------------
3バンドグラフィックイコライザー

低音域、中音域、高音域でイコライジング処理をします。

上記想定により窓関数は使用していません。


LET SIZE=2^23
DIM LIN(SIZE),RIN(SIZE)
DIM LI(SIZE),RI(SIZE)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET NN=2^INT(LOG2(NUM)+1) !'FFTポイント数
INPUT PROMPT "低音 db=":DB
LET GAIN1=10^(DB/20)
INPUT PROMPT "中音 db=":DB
LET GAIN2=10^(DB/20)
INPUT PROMPT "高音 db=":DB
LET GAIN3=10^(DB/20)
CALL FFTFILTER(NN,SAMPLINGFREQ,SAMPLINGFREQ/12,SAMPLINGFREQ*3/12,SAMPLINGFREQ*5/12,SAMPLINGFREQ/12,GAIN1,GAIN2,GAIN3,LIN,LI)
IF CHANNEL=2 THEN CALL FFTFILTER(NN,SAMPLINGFREQ,SAMPLINGFREQ/12,SAMPLINGFREQ*3/12,SAMPLINGFREQ*5/12,SAMPLINGFREQ/12,GAIN1,GAIN2,GAIN3,RIN,RI)
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN)
END

EXTERNAL SUB FFTFILTER(NN,SAMPLINGFREQ,FREQ1,FREQ2,FREQ3,WIDTH,GAIN1,GAIN2,GAIN3,FR(),FI())
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
IF FREQ1<SAMPLINGFREQ/2 THEN
LET F=FREQ1/SAMPLINGFREQ*NN
LET FS=F-WIDTH
LET FE=F+WIDTH
IF FS<0 THEN LET FS=0
IF FE>SAMPLINGFREQ/2*NN THEN LET FE=SAMPLINGFREQ/2*NN
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN1
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN1
NEXT I
END IF
IF FREQ2<SAMPLINGFREQ/2 THEN
LET F=FREQ2/SAMPLINGFREQ*NN
LET FS=F-WIDTH
LET FE=F+WIDTH
IF FS<0 THEN LET FS=0
IF FE>SAMPLINGFREQ/2*NN THEN LET FE=SAMPLINGFREQ/2*NN
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN2
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN2
NEXT I
END IF
IF FREQ3<SAMPLINGFREQ/2 THEN
LET F=FREQ3/SAMPLINGFREQ*NN
LET FS=F-WIDTH
LET FE=F+WIDTH
IF FS<0 THEN LET FS=0
IF FE>SAMPLINGFREQ/2*NN THEN LET FE=SAMPLINGFREQ/2*NN
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN3
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN3
NEXT I
END IF
CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),FR,FI)
END SUB
以下略
---------------------------------------------------------------------------
8バンドグラフィックイコライザー

100Hz,200Hz,400Hz,800Hz,1600Hz,3200Hz,6400Hz,12800Hz付近の増減によるイコライジング処理をします。

※6400Hzはサンプリング周波数16000Hz以上、12800Hzはサンプリング周波数32000Hz以上の時に有効です。

上記想定により窓関数は使用していません。


LET SIZE=2^23
DIM LIN(SIZE),RIN(SIZE)
DIM LI(SIZE),RI(SIZE)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET NN=2^INT(LOG2(NUM)+1) !'FFTポイント数
INPUT PROMPT "100Hz db=":DB ! -10db～10db
LET GAIN1=10^(DB/20)
INPUT PROMPT "200Hz db=":DB
LET GAIN2=10^(DB/20)
INPUT PROMPT "400Hz db=":DB
LET GAIN3=10^(DB/20)
INPUT PROMPT "800Hz db=":DB
LET GAIN4=10^(DB/20)
INPUT PROMPT "1600Hz db=":DB
LET GAIN5=10^(DB/20)
INPUT PROMPT "3200Hz db=":DB
LET GAIN6=10^(DB/20)
IF SAMPLINGFREQ>=16000 THEN
INPUT PROMPT "6400Hz db=":DB
LET GAIN7=10^(DB/20)
IF SAMPLINGFREQ>=32000 THEN
INPUT PROMPT "12800Hz db=":DB
LET GAIN8=10^(DB/20)
END IF
END IF
CALL FFTFILTER(NN,SAMPLINGFREQ,100,200,400,800,1600,3200,6400,12800,GAIN1,GAIN2,GAIN3,GAIN4,GAIN5,GAIN6,GAIN7,GAIN8,LIN,LI)
IF CHANNEL=2 THEN CALL FFTFILTER(NN,SAMPLINGFREQ,100,200,400,800,1600,3200,6400,12800,GAIN1,GAIN2,GAIN3,GAIN4,GAIN5,GAIN6,GAIN7,GAIN8,RIN,RI)
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN)
END

EXTERNAL SUB FFTFILTER(NN,SAMPLINGFREQ,FREQ1,FREQ2,FREQ3,FREQ4,FREQ5,FREQ6,FREQ7,FREQ8,GAIN1,GAIN2,GAIN3,GAIN4,GAIN5,GAIN6,GAIN7,GAIN8,FR(),FI())
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
IF FREQ1<SAMPLINGFREQ/2 THEN
LET F=FREQ1/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN1
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN1
NEXT I
END IF
IF FREQ2<SAMPLINGFREQ/2 THEN
LET F=FREQ2/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN2
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN2
NEXT I
END IF
IF FREQ3<SAMPLINGFREQ/2 THEN
LET F=FREQ3/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN3
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN3
NEXT I
END IF
IF FREQ4<SAMPLINGFREQ/2 THEN
LET F=FREQ4/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN4
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN4
NEXT I
END IF
IF FREQ5<SAMPLINGFREQ/2 THEN
LET F=FREQ5/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN5
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN5
NEXT I
END IF
IF FREQ6<SAMPLINGFREQ/2 THEN
LET F=FREQ6/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN6
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN6
NEXT I
END IF
IF FREQ7<SAMPLINGFREQ/2 THEN
LET F=FREQ7/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN7
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN7
NEXT I
END IF
IF FREQ8<SAMPLINGFREQ/2 THEN
LET F=FREQ8/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN8
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN8
NEXT I
END IF
CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),FR,FI)
END SUB

以下略
---------------------------------------------------------------------------
パラメトリックイコライザー

中心周波数からバンド幅域内の周波数帯を増減させます。

上記想定により窓関数は使用していません。


LET SIZE=2^23
DIM LIN(SIZE),RIN(SIZE)
DIM LI(SIZE),RI(SIZE)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET NN=2^INT(LOG2(NUM)+1) !'FFTポイント数
INPUT PROMPT "中心周波数 (0～"&STR$(SAMPLINGFREQ/2)&")=":FREQ
INPUT PROMPT "バンド幅(0～"&STR$(SAMPLINGFREQ/4)&")=":WIDTH
INPUT PROMPT "増幅率(-10～10)=":DB
LET GAIN=10^(DB/20)
CALL FFTFILTER(NN,SAMPLINGFREQ,FREQ,WIDTH,GAIN,LIN,LI)
IF CHANNEL=2 THEN CALL FFTFILTER(NN,SAMPLINGFREQ,FREQ,WIDTH,GAIN,RIN,RI)
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,16,LIN,RIN)
END

EXTERNAL SUB FFTFILTER(N,SAMPLINGFREQ,FREQ,WIDTH,GAIN,RR(),II()) ! パラメトリックEQ
CALL CDFT(2*N,COS(PI/N),SIN(PI/N),RR,II)
LET FS=(FREQ-WIDTH/2)/SAMPLINGFREQ*N
LET FE=(FREQ+WIDTH/2)/SAMPLINGFREQ*N
IF FS<0 THEN LET FS=0
IF FE>SAMPLINGFREQ/2*N THEN LET FE=SAMPLINGFREQ/2*N
FOR I=FS TO FE
LET RR(I)=RR(I)*GAIN
NEXT I
FOR I=N-FE TO N-FS
LET RR(I)=RR(I)*GAIN
NEXT I
CALL CDFT(2*N,COS(-PI/N),SIN(-PI/N),RR,II)
END SUB

以下略
---------------------------------------------------------------
16バンドグラフィックイコライザー

LET SIZE=44100*200
LET NN=1024 !'FFTポイント数
LET SFT=NN/8
LET M=16
DIM LIN(SIZE),RIN(SIZE),LOUT(SIZE),ROUT(SIZE)
DIM LR(NN),LI(NN),RR(NN),RI(NN),FREQ(M),GAIN(M)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
FOR K=1 TO M
READ FREQ(K)
IF FREQ(K)<SAMPLINGFREQ/2 THEN
INPUT PROMPT STR$(FREQ(K))&"Hz db=":DB
LET GAIN(K)=10^(DB/20)
END IF
NEXT K
FOR K=1 TO NUM-NN STEP SFT
MAT LI=ZER
MAT RI=ZER
MAT LR=ZER
MAT RR=ZER
FOR J=0 TO NN-1
LET LR(J+1)=LIN(K+J)*W(J,NN) ! 窓関数を掛ける
LET RR(J+1)=RIN(K+J)*W(J,NN)
NEXT J
CALL FFTFILTER(NN,SAMPLINGFREQ,FREQ,GAIN,M,LR,LI)
IF CHANNEL=2 THEN CALL FFTFILTER(NN,SAMPLINGFREQ,FREQ,GAIN,M,RR,RI)
FOR J=0 TO NN-1
LET LOUT(K+J)=LOUT(K+J)+LR(J+1)/4
LET ROUT(K+J)=ROUT(K+J)+RR(J+1)/4
NEXT J
NEXT K
FOR K=1 TO NUM
LET LMAX=MAX(LMAX,ABS(LOUT(K)))
LET RMAX=MAX(RMAX,ABS(ROUT(K)))
NEXT K
MAT LOUT=(1/LMAX)*LOUT
IF CHANNEL=2 THEN MAT ROUT=(1/RMAX)*ROUT
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN)
DATA 50,100,200,400,600,800,1000,1200,1600,2000,2500,3000,4000,8000,12000,16000
END

EXTERNAL SUB FFTFILTER(NN,SAMPLINGFREQ,FREQ(),GAIN(),M,FR(),FI())
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),FR,FI)
FOR K=1 TO M
IF FREQ(K)<SAMPLINGFREQ/2 THEN
LET F=FREQ(K)/SAMPLINGFREQ*NN
LET FS=F-1
LET FE=F+1
FOR I=FS TO FE
LET FR(I)=FR(I)*GAIN(K)
NEXT I
FOR I=NN-FE TO NN-FS
LET FR(I)=FR(I)*GAIN(K)
NEXT I
END IF
NEXT K
CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),FR,FI)
END SUB

以下略
---------------------------------------------------------------
ヘリウムボイス(ドナルドダック効果、ボイスチェンジャー)

https://nettyukobo.com/helium_voice/
https://ja.wikipedia.org/wiki/音速

https://www.balloonworld.jp/column/voice/



OPTION BASE 0
LET SIZE=44100*200
DIM LIN(SIZE),RIN(SIZE),LOUT(SIZE),ROUT(SIZE)
FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LIN,RIN,ERR)
IF ERR>0 THEN STOP
LET NN=2^10 !'FFTポイント数 短時間FFT
LET SFT=NN/8
LET LL=26
LET EPS=1E-15
! 音速(1気圧0度) 空気=331m/s 水素=1269m/s ヘリウム=970m/s 窒素=337m/s 酸素=317m/s 塩素=205m/s アルゴン=319m/s 二酸化炭素=258m/s
LET RATE=2 ! 空気中の音速との速度比 ヘリウム/空気=2.93 水素/空気=3.83 塩素/空気=0.61 二酸化炭素/空気=0.77
DIM LR(NN),LI(NN),RR(NN),RI(NN),LLOG(NN),RLOG(NN),LTMP(NN),RTMP(NN)
DIM LCL(NN),LCH(NN),RCL(NN),RCH(NN),L(NN),R(NN),LAMP(NN),RAMP(NN)
FOR K=1 TO NUM-NN STEP SFT
MAT LI=ZER
MAT RI=ZER
MAT LR=ZER
MAT RR=ZER
FOR J=0 TO NN-1
IF K+J>=1 AND K+J<=NUM THEN
LET LR(J+1)=LIN(K+J)*W(J,NN) ! 窓関数をかける
LET RR(J+1)=RIN(K+J)*W(J,NN)
END IF
NEXT J
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),LR,LI)
IF CHANNEL=2 THEN CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),RR,RI)
FOR I=1 TO NN
LET LLOG(I)=LOG(SQR(LR(I)^2+LI(I)^2)+EPS)
IF CHANNEL=2 THEN LET RLOG(I)=LOG(SQR(RR(I)^2+RI(I)^2)+EPS)
NEXT I
FOR I=1 TO NN
IF LR(I)=0 AND LI(I)=0 THEN LET LI(I)=0 ELSE LET LI(I)=ANGLE(LR(I),LI(I))
IF RR(I)=0 AND RI(I)=0 THEN LET RI(I)=0 ELSE LET RI(I)=ANGLE(RR(I),RI(I))
NEXT I
MAT LTMP=ZER
MAT RTMP=ZER
CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),LLOG,LTMP)
IF CHANNEL=2 THEN CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),RLOG,RTMP)
MAT LCL=LLOG
MAT LCH=LLOG
MAT RCL=RLOG
MAT RCH=RLOG
FOR I=LL TO NN-LL
LET LCL(I)=0
LET RCL(I)=0
NEXT I
FOR I=1 TO LL
LET LCH(I)=0
LET RCH(I)=0
NEXT I
FOR I=NN-LL TO NN
LET LCH(I)=0
LET RCH(I)=0
NEXT I
MAT LTMP=ZER
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),LCL,LTMP)
MAT LTMP=ZER
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),LCH,LTMP)
IF CHANNEL=2 THEN
MAT RTMP=ZER
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),RCL,RTMP)
MAT RTMP=ZER
CALL CDFT(2*NN,COS(PI/NN),SIN(PI/NN),RCH,RTMP)
END IF
MAT LAMP=(-1000)*CON
MAT RAMP=(-1000)*CON
FOR I=1 TO NN/2
LET CN=IP(I/RATE)
LET CR=FP(I/RATE)
IF CN<NN/2 THEN
LET LAMP(I)=(1-CR)*LCL(CN)+CR*LCL(CN+1)
LET RAMP(I)=(1-CR)*RCL(CN)+CR*RCL(CN+1)
END IF
NEXT I
FOR I=1 TO NN
LET L(I)=EXP(LAMP(I)+LCH(I))
LET R(I)=EXP(RAMP(I)+RCH(I))
NEXT I
FOR I=1 TO NN
LET LI(I)=L(I)*SIN(LI(I))
LET L(I)=L(I)*COS(L(I))
LET RI(I)=R(I)*SIN(RI(I))
LET R(I)=R(I)*COS(R(I))
NEXT I
CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),L,LI)
IF CHANNEL=2 THEN CALL CDFT(2*NN,COS(-PI/NN),SIN(-PI/NN),R,RI)
FOR J=0 TO NN-1
LET LOUT(K+J)=LOUT(K+J)+L(J+1)/4
LET ROUT(K+J)=ROUT(K+J)+R(J+1)/4
NEXT J
NEXT K
FOR K=1 TO NUM
LET LMAX=MAX(LMAX,ABS(LOUT(K)))
LET RMAX=MAX(RMAX,ABS(ROUT(K)))
NEXT K
MAT LOUT=(1/LMAX)*LOUT
IF CHANNEL=2 THEN MAT ROUT=(1/RMAX)*ROUT
FILE GETSAVENAME F$,"保存 WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL WRITEWAV(F$,NUM,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LOUT,ROUT)
END

以下略

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:09:40

高速なFFTライブラリーFFTWを使用したFFTルーチンも用意した。

https://www.fftw.org/

VC++2022でコンパイルしました。32bit版のみです。

実行には別途libfftw3-3.dllが必要です。

サンプル数2^24個のFFT演算が約2分程でした。

LET N=2^20
DIM RR(N),II(N)
FOR I=1 TO N
LET RR(I)=I
LET II(I)=0
NEXT I
CALL FFTW(N,RR,II,0) ! FFT
FOR I=1 TO 100
PRINT RR(I);II(I)
NEXT I
CALL FFTW(N,RR,II,1) ! IFFT
FOR I=1 TO 100
PRINT RR(I);II(I)
NEXT I
END

EXTERNAL SUB FFTW(N,RR(),II(),SW)
OPTION CHARACTER BYTE
LET RE$=REPEAT$(CHR$(0),N*8)
LET IM$=REPEAT$(CHR$(0),N*8)
FOR I=0 TO N-1
LET RE$(8*I+1:8*I+8)=PACKDBL$(RR(I+1)) !'実部
LET IM$(8*I+1:8*I+8)=PACKDBL$(II(I+1)) !'虚部
NEXT I
CALL FFTW_(N,RE$,IM$,SW)
FOR I=0 TO N-1
LET RR(I+1)=UNPACKDBL(RE$(8*I+1:8*I+8))
LET II(I+1)=UNPACKDBL(IM$(8*I+1:8*I+8))
IF SW<>0 THEN
LET RR(I+1)=RR(I+1)/N
LET II(I+1)=II(I+1)/N
END IF
NEXT I

SUB FFTW_(N,RE$,IM$,SW) ! SW=0...FFT SW<>0...IFFT
OPTION CHARACTER BYTE
ASSIGN ".\DLL\fftw.dll","fftwsub"
END SUB
END SUB
-----------------------------------------------------------------------
fftw.c


#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <fftw3.h>

__declspec(dllexport) void fftwsub(int n,double *re,double *im,int sw)
{
int i;
fftw_complex *a, *b;
a = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
b = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
if (a==NULL || b==NULL) exit(1);

fftw_plan plan;
for (i=0; i<n; i++)
{
a[i][0]=re[i];
a[i][1]=im[i];
}
if (sw==0)
{
plan = fftw_plan_dft_1d( n, a, b, FFTW_FORWARD, FFTW_ESTIMATE);
}
else
{
plan = fftw_plan_dft_1d( n, a, b, FFTW_BACKWARD, FFTW_ESTIMATE);
}
fftw_execute(plan);
for (i=0; i<n; i++) {
re[i]=b[i][0];
im[i]=b[i][1];
}
if(plan) fftw_destroy_plan(plan);
fftw_free(a);
fftw_free(b);
}

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:11:08

高速なfftライブラリーによるfftルーチン

https://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html


VC++2022でコンパイルしました。32bit版のみです。

サンプル数2^24個のFFT演算が100秒あまりでした。


LET N=2^16
DIM RR(N),II(N)
FOR I=1 TO N
LET RR(I)=I
LET II(I)=0
NEXT I
CALL CDFT(N,RR,II,1) ! FFT
FOR I=1 TO 100
PRINT RR(I);II(I)
NEXT I
CALL CDFT(N,RR,II,-1) ! IFFT
FOR I=1 TO 100
PRINT RR(I);II(I)
NEXT I
END

EXTERNAL SUB CDFT(N,RR(),II(),SW)
OPTION CHARACTER BYTE
LET A$=REPEAT$(CHR$(0),N*8*2)
LET IP$=REPEAT$(CHR$(0),(2+SQR(N))*4)
LET W$=REPEAT$(CHR$(0),N/2*8)
FOR I=0 TO N-1
LET A$(2*8*I+1:2*8*I+8)=PACKDBL$(RR(I+1)) !RE
LET A$(2*8*I+9:2*8*I+16)=PACKDBL$(II(I+1)) !IM
NEXT I
CALL CDFT_(2*N,SW,A$,IP$,W$) ! SW=1...FFT SW=-1...IFFT
FOR I=0 TO N-1
LET RR(I+1)=UNPACKDBL(A$(2*8*I+1:2*8*I+8))
LET II(I+1)=UNPACKDBL(A$(2*8*I+9:2*8*I+16))
IF SW=-1 THEN
LET RR(I+1)=RR(I+1)/N
LET II(I+1)=II(I+1)/N
END IF
NEXT I

SUB CDFT_(N,ISGN,A$,IP$,W$)
IF BITAND(N,N-1)<>0 THEN STOP
ASSIGN ".\DLL\fftsg.dll","cdft"
END SUB

!SUB CDFT_(N,ISGN,A$,IP$,W$)
! IF BITAND(N,N-1)<>0 THEN STOP
! ASSIGN ".\DLL\fft8g.dll","cdft"
!END SUB

!SUB CDFT(N,ISGN,A$,IP$,W$)
! OPTION CHARACTER BYTE
! IF BITAND(N,N-1)<>0 THEN STOP
! ASSIGN ".\DLL\fft4g.dll","cdft"
!END SUB
END SUB
----------------------------------------------------------------------------
fftsg.c


__declspec(dllexport) void cdft(int n, int isgn, double *a, int *ip, double *w)
{
void makewt(int nw, int *ip, double *w);
void cftfsub(int n, double *a, int *ip, int nw, double *w);
void cftbsub(int n, double *a, int *ip, int nw, double *w);
int nw;

nw = ip[0];
if (n > (nw << 2)) {
nw = n >> 2;
makewt(nw, ip, w);
}
if (isgn >= 0) {
cftfsub(n, a, ip, nw, w);
} else {
cftbsub(n, a, ip, nw, w);
}
}


以下略

Re: FFTフィルタ処理 - しばっち

2024/10/13 (Sun) 08:12:29

cによるDFTのDLL版です。OPENMPのマルチスレッド化による力技で高速化しています。(※CPU使用率が100%になります)
数千個程度なら個数Nに制約なしで実行できます。

VC++2022でコンパイルしました。32bit版のみです。

実行にはVCのOPEN MPランタイムvcomp140.dllが別途必要です。


EXTERNAL SUB DFT(N,RR(),II())
OPTION CHARACTER BYTE
LET AR$=REPEAT$(CHR$(0),N*8)
LET AI$=REPEAT$(CHR$(0),N*8)
FOR I=0 TO N-1
LET AR$(8*I+1:8*I+8)=PACKDBL$(RR(I+1)) !RE
LET AI$(8*I+1:8*I+8)=PACKDBL$(II(I+1)) !IM
NEXT I
CALL DFT_(N,AR$,AI$)
FOR I=0 TO N-1
LET RR(I+1)=UNPACKDBL(AR$(8*I+1:8*I+8))
LET II(I+1)=UNPACKDBL(AI$(8*I+1:8*I+8))
NEXT I

SUB DFT_(N,AR$,AI$)
OPTION CHARACTER BYTE
ASSIGN ".\DLL\dft.dll","dft"
END SUB
END SUB

EXTERNAL SUB IDFT(N,RR(),II())
OPTION CHARACTER BYTE
LET AR$=REPEAT$(CHR$(0),N*8)
LET AI$=REPEAT$(CHR$(0),N*8)
FOR I=0 TO N-1
LET AR$(8*I+1:8*I+8)=PACKDBL$(RR(I+1)) !RE
LET AI$(8*I+1:8*I+8)=PACKDBL$(II(I+1)) !IM
NEXT I
CALL IDFT_(N,AR$,AI$)
FOR I=0 TO N-1
LET RR(I+1)=UNPACKDBL(AR$(8*I+1:8*I+8))
LET II(I+1)=UNPACKDBL(AI$(8*I+1:8*I+8))
NEXT I

SUB IDFT_(N,AR$,AI$)
ASSIGN ".\DLL\dft.dll","idft"
END SUB
END SUB

!EXTERNAL SUB DFT2D(M,N,AR(,),AI(,),SW) ! SW=1...DFT SW=-1...IDFT
!OPTION CHARACTER BYTE
!LET AR$=REPEAT$(CHR$(0),M*N*8)
!LET AI$=REPEAT$(CHR$(0),M*N*8)
!FOR J=0 TO N-1
! FOR I=0 TO M-1
! LET AR$(8*(J*M+I)+1:8*(J*M+I)+8)=PACKDBL$(AR(I+1,J+1))
! LET AI$(8*(J*M+I)+1:8*(J*M+I)+8)=PACKDBL$(AI(I+1,J+1))
! NEXT I
!NEXT J
!CALL DFT2(M,N,SW,AR$,AI$)
!FOR J=0 TO N-1
! FOR I=0 TO M-1
! LET AR(I+1,J+1)=UNPACKDBL(AR$(8*(J*M+I)+1:8*(J*M+I)+8))
! LET AI(I+1,J+1)=UNPACKDBL(AI$(8*(J*M+I)+1:8*(J*M+I)+8))
! NEXT I
!NEXT J
!
!SUB DFT2(M,N,FLAG,AR$,AI$)
! ASSIGN ".\DLL\dft.dll","dft2"
!END SUB
!END SUB

!EXTERNAL SUB DFT2D(M,N,AR(,),AI(,))
!OPTION CHARACTER BYTE
!LET AR$=REPEAT$(CHR$(0),M*8)
!LET AI$=REPEAT$(CHR$(0),M*8)
!FOR J=0 TO N-1
! FOR I=0 TO M-1
! LET AR$(8*I+1:8*I+8)=PACKDBL$(AR(I,J))
! LET AI$(8*I+1:8*I+8)=PACKDBL$(AI(I,J))
! NEXT I
! CALL DFT(M,AR$,AI$)
! FOR I=0 TO M-1
! LET AR(I,J)=UNPACKDBL(AR$(8*I+1:8*I+8))
! LET AI(I,J)=UNPACKDBL(AI$(8*I+1:8*I+8))
! NEXT I
!NEXT J
!LET AR$=REPEAT$(CHR$(0),N*8)
!LET AI$=REPEAT$(CHR$(0),N*8)
!FOR J=0 TO M-1
! FOR I=0 TO N-1
! LET AR$(8*I+1:8*I+8)=PACKDBL$(AR(J,I))
! LET AI$(8*I+1:8*I+8)=PACKDBL$(AI(J,I))
! NEXT I
! CALL DFT(N,AR$,AI$)
! FOR I=0 TO N-1
! LET AR(J,I)=UNPACKDBL(AR$(8*I+1:8*I+8))
! LET AI(J,I)=UNPACKDBL(AI$(8*I+1:8*I+8))
! NEXT I
!NEXT J
!END SUB

!EXTERNAL SUB IDFT2D(M,N,AR(,),AI(,))
!OPTION CHARACTER BYTE
!LET AR$=REPEAT$(CHR$(0),M*8)
!LET AI$=REPEAT$(CHR$(0),M*8)
!FOR J=0 TO N-1
! FOR I=0 TO M-1
! LET AR$(8*I+1:8*I+8)=PACKDBL$(AR(I,J))
! LET AI$(8*I+1:8*I+8)=PACKDBL$(AI(I,J))
! NEXT I
! CALL IDFT(M,AR$,AI$)
! FOR I=0 TO M-1
! LET AR(I,J)=UNPACKDBL(AR$(8*I+1:8*I+8))
! LET AI(I,J)=UNPACKDBL(AI$(8*I+1:8*I+8))
! NEXT I
!NEXT J
!LET AR$=REPEAT$(CHR$(0),N*8)
!LET AI$=REPEAT$(CHR$(0),N*8)
!FOR J=0 TO M-1
! FOR I=0 TO N-1
! LET AR$(8*I+1:8*I+8)=PACKDBL$(AR(J,I))
! LET AI$(8*I+1:8*I+8)=PACKDBL$(AI(J,I))
! NEXT I
! CALL IDFT(N,AR$,AI$)
! FOR I=0 TO N-1
! LET AR(J,I)=UNPACKDBL(AR$(8*I+1:8*I+8))
! LET AI(J,I)=UNPACKDBL(AI$(8*I+1:8*I+8))
! NEXT I
!NEXT J
!END SUB
------------------------------------------------------------------------------
dft.c



#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#define PI 3.14159265358979323846

__declspec(dllexport) void dft(int m, double *xr, double *xi)
{
double *rr,*ii;
double p;
int i,j;
rr = malloc(m * sizeof(double));
ii = malloc(m * sizeof(double));
if (rr == NULL || ii == NULL) exit(1);
#pragma omp parallel for
for (i=0; i<m; i++) {
rr[i]=0;
ii[i]=0;
}
p = 2*PI/(double)m;
#pragma omp parallel for private(i,j)
for (j=0; j<m; j++) {
for (i=0; i<m; i++) {
rr[j] += xr[i] * cos(p * j * i) - xi[i] * sin(p * j * i);
ii[j] += xr[i] * sin(p * j * i) + xi[i] * cos(p * j * i);
}
}
#pragma omp parallel for
for(i = 0; i<m; i++) {
xr[i] = rr[i];
xi[i] = ii[i];
}
free(rr);
free(ii);
}

__declspec(dllexport) void idft(int m, double *xr, double *xi)
{
double *rr,*ii;
double p;
int i,j;
rr = malloc(m * sizeof(double));
ii = malloc(m * sizeof(double));
if (rr == NULL || ii == NULL) exit(1);

p = 2*PI/(double)m;
#pragma omp parallel for
for (i=0; i<m; i++) {
rr[i]=0;
ii[i]=0;
}
#pragma omp parallel for private(i,j)
for (j=0; j<m; j++) {
for (i=0; i<m; i++) {
rr[j] += xr[i] * cos(-p * j * i) - xi[i] * sin(-p * j * i);
ii[j] += xr[i] * sin(-p * j * i) + xi[i] * cos(-p * j * i);
}
}
#pragma omp parallel for
for(i = 0; i<m; i++) {
xr[i] = rr[i]/(double)m;
xi[i] = ii[i]/(double)m;
}
free(rr);
free(ii);
}

__declspec(dllexport) void dft2(int m,int n,int flag,double *rr, double *ii)
{
void idft(int m, double *xr, double *xi);
void dft(int m, double *xr, double *xi);

double *xr,*xi;
int x,y;
xr = malloc(m * sizeof(double));
xi = malloc(m * sizeof(double));
if (rr == NULL || ii == NULL) return;
for(y = 0; y<n; y++) {
#pragma omp parallel for
for(x = 0; x<m; x++) {
xr[x] = rr[y*m+x];
xi[x] = ii[y*m+x];
}
if(flag>0) dft(m, xr, xi);
else idft(m, xr, xi);
#pragma omp parallel for
for(x = 0; x<m; x++) {
rr[y*m+x] = xr[x];
ii[y*m+x] = xi[x];
}
}
free(xr);
free(xi);
xr = malloc(n * sizeof(double));
xi = malloc(n * sizeof(double));

for(x = 0; x<m; x++) {
#pragma omp parallel for
for(y = 0; y<n; y++) {
xr[y] = rr[y*m+x];
xi[y] = ii[y*m+x];
}
if(flag>0) dft (n, xr, xi);
else idft (n, xr, xi);
#pragma omp parallel for
for(y = 0; y<n; y++) {
rr[y*m+x] = xr[y];
ii[y*m+x] = xi[y];
}
}
free(xr);
free(xi);
}

Re: 音声信号処理 - しばっち

2024/10/13 (Sun) 07:57:10

FOR I=1 TO NUM
LET LOUT(I)=F(LIN(I))*.5 ! 非線形変換
LET ROUT(I)=F(RIN(I))*.5
NEXT

EXTERNAL FUNCTION F(X)
LET F=SGN(SIN(ASIN(X))) ! 方形波
!LET F=ASIN(X)/PI*2 ! 三角波
!LET F=SIN(PI*ASIN(X))
!LET F=SIN(2*PI*X)
END FUNCTION
---------------------------------------------------------------
LET R=.4
FOR I=1 TO NUM
LET LOUT(I)=SGN(LIN(I))*ABS(LIN(I))^(1/R) ! 非線形変換
LET ROUT(I)=SGN(RIN(I))*ABS(RIN(I))^(1/R)
NEXT I
---------------------------------------------------------------
RANDOMIZE
LET M=RND*4
LET N=RND*4
SET WINDOW -1,1,-1,1
DRAW GRID(.2,.2)
FOR X=-1 TO 1 STEP 1/32
PLOT LINES:X,INTERPOLATION(X,-1,1,-1,1,M,N);
NEXT X
FOR I=1 TO NUM
LET LOUT(I)=INTERPOLATION(LIN(I),-1,1,-1,1,M,N) ! 非線形変換
LET ROUT(I)=INTERPOLATION(RIN(I),-1,1,-1,1,M,N)
NEXT I

EXTERNAL FUNCTION INTERPOLATION(X,XMIN,XMAX,YMIN,YMAX,M,N) ! XMIN<=X<=XMAX
LET X=MIN(XMAX,MAX(X,XMIN))
LET INTERPOLATION=(YMAX-YMIN)/(SGN(XMAX)*ABS(XMAX)^M-XMIN)^N*(SGN(X)*ABS(X)^M-XMIN)^N+YMIN
END FUNCTION
---------------------------------------------------------------
RANDOMIZE
LET M=RND*5
LET N=RND*5
LET P=RND*5
SET WINDOW -1,1,-1,1
DRAW GRID(.2,.2)
FOR X=-1 TO 1 STEP 1/32
PLOT LINES:X,INTERPOLATION(X,-1,1,-1,1,M,N,P);
NEXT X
FOR I=1 TO NUM
LET LOUT(I)=INTERPOLATION(LIN(I),-1,1,-1,1,M,N,P) ! 非線形変換
LET ROUT(I)=INTERPOLATION(RIN(I),-1,1,-1,1,M,N,P)
NEXT I

EXTERNAL FUNCTION INTERPOLATION(X,XMIN,XMAX,YMIN,YMAX,M,N,P) ! XMIN<=X<=XMAX
LET X=MIN(XMAX,MAX(X,XMIN))
LET T=(X-XMIN)/(XMAX-XMIN)
LET INTERPOLATION=YMIN*(1-T^M)^N+YMAX*T^P
END FUNCTION
---------------------------------------------------------------
LET N=10
DIM X(N),Y(N)
SET WINDOW -1,1,-1,1
DRAW GRID(.2,.2)
FOR I=1 TO N
READ X(I),Y(I)
PLOT LINES:X(I),Y(I);
NEXT I
DATA -1,1
DATA -.6,-1
DATA -.6,1
DATA -.2,-1
DATA -.2,1
DATA .2,-1
DATA .2,1
DATA .6,-1
DATA .6,1
DATA 1,-1
FOR I=1 TO NUM
LET LOUT(I)=F(LIN(I),N,X,Y) ! 線形変換
LET ROUT(I)=F(RIN(I),N,X,Y)
NEXT I

EXTERNAL FUNCTION F(XX,N,X(),Y()) !
FOR I=1 TO N-1
IF X(I)<=XX AND X(I+1)>=XX THEN
LET F=(Y(I+1)-Y(I))/(X(I+1)-X(I))*(XX-X(I))+Y(I)
EXIT FUNCTION
END IF
NEXT I
LET F=XX
END FUNCTION
---------------------------------------------------------------
DIM LUT(255)
SET WINDOW -1,1,-1,1
DRAW GRID(.2,.2)
FOR I=1 TO 255
READ LUT(I) ! ルックアップテーブル
PLOT LINES:2*I/255-1,LUT(I);
NEXT I
DATA -7.85426707962836E-3 ,-1.57141319943974E-2 ,-.023461508734477 ,-3.09781144044374E-2 ,-.038147213502212 ,-4.48553661866093E-2 ,-5.09941552129774E-2 ,-5.64618661668915E-2 ,-6.11650959381815E-2 ,-6.50202650522969E-2 ,-6.79550105078911E-2 ,-6.99094371484411E-2 ,-7.08372073065202E-2 ,-7.07064504819803E-2 ,-6.95004771250992E-2 ,-6.72182831637032E-2
DATA -6.38748347064068E-2 ,-5.95011253358703E-2 ,-5.41440015367571E-2 ,-4.78657550406853E-2 ,-4.07434841707331E-2 ,-3.28682295854111E-2 ,-.024343893110078 ,-1.52859515560646E-2 ,-5.81998051829323E-3 , 3.91999393392034E-3 , 1.37932950439433E-2 , 2.36545235399083E-2 , 3.33555889229312E-2 , 4.27478167486482E-2 , 5.16841032773653E-2 , 6.00210877438072E-2
DATA 6.76213117712769E-2 , 7.43553351423367E-2 , 8.01037772486475E-2 , 8.47592540804446E-2 , 8.82281815794658E-2 , 9.04324175603358E-2 , 9.13107161909428E-2 , .090819971192986 , 8.89362264549138E-2 , 8.56554356109361E-2 , 8.09939552968595E-2 , 7.49887602068642E-2 , .06769737170178 , .059197495512273 , 4.95863679901577E-2 , 3.89798143361773E-2
DATA 2.75110262199412E-2 , 1.53290701533799E-2 , 2.59714182913288E-3 ,-1.05094146635328E-2 ,-2.38053003322234E-2 ,-3.70976418213195E-2 ,-5.01885303826588E-2 ,-6.28776880496255E-2 ,-7.49652275556813E-2 ,-8.62544707216416E-2 ,-9.65547886636063E-2 ,-.105684426269286 ,-.113473272972373 ,-.119765541934187 ,-.124422320324079 ,-.127323954473516
DATA -.128372235255002 ,-.127492351091231 ,-.124634578510639 ,-.11977568310519 ,-.112920007081167 ,-.10410022328474 ,-9.33777395869012E-2 ,-8.08427417780899E-2 ,-6.66138675988017E-2 ,-5.08375091625488E-2 ,-3.36867457531462E-2 ,-1.53599137390094E-2 , 3.9211749184947E-3 , .023913349439273 , 4.43545116607216E-2 , 6.49663572269623E-2
DATA 8.54573576000593E-2 , .105525970261009 , .124864041138839 , .143160360394923 , .160104330237978 , .17538970148906 , .188718334184817 , .199803936624569 , .208375736948562 , .214182041590269 , .216993635776524 , .216606982650025 , .212847179546406 , .205570632452343 , .194667412674888 , .180063263231421
DATA .161721226384014 , .139642868045584 , .113869079426558 , 8.44804412139577E-2 , 5.15971407199358E-2 , 1.53784377407642E-2 ,-2.39783197353682E-2 ,-6.62391234096138E-2 ,-.111135164198182 ,-.158365106184254 ,-.207597752679449 ,-.258475077219701 ,-.310615587763502 ,-.363617988163652 ,-.417065097200656 ,-.470527982145922
DATA -.52357026100908 ,-.575752525353269 ,-.626636833865371 ,-.675791225768264 ,-.722794202674842 ,-.767239127617103 ,-.808738490738177 ,-.846927992503372 ,-.881470397252895 ,-.912059112461055 ,-.938421452154413 ,-.960321546537559 ,-.97756286393654 ,-.989990315647039 ,-.997491919112793 , 1
DATA .997491919112793 , .989990315647039 , .97756286393654 , .960321546537559 , .938421452154413 , .912059112461055 , .881470397252895 , .846927992503372 , .808738490738177 , .767239127617103 , .722794202674842 , .675791225768264 , .626636833865371 , .575752525353269 , .52357026100908 , .470527982145922
DATA .417065097200656 , .363617988163652 , .310615587763502 , .258475077219701 , .207597752679449 , .158365106184254 , .111135164198182 , 6.62391234096138E-2 , 2.39783197353682E-2 ,-1.53784377407642E-2 ,-5.15971407199358E-2 ,-8.44804412139577E-2 ,-.113869079426558 ,-.139642868045584 ,-.161721226384014 ,-.180063263231421
DATA -.194667412674888 ,-.205570632452343 ,-.212847179546406 ,-.216606982650025 ,-.216993635776524 ,-.214182041590269 ,-.208375736948562 ,-.199803936624569 ,-.188718334184817 ,-.17538970148906 ,-.160104330237978 ,-.143160360394923 ,-.124864041138839 ,-.105525970261009 ,-8.54573576000593E-2 ,-6.49663572269623E-2
DATA -4.43545116607216E-2 ,-.023913349439273 ,-3.9211749184947E-3 , 1.53599137390094E-2 , 3.36867457531462E-2 , 5.08375091625488E-2 , 6.66138675988017E-2 , 8.08427417780899E-2 , 9.33777395869012E-2 , .10410022328474 , .112920007081167 , .11977568310519 , .124634578510639 , .127492351091231 , .128372235255002 , .127323954473516
DATA .124422320324079 , .119765541934187 , .113473272972373 , .105684426269286 , 9.65547886636063E-2 , 8.62544707216416E-2 , 7.49652275556813E-2 , 6.28776880496255E-2 , 5.01885303826588E-2 , 3.70976418213195E-2 , 2.38053003322234E-2 , 1.05094146635328E-2 ,-2.59714182913288E-3 ,-1.53290701533799E-2 ,-2.75110262199412E-2 ,-3.89798143361773E-2
DATA -4.95863679901577E-2 ,-.059197495512273 ,-.06769737170178 ,-7.49887602068642E-2 ,-8.09939552968595E-2 ,-8.56554356109361E-2 ,-8.89362264549138E-2 ,-.090819971192986 ,-9.13107161909428E-2 ,-9.04324175603358E-2 ,-8.82281815794658E-2 ,-8.47592540804446E-2 ,-8.01037772486475E-2 ,-7.43553351423367E-2 ,-6.76213117712769E-2 ,-6.00210877438072E-2
DATA -5.16841032773653E-2 ,-4.27478167486482E-2 ,-3.33555889229312E-2 ,-2.36545235399083E-2 ,-1.37932950439433E-2 ,-3.91999393392034E-3 , 5.81998051829323E-3 , 1.52859515560646E-2 , .024343893110078 , 3.28682295854111E-2 , 4.07434841707331E-2 , 4.78657550406853E-2 , 5.41440015367571E-2 , 5.95011253358703E-2 , 6.38748347064068E-2 , 6.72182831637032E-2
DATA 6.95004771250992E-2 , 7.07064504819803E-2 , 7.08372073065202E-2 , 6.99094371484411E-2 , 6.79550105078911E-2 , 6.50202650522969E-2 , 6.11650959381815E-2 , 5.64618661668915E-2 , 5.09941552129774E-2 , 4.48553661866093E-2 , .038147213502212 , 3.09781144044374E-2 , .023461508734477 , 1.57141319943974E-2 , 7.85426707962836E-3
FOR I=1 TO NUM
LET LOUT(I)=LUT(LIN(I)*127+128)
LET ROUT(I)=LUT(RIN(I)*127+128)
NEXT I
---------------------------------------------------------------
DATA -.997198879551821 ,-.994397759103641 ,-.991596638655462 ,-.988795518207283 ,-.985994397759104 ,-.983193277310924 ,-.980392156862745 ,-.977591036414566 ,-.974789915966387 ,-.971988795518207 ,-.969187675070028 ,-.966386554621849 ,-.963585434173669 ,-.96078431372549 ,-.957983193277311 ,-.955182072829132
DATA -.952380952380953 ,-.949579831932773 ,-.946778711484594 ,-.943977591036415 ,-.941176470588235 ,-.938375350140056 ,-.935574229691877 ,-.932773109243698 ,-.929971988795518 ,-.927170868347339 ,-.92436974789916 ,-.92156862745098 ,-.918767507002801 ,-.915966386554622 ,-.913165266106443 ,-.910364145658263
DATA -.907563025210084 ,-.904761904761905 ,-.901960784313726 ,-.899159663865546 ,-.896358543417367 ,-.893557422969188 ,-.890756302521008 ,-.887955182072829 ,-.88515406162465 ,-.882352941176471 ,-.879551820728291 ,-.876750700280112 ,-.873949579831933 ,-.871148459383753 ,-.868347338935574 ,-.865546218487395
DATA -.862745098039216 ,-.859943977591036 ,-.857142857142857 ,-.854341736694678 ,-.851540616246499 ,-.848739495798319 ,-.84593837535014 ,-.843137254901961 ,-.840336134453782 ,-.837535014005602 ,-.834733893557423 ,-.831932773109244 ,-.829131652661064 ,-.826330532212885 ,-.823529411764706 ,-.820728291316527
DATA -.817927170868347 ,-.815126050420168 ,-.812324929971989 ,-.80952380952381 ,-.80672268907563 ,-.803921568627451 ,-.801120448179272 ,-.798319327731092 ,-.795518207282913 ,-.792717086834734 ,-.789915966386555 ,-.787114845938375 ,-.784313725490196 ,-.781512605042017 ,-.778711484593838 ,-.775910364145658
DATA -.773109243697479 ,-.7703081232493 ,-.76750700280112 ,-.764705882352941 ,-.761904761904762 ,-.759103641456583 ,-.756302521008403 ,-.753501400560224 ,-.750700280112045 ,-.747899159663865 ,-.745098039215686 ,-.742296918767507 ,-.739495798319328 ,-.736694677871148 ,-.733893557422969 ,-.73109243697479
DATA -.728291316526611 ,-.725490196078431 ,-.722689075630252 ,-.719887955182073 ,-.717086834733894 ,-.714285714285714 ,-.686274509803921 ,-.658263305322129 ,-.630252100840336 ,-.602240896358543 ,-.574229691876751 ,-.546218487394958 ,-.518207282913165 ,-.490196078431372 ,-.46218487394958 ,-.434173669467787
DATA -.406162464985995 ,-.378151260504202 ,-.350140056022409 ,-.322128851540616 ,-.294117647058824 ,-.266106442577031 ,-.238095238095238 ,-.210084033613445 ,-.182072829131653 ,-.15406162464986 ,-.126050420168067 ,-9.80392156862746E-2 ,-7.00280112044819E-2 ,-4.20168067226889E-2 ,-1.40056022408963E-2 , 1.40056022408963E-2
DATA 4.20168067226889E-2 , 7.00280112044815E-2 , 9.80392156862741E-2 , .126050420168068 , .15406162464986 , .182072829131653 , .210084033613445 , .238095238095238 , .266106442577031 , .294117647058823 , .322128851540616 , .350140056022409 , .378151260504202 , .406162464985995 , .434173669467787 , .46218487394958
DATA .490196078431372 , .518207282913165 , .546218487394958 , .574229691876751 , .602240896358544 , .630252100840336 , .658263305322129 , .686274509803921 , .714285714285714 , .717086834733893 , .719887955182073 , .722689075630252 , .725490196078431 , .728291316526611 , .73109243697479 , .733893557422969
DATA .736694677871148 , .739495798319328 , .742296918767507 , .745098039215686 , .747899159663866 , .750700280112045 , .753501400560224 , .756302521008403 , .759103641456583 , .761904761904762 , .764705882352941 , .76750700280112 , .7703081232493 , .773109243697479 , .775910364145658 , .778711484593838
DATA .781512605042017 , .784313725490196 , .787114845938375 , .789915966386555 , .792717086834734 , .795518207282913 , .798319327731092 , .801120448179272 , .803921568627451 , .80672268907563 , .80952380952381 , .812324929971989 , .815126050420168 , .817927170868347 , .820728291316527 , .823529411764706
DATA .826330532212885 , .829131652661064 , .831932773109244 , .834733893557423 , .837535014005602 , .840336134453782 , .843137254901961 , .84593837535014 , .848739495798319 , .851540616246499 , .854341736694678 , .857142857142857 , .859943977591036 , .862745098039216 , .865546218487395 , .868347338935574
DATA .871148459383753 , .873949579831933 , .876750700280112 , .879551820728291 , .882352941176471 , .88515406162465 , .887955182072829 , .890756302521008 , .893557422969188 , .896358543417367 , .899159663865546 , .901960784313726 , .904761904761905 , .907563025210084 , .910364145658263 , .913165266106443
DATA .915966386554622 , .918767507002801 , .92156862745098 , .92436974789916 , .927170868347339 , .929971988795518 , .932773109243698 , .935574229691877 , .938375350140056 , .941176470588235 , .943977591036415 , .946778711484594 , .949579831932773 , .952380952380953 , .955182072829132 , .957983193277311
DATA .96078431372549 , .963585434173669 , .966386554621849 , .969187675070028 , .971988795518207 , .974789915966387 , .977591036414566 , .980392156862745 , .983193277310924 , .985994397759104 , .988795518207283 , .991596638655462 , .994397759103641 , .997198879551821 , 1
---------------------------------------------------------------
LET THRESHOLD=.2
FOR K=1 TO NUM
IF ABS(LIN(K))<THRESHOLD THEN ! コンプレッサー
LET LOUT(K)=SGN(LIN(K))*INTERPOLATE(ABS(LIN(K)),0,THRESHOLD,0,.8)
ELSE
LET LOUT(K)=SGN(LIN(K))*INTERPOLATE(ABS(LIN(K)),THRESHOLD,1,.8,1)
END IF
IF ABS(RIN(K))<THRESHOLD THEN
LET ROUT(K)=SGN(RIN(K))*INTERPOLATE(ABS(RIN(K)),0,THRESHOLD,0,.8)
ELSE
LET ROUT(K)=SGN(RIN(K))*INTERPOLATE(ABS(RIN(K)),THRESHOLD,1,.8,1)
END IF
NEXT K

EXTERNAL FUNCTION INTERPOLATE(X,XMIN,XMAX,YMIN,YMAX) ! XMIN<=X<=XMAX
LET X=MIN(XMAX,MAX(X,XMIN))
LET INTERPOLATE=(X-XMIN)*(YMAX-YMIN)/(XMAX-XMIN) + YMIN
END FUNCTION
---------------------------------------------------------------
LET N=32 ! 次数(偶数)
DIM COEF(0 TO N)
LET FILTER=2
LET FREQ=.5
FOR K=1 TO NUM ! オートワウ、ワウワウ
LET FE=2000+800*SIN(FREQ/SAMPLINGFREQ*K*2*PI) ! LFO(低周波発振器)
SELECT CASE FILTER
CASE 0
CALL GETCOEFFICIENTLPF(FE,SAMPLINGFREQ,N,COEF)
CASE 1
CALL GETCOEFFICIENTHPF(FE,SAMPLINGFREQ,N,COEF)
CASE 2
CALL GETCOEFFICIENTBPF(FE,FE+500,SAMPLINGFREQ,N,COEF)
END SELECT
LET L=0
LET R=0
FOR J=0 TO N
IF K+J<=NUM THEN
LET L=L+LIN(K+J)*COEF(J)
LET R=R+RIN(K+J)*COEF(J)
END IF
NEXT J
LET LOUT(K)=L
LET ROUT(K)=R
NEXT K
---------------------------------------------------------------
画像処理でいうガラス加工


LET N=64
FOR K=1 TO NUM STEP N
FOR I=0 TO N-1
LET LOUT(K+I)=LIN(K+N-I) ! 入れ替え
LET ROUT(K+I)=RIN(K+N-I)
NEXT I
NEXT K
---------------------------------------------------------------
画像処理でいう点画


RANDOMIZE
LET N=32
FOR K=1 TO NUM
DO
LET P=INT((RND-.5)*N)
LOOP UNTIL K+P<=NUM AND K+P>=1
LET LOUT(K)=LIN(K+P) ! 乱数入れ替え
LET ROUT(K)=RIN(K+P)
NEXT K
---------------------------------------------------------------
LET SIZE=44100*180
DIM LIN1(SIZE),RIN1(SIZE),LOUT(2*SIZE),ROUT(2*SIZE)
DIM LIN2(SIZE),RIN2(SIZE)
FILE GETOPENNAME F$,"読み込み1　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM1,CHANNEL,SAMPLINGFREQ1,LIN1,RIN1,ERR)
IF ERR>0 THEN STOP
FILE GETOPENNAME F$,"読み込み2　WAVファイル|*.WAV"
IF F$="" THEN STOP
CALL READWAV(F$,NUM2,CHANNEL,SAMPLINGFREQ2,LIN2,RIN2,ERR)
IF ERR>0 THEN STOP
IF SAMPLINGFREQ1<>SAMPLINGFREQ2 THEN
PRINT "サンプリング周波数が違います"
STOP
END IF
FOR K=1 TO NUM1 ! ファイル連結
LET N=N+1
LET LOUT(N)=LIN1(K)
LET ROUT(N)=RIN1(K)
NEXT K
FOR K=1 TO NUM2
LET N=N+1
LET LOUT(N)=LIN2(K)
LET ROUT(N)=RIN2(K)
NEXT K
---------------------------------------------------------------
https://ja.wikipedia.org/wiki/平均
https://ja.wikipedia.org/wiki/ヘルダー平均


IF SAMPLINGFREQ1<>SAMPLINGFREQ2 THEN
PRINT "サンプリング周波数が違います"
STOP
END IF
LET NN=MAX(NUM1,NUM2)
LET MODE=0
FOR K=1 TO NN
IF K<=NUM1 AND K<=NUM2 THEN ! ファイルミックス(ファイル合成)
SELECT CASE MODE
CASE 0
LET LOUT(K)=(LIN1(K)+LIN2(K))/2 ! 相加平均
LET ROUT(K)=(RIN1(K)+RIN2(K))/2
CASE 1
LET LOUT(K)=MAX(LIN1(K),LIN2(K))
LET ROUT(K)=MAX(RIN1(K),RIN2(K))
CASE 2
LET LOUT(K)=MIN(LIN1(K),LIN2(K))
LET ROUT(K)=MIN(RIN1(K),RIN2(K))
CASE 3
LET LOUT(K)=(LIN1(K)-LIN2(K))
LET ROUT(K)=(RIN1(K)-RIN2(K))
CASE 4
LET LOUT(K)=SGN_(LIN1(K))*SGN_(LIN2(K))*SQR(ABS(LIN1(K)*LIN2(K))) ! 相乗平均
LET ROUT(K)=SGN_(RIN1(K))*SGN_(RIN2(K))*SQR(ABS(RIN1(K)*RIN2(K)))
CASE 5
IF LIN1(K)=0 OR LIN2(K)=0 THEN
LET LOUT(K)=0
ELSE
LET LOUT(K)=SGN_(LIN1(K))*SGN_(LIN2(K))/(ABS(1/LIN1(K))+ABS(1/LIN2(K))) ! 調和平均 1/(1/M+1/N)=N*M/(N+M)
END IF
IF RIN1(K)=0 OR RIN2(K)=0 THEN
LET ROUT(K)=0
ELSE
LET ROUT(K)=SGN_(RIN1(K))*SGN_(RIN2(K))/(ABS(1/RIN1(K))+ABS(1/RIN2(K)))
END IF
CASE 6
LET LOUT(K)=SGN_(LIN1(K))*SGN_(LIN2(K))*SQR(LIN1(K)^2+LIN2(K)^2)/2 ! 2乗平均
LET ROUT(K)=SGN_(RIN1(K))*SGN_(RIN2(K))*SQR(RIN1(K)^2+RIN2(K)^2)/2
CASE 7
LET M=2.5
LET N=3.5
LET LOUT(K)=SGN_(LIN1(K))*SGN_(LIN2(K))*(ABS(LIN1(K))^M+ABS(LIN2(K))^N)^(2/(M+N)) ! ミンコフスキー距離(ヘルダー平均)
LET ROUT(K)=SGN_(RIN1(K))*SGN_(RIN2(K))*(ABS(RIN1(K))^M+ABS(RIN2(K))^N)^(2/(M+N))

! LET LOUT(K)=SGN_(LIN1(K))*SGN_(LIN2(K))*(ABS(LIN1(K))^(1/M)+ABS(LIN2(K))^(1/N))^(2*N*M/(M+N))
! LET ROUT(K)=SGN_(RIN1(K))*SGN_(RIN2(K))*(ABS(RIN1(K))^(1/M)+ABS(RIN2(K))^(1/N))^(2*N*M/(M+N))
CASE 8
LET LOUT(K)=PXOR(ABS(LIN1(K)),ABS(LIN2(K)))*SGN_(LIN1(K))*SGN_(LIN2(K)) ! XOR
LET ROUT(K)=PXOR(ABS(RIN1(K)),ABS(RIN2(K)))*SGN_(RIN1(K))*SGN_(RIN2(K))
CASE 9
LET LOUT(K)=PEQV(ABS(LIN1(K)),ABS(LIN2(K)))*SGN_(LIN1(K))*SGN_(LIN2(K)) ! EQV
LET ROUT(K)=PEQV(ABS(RIN1(K)),ABS(RIN2(K)))*SGN_(RIN1(K))*SGN_(RIN2(K))
CASE 10
LET LOUT(K)=PIMP(ABS(LIN1(K)),ABS(LIN2(K)))*SGN_(LIN1(K))*SGN_(LIN2(K)) ! IMP
LET ROUT(K)=PIMP(ABS(RIN1(K)),ABS(RIN2(K)))*SGN_(RIN1(K))*SGN_(RIN2(K))
CASE 11
LET LOUT(K)=PNOR(ABS(LIN1(K)),ABS(LIN2(K)))*SGN_(LIN1(K))*SGN_(LIN2(K)) ! NOR
LET ROUT(K)=PNOR(ABS(RIN1(K)),ABS(RIN2(K)))*SGN_(RIN1(K))*SGN_(RIN2(K))
CASE 12
LET LOUT(K)=PNAND(ABS(LIN1(K)),ABS(LIN2(K)))*SGN_(LIN1(K))*SGN_(LIN2(K)) ! NAND
LET ROUT(K)=PNAND(ABS(RIN1(K)),ABS(RIN2(K)))*SGN_(RIN1(K))*SGN_(RIN2(K))
CASE 13
LET LOUT(K)=PNIMP(ABS(LIN1(K)),ABS(LIN2(K)))*SGN_(LIN1(K))*SGN_(LIN2(K)) ! NIMP
LET ROUT(K)=PNIMP(ABS(RIN1(K)),ABS(RIN2(K)))*SGN_(RIN1(K))*SGN_(RIN2(K))
END SELECT
ELSEIF K<=NUM2 AND K>NUM1 THEN
LET LOUT(K)=LIN2(K)
LET ROUT(K)=RIN2(K)
ELSEIF K<=NUM1 AND K>NUM2 THEN
LET LOUT(K)=LIN1(K)
LET ROUT(K)=RIN1(K)
END IF
NEXT K
---------------------------------------------------------------
IF CHANNEL<>2 THEN
PRINT "ステレオではありません"
STOP
END IF
CALL WRITEWAV("左チャンネル.wav",NUM,SAMPLINGFREQ,16,LOUT) ! ステレオ分解　モノラル出力
CALL WRITEWAV("右チャンネル.wav",NUM,SAMPLINGFREQ,16,ROUT)
END
---------------------------------------------------------------
INPUT PROMPT "分割点(秒)=":PLAYTIME
LET NN=SAMPLINGFREQ*PLAYTIME
FOR K=1 TO NN
LET LOUT(K)=LIN(K)
LET ROUT(K)=RIN(K)
NEXT K
CALL WRITEWAV("分割1.wav",NN,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LOUT,ROUT) ! ファイル分割
MAT LOUT=ZER
MAT ROUT=ZER
FOR K=NN+1 TO NUM
LET N=N+1
LET LOUT(N)=LIN(K)
LET ROUT(N)=RIN(K)
NEXT K
CALL WRITEWAV("分割2.wav",N,CHANNEL,SAMPLINGFREQ,SAMPLEBIT,LOUT,ROUT)
END

Re: コンソールアプリ - しばっち

2024/10/13 (Sun) 07:40:51

C/C++ライブラリーによるコンソールアプリでwavファイルからcsvファイルを書き出します。

https://www.un4seen.com/

これはDLL化には失敗した(初期化に失敗。原因不明)ものの、wavファイル以外にもmp3を直接読み込めるので作成しました。
その他mod形式と呼ばれるitやs3mファイル等にも対応しているようですがよくわかりません。

実行には別途BASS.DLLが必要です。

csvファイルなのでエクセル等でも読み込めます。
wav2csv.exeはwin32版です。win2csv64.exeはwin64版になります。

コマンドプロンプトから

wav2csv test.wav 標準出力へ書き出し。
wav2csv input.wav sample.csv wavファイルからcsvファイルへ書き出し。
wav2csv input.mp3 sample.csv 3 5 mp3ファイルの3秒から5秒までを切り出してcsvファイル化

-1～1の正規値を出力します。
マルチチャンネル出力可能ですが拡張フォーマットには対応していないようです。

VC++2022でコンパイルしました。
下記からダウンロードしてください。(wave.zip)

https://66.gigafile.nu/0121-d5d86eb73a1f6ce256c440f69c9462329
ダウンロード期限：2025年1月21日(火)
パスワード:設定していません



----------------------------------------------------------------------------------------
wav2csv.cpp


#include <cstdlib>
#include <cstdio>
#include "bass.h"
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <iomanip>

#pragma comment(lib,"bass.lib")
using namespace std;

int main(int argc, char *argv[])
{
BASS_CHANNELINFO info;
DWORD chan,bpf,p;
QWORD pos;
double secs;
int n=0,no,samplebit,i,j,flg=0;
double starttime=-1.0,endtime=-1.0;
char *wavname,*csvname;

if (argc==1) return 1;
if (argc>=2) wavname=argv[1];
if (argc>=3) {
csvname=argv[2];
flg=1;
}
if (argc>=4) starttime=atof(argv[3]);
if (argc==5) endtime=atof(argv[4]);

if (!BASS_Init(-1, 44100, 0, 0, NULL)) return 1;

chan = BASS_StreamCreateFile(FALSE, wavname, 0, 0, BASS_STREAM_DECODE);
if (!chan && BASS_ErrorGetCode() == BASS_ERROR_FILEFORM) {
chan=BASS_MusicLoad(FALSE,wavname,0,0,BASS_MUSIC_DECODE | BASS_MUSIC_RAMPS | BASS_MUSIC_PRESCAN,0);
}
if(!chan) {
BASS_Free();
return 1;
}
BASS_ChannelGetInfo(chan,&info);
samplebit=(info.flags & BASS_SAMPLE_8BITS ? 8 : 16);
bpf=info.chans*(info.flags & BASS_SAMPLE_8BITS ? 1 : 2);
pos=BASS_ChannelGetLength(chan,BASS_POS_BYTE);
if(pos==0) {
BASS_Free();
return 1;
} else {
secs=BASS_ChannelBytes2Seconds(chan,pos);
}
if(starttime>=0.0) if (!BASS_ChannelSetPosition(chan, starttime * bpf, BASS_POS_BYTE))
{
BASS_Free();
return 1;
}
if(starttime<endtime && endtime>0.0) if (!BASS_ChannelSetPosition(chan, endtime * bpf, BASS_POS_END))
{
BASS_Free();
return 1;
}
if (flg==0) {
cout << setprecision(15) << info.freq << endl; // サンプリング周波数
cout << setprecision(15) << samplebit << endl; // サンプルビット
cout << setprecision(15) << info.chans << endl; // チャンネル数
cout << setprecision(15) << pos/bpf << endl; // データ数
} else {
ofstream ofs(csvname);
ofs.precision(15);
ofs << info.freq << endl; // サンプリング周波数
ofs << samplebit << endl; // サンプルビット
ofs << info.chans << endl; // チャンネル数
ofs << pos/bpf << endl; // データ数
}
while (1) {
short buf[10000];
ofstream ofs(csvname, std::ios::out | std::ios::app);
ofs.precision(15);

int c=BASS_ChannelGetData(chan,buf,sizeof(buf));
if (c==-1) break;
for (i=0; i<c; i+=info.chans)
for(j=0; j<info.chans; j++) {
if (flg==0) cout << setprecision(15) << (double)buf[i+j]/pow(2.0,samplebit-1);
else {
ofs << (double)buf[i+j]/pow(2.0,samplebit-1);
}
if (j<info.chans-1) {
if (flg==0) cout << ",";
else ofs << ",";
}
else {
if (flg==0) cout << endl;
else ofs << endl;
}
}
}
BASS_Free();
return 0;
}
----------------------------------------------------------------------------------------
使用例

下記のようにするとGUI化できます。


FILE GETOPENNAME F$,"読み込み　対応ファイル|*.wav;*.aiff;*.mp3;*.mp2;*.mp1;*.ogg;*.mo3;*.it;*.s3m;*.mtm;*.mod;*.umx;*.xm"
IF F$="" THEN STOP
LET TEMP$="temporary.csv"
LET STARTTIME=-1 ! 負数では無効。STARTTIME>=0で有効。秒数指定
LET ENDTIME=-1 ! 負数で無効。ENDTIME>0 STARTTIME<ENDTIME で有効
EXECUTE "wav2csv.exe" WITH (F$,TEMP$) ! パスを指定。
!!EXECUTE "wav2csv.exe" WITH (F$,TEMP$,STARTTIME,ENDTIME)
END

Re: コンソールアプリ - しばっち

2024/10/13 (Sun) 07:41:57

C/C++ライブラリーによるコンソールアプリはwavファイルからcsvファイルを書き出します。

https://github.com/adamstark/AudioFile

ヘッダーライブラリーです。ヘッダーファイルをインクルードするだけで使えます。
上記URLからZIPファイルをダウンロードしてVC++2022で下記cppソースからコンパイルできます。

コマンドプロンプトから

wave2csv sample.wav sample.csv

のように打ち込んでください。


wave2csv.exeはwin32版です。wave2csv64.exeはwin64版になります。

-1～1の正規値を出力します。


VC++2022でコンパイルしました。

マルチチャンネル出力可能ですが拡張フォーマットには対応していないようです。


----------------------------------------------------------------------------------------
wave2csv.cpp



#include "AudioFile.h"
#include <iostream>
#include <fstream>

using namespace std;


int main(int argc,char **argv)
{
AudioFile<double> audioFile;
int samplingfreq,channel,samplebit,num,i,ch;
char *wavname,*csvname;

if (argc!=3) return 1;
wavname=argv[1];
csvname=argv[2];

ofstream csv(csvname,ios::out);

if(!audioFile.load(wavname)) return 1;

samplingfreq = audioFile.getSampleRate();
samplebit = audioFile.getBitDepth();
channel = audioFile.getNumChannels();
num = audioFile.getNumSamplesPerChannel();

ofstream ofs(csvname);
if(!ofs) return 1;
ofs << samplingfreq << endl; // サンプリング周波数
ofs << samplebit << endl; // サンプルビット
ofs << channel << endl; // チャンネル数
ofs << num << endl; // データ数

for (i = 0; i < num; i++) {
for (ch=0; ch<channel; ch++) {
ofs << audioFile.samples[ch][i];
if (ch<channel-1) ofs << ",";
else ofs << endl;
}
}
return 0;
}
----------------------------------------------------------------------------------------
使用例

次のようにするとGUI化できます。



FILE GETOPENNAME F$,"読み込み　WAVファイル|*.WAV"
IF F$="" THEN STOP
FILE GETSAVENAME G$,"保存 CSVファイル|*.CSV"
IF G$="" THEN STOP
IF POS(LCASE$(G$),".csv")=0 THEN LET G$=G$&".csv"
EXECUTE "wave2csv.exe" WITH (F$,G$)
END