shimauma3sanのブログ

畳み込みの計算をFFTを使って効率よく（計算回数を減らす）計算する方法を説明します。　畳み込み(離散畳み込み）、フーリエ変換（離散フーリエ変換:DFT）、フーリエ逆変換（離散フーリエ逆変換：IDFT)は、知っているものとします。　（以降、フーリエ変換、フーリエ逆変換という言葉を使います。）　
　説明は、数式処理ソフト maxima での計算例で示します。( maxima はフリーソフトです。)

畳み込み、フーリエ変換の関係式の整理

$x_n, h_n, y_n$ のフーリエ変換を、それぞれ $X_k, H_k, Y_k$ とすると、フーリエ変換（左側）、フーリエ逆変換（右側）は次のようになる。

${ \displaystyle \begin{eqnarray} X_k = \sum_{n=0}^{N-1}x_ne^{-j\frac{2\pi kn}N}　~~~ x_n = \frac{1}N\sum_{k=0}^{N-1}X_ke^{j\frac{2\pi kn}N} \tag{1} \\ H_k = \sum_{n=0}^{N-1}x_ne^{-j\frac{2\pi kn}N}　~~~ h_n = \frac{1}N\sum_{k=0}^{N-1}X_ke^{j\frac{2\pi kn}N} \tag{2} \\ Y_k = \sum_{n=0}^{N-1}x_ne^{-j\frac{2\pi kn}N}　~~~ y_n = \frac{1}N\sum_{k=0}^{N-1}X_ke^{j\frac{2\pi kn}N} \tag{3} \end{eqnarray} }$

$x_n$ が、フィルタ $h_n$ に入力された時の出力を $y_n$ とすると、 $y_n$ は、 $x_n$ と $h_n$ の畳み込みで表される。
${ \displaystyle \begin{eqnarray} y_n &= x_n * h_n \tag{4}\\ \end{eqnarray} }$
リアルタイムに(4)式の畳み込みを、効率よく（計算量を減らす）行える、フーリエ変換・逆変換を使う手法、重量保留法と重畳加算法を説明する。以下関係式を整理する。

(4) 式の両辺をフーリエ変換する。

${ \displaystyle \begin{eqnarray} \mathcal{F} (y_n) &= \mathcal{F}(x_n * h_n) \tag{5} \\ \end{eqnarray} }$
畳み込みのフーリエ変換は、それぞれのフーリエ変換の積。

${ \displaystyle \begin{eqnarray} \mathcal{F} (y_n) &= X_k H_k \tag{6}\\ \end{eqnarray} }$
(6)式の両辺をそれぞれ逆フーリエ変換する。

${ \displaystyle \begin{eqnarray} y_n &= \mathcal{F}^{-1}(X_k H_k) \tag{7} \\ \end{eqnarray} }$

畳み込みをフーリエ変換、フーリエ逆変換を使って計算する

　次のような4点の入力系列 $x_n$ と、4点のフィルタ $h_n$ の畳み込みを行う。

　 $x_n = (x_0, x_1, x_2, x_3, 0, 0, 0, 0)$
　 $h_n = (h_0, h_1, h_2, h_3, 0, 0 ,0, 0)$

数式処理ソフト maxima で使う。maxima コードの説明をしよう。フーリエ変換フーリエ逆変換を使って、入力 $x_n$ とフィルタ $h_n$ の畳み込みを行う。 $x_n$ と $h_n$ は、それぞれ4点のデータで、残りはゼロ埋めをしている。入力4点とフィルタ4点のデータを畳み込むと7点のデータとなる。maxima では、それぞれ $x_n, h_n$ を定義して、それぞれをフーリエ変換する。それぞれのフーリエ変換の結果を掛算し( $H_k X_k$ )、その結果を逆フーリエ変換(7式)すると畳み込みが求まる。（maxima コードのコメント参照）

maxima コード

N : 8;
c : (2.0 * %pi / N)$;
ef(kn) := cos(c*kn) - %i * sin(c*kn)$;
F : matrix([ef(0), ef(0), ef(0), ef(0),  ef(0), ef(0), ef(0), ef(0)],
           [ef(0), ef(1), ef(2), ef(3),  ef(4), ef(5), ef(6), ef(7)],
           [ef(0), ef(2), ef(4), ef(6),  ef(8), ef(10), ef(12), ef(14)],
           [ef(0), ef(3), ef(6), ef(9),  ef(12), ef(15), ef(18), ef(21)],
           [ef(0), ef(4), ef(8), ef(12),  ef(16), ef(20), ef(24), ef(28)],
           [ef(0), ef(5), ef(10), ef(15),  ef(20), ef(25), ef(30), ef(35)],
           [ef(0), ef(6), ef(12), ef(18),  ef(24), ef(30), ef(36), ef(42)],
           [ef(0), ef(7), ef(14), ef(21),  ef(28), ef(35), ef(42), ef(49)] )$;
FF : trigsimp(F)$;

eg(kn) := cos(c*kn) + %i * sin(c*kn)$;
I : matrix([eg(0), eg(0), eg(0), eg(0),  eg(0), eg(0), eg(0), eg(0)],
           [eg(0), eg(1), eg(2), eg(3),  eg(4), eg(5), eg(6), eg(7)],
           [eg(0), eg(2), eg(4), eg(6),  eg(8), eg(10), eg(12), eg(14)],
           [eg(0), eg(3), eg(6), eg(9),  eg(12), eg(15), eg(18), eg(21)],
           [eg(0), eg(4), eg(8), eg(12),  eg(16), eg(20), eg(24), eg(28)],
           [eg(0), eg(5), eg(10), eg(15),  eg(20), eg(25), eg(30), eg(35)],
           [eg(0), eg(6), eg(12), eg(18),  eg(24), eg(30), eg(36), eg(42)],
           [eg(0), eg(7), eg(14), eg(21),  eg(28), eg(35), eg(42), eg(49)] )$;
II : trigsimp(I)$;

x : transpose(matrix([x0, x1, x2, x3,  0, 0, 0, 0]));   /* input x_n */

X : FF . x;                                             /* X = DFT(x) */

h : transpose(matrix([h0, h1, h2, h3,  0, 0, 0, 0]));   /* filter h_n */
H : FF.h;                                               /* H = DFT(h) */
Y : H*X;                                                /* Y = H*X  */
y : ratsimp(ratsimp((II.Y)/N));                         /* y = IDFT(Y) */

$y_n$ を、maxima は次のように出力する。
f:id:shimauma3san:20200716144613j:plain

重畳保留法(overlap save method)

重畳保留法は、 $x_n$ を、4点づつ次のように取る。
　_ $x0_n = (~~0,~~0, ~~0, ~~0, ~~x_0, x_1, x_2, x_3)$
　_ $x1_n = (x_0, x_1, x_2, x_3,~x_4, x_5, x_6, x_7)$
　_ $x2_n = (x_4, x_5, x_6, x_7, ~x_8, x_9, x_{10}, x_{11})$
　_ $x3_n = (x_8, x_9, ......)$
$h_n$ は同じである。
　 $h_n = (h_0, h_1, h_2, h_3, 0, 0 ,0, 0)$
それぞれの_ $x_n$ のフーリエ変換と、 $h_n$ のフーリエ変換を掛け算する。　その結果の逆フーリエ変換を取り、それぞれの $y_n$ の index 4 から 7 までを取って並べたものが、畳み込みになっている。maxima コードと結果を示す。結果の黄色の部分を使い、その他は捨てる。

　maxima コード

/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 19.05.7 ] */
/* [wxMaxima: input   start ] */
N : 8;
c : (2.0 * %pi / N)$;
ef(kn) := cos(c*kn) - %i * sin(c*kn)$;
F : matrix([ef(0), ef(0), ef(0), ef(0),  ef(0), ef(0), ef(0), ef(0)],
           [ef(0), ef(1), ef(2), ef(3),  ef(4), ef(5), ef(6), ef(7)],
           [ef(0), ef(2), ef(4), ef(6),  ef(8), ef(10), ef(12), ef(14)],
           [ef(0), ef(3), ef(6), ef(9),  ef(12), ef(15), ef(18), ef(21)],
           [ef(0), ef(4), ef(8), ef(12),  ef(16), ef(20), ef(24), ef(28)],
           [ef(0), ef(5), ef(10), ef(15),  ef(20), ef(25), ef(30), ef(35)],
           [ef(0), ef(6), ef(12), ef(18),  ef(24), ef(30), ef(36), ef(42)],
           [ef(0), ef(7), ef(14), ef(21),  ef(28), ef(35), ef(42), ef(49)] )$;
FF : trigsimp(F)$;

eg(kn) := cos(c*kn) + %i * sin(c*kn)$;
I : matrix([eg(0), eg(0), eg(0), eg(0),  eg(0), eg(0), eg(0), eg(0)],
           [eg(0), eg(1), eg(2), eg(3),  eg(4), eg(5), eg(6), eg(7)],
           [eg(0), eg(2), eg(4), eg(6),  eg(8), eg(10), eg(12), eg(14)],
           [eg(0), eg(3), eg(6), eg(9),  eg(12), eg(15), eg(18), eg(21)],
           [eg(0), eg(4), eg(8), eg(12),  eg(16), eg(20), eg(24), eg(28)],
           [eg(0), eg(5), eg(10), eg(15),  eg(20), eg(25), eg(30), eg(35)],
           [eg(0), eg(6), eg(12), eg(18),  eg(24), eg(30), eg(36), eg(42)],
           [eg(0), eg(7), eg(14), eg(21),  eg(28), eg(35), eg(42), eg(49)] )$;
II : trigsimp(I)$;

_x0 : transpose(matrix( [ 0,  0,   0,   0,   x0,  x1,  x2,  x3] )); /* input, a part of x_n */
_x1 : transpose(matrix( [x0, x1,  x2,  x3,   x4,  x5,  x6,  x7] ));
_x2 : transpose(matrix( [x4, x5,  x6,  x7,   x8,  x9, x10, x11] ));
_x3 : transpose(matrix( [x8, x9, x10, x11,  x12, x13, x14, x15] ));

X0 : FF._x0;                             /* X0 = DFT(_x0) */
X1 : FF._x1;
X2 : FF._x2;
X3 : FF._x3;

/****
x0_ : ratsimp(ratsimp(( II.X0 )/N));
x1_ : ratsimp(ratsimp(( II.X1 )/N));
x2_ : ratsimp(ratsimp(( II.X2 )/N));
x3_ : ratsimp(ratsimp(( II.X3 )/N));
*****/

_h : transpose(matrix( [h0, h1,h2, h3, 0, 0, 0, 0] ));  /* filter, a part of h_n */
H : FF._h;                           /* H = DFT(_h) */

/****
h_ : ratsimp(ratsimp((II.H)/N));
***/


Y0 : H*X0;                           /* Y0 = H * X0 */
Y1 : H*X1;
Y2 : H*X2;
Y3 : H*X3;


y0_ : ratsimp(ratsimp(( II.Y0 )/N));  /* y0_ = IDFT(Y0) */
y1_ : ratsimp(ratsimp(( II.Y1 )/N));
y2_ : ratsimp(ratsimp(( II.Y2 )/N));
y3_ : ratsimp(ratsimp(( II.Y3 )/N));



/* [wxMaxima: input   end   ] */



/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 19.05.7"$

maxima 出力

f:id:shimauma3san:20200716150208j:plain

重畳保留法でフィルタ系列をブロックに分割する方法(overlap save method for split block processing)

$x_n$ と $h_n$ を次のように取る。この方法では、 $h_n$ も分割する。
　_ $x0_n = (~~0,~~0, ~~0, ~~0, ~~x_0, x_1, x_2, x_3)$
　_ $x1_n = (x_0, x_1, x_2, x_3,~x_4, x_5, x_6, x_7)$
　_ $x2_n = (x_4, x_5, x_6, x_7, ~x_8, x_9, x_{10}, x_{11})$
　_ $x3_n = (x_8, x_9, ......)$

　_ $h0_n = (h_0, h_1, ~0, ~0, ~~0, ~0, ~0, ~0, ~0)$
　_ $h1_n = (h_2, h_3, ~0, ~0, ~~0, ~0, ~0, ~0, ~0)$

結果は、maxima 出力と、その次の図をみてほしい。

maxima コード

/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 19.05.7 ] */
/* [wxMaxima: input   start ] */
N : 8;
c : (2.0 * %pi / N)$;
ef(kn) := cos(c*kn) - %i * sin(c*kn)$;
F : matrix([ef(0), ef(0), ef(0), ef(0),  ef(0), ef(0), ef(0), ef(0)],
           [ef(0), ef(1), ef(2), ef(3),  ef(4), ef(5), ef(6), ef(7)],
           [ef(0), ef(2), ef(4), ef(6),  ef(8), ef(10), ef(12), ef(14)],
           [ef(0), ef(3), ef(6), ef(9),  ef(12), ef(15), ef(18), ef(21)],
           [ef(0), ef(4), ef(8), ef(12),  ef(16), ef(20), ef(24), ef(28)],
           [ef(0), ef(5), ef(10), ef(15),  ef(20), ef(25), ef(30), ef(35)],
           [ef(0), ef(6), ef(12), ef(18),  ef(24), ef(30), ef(36), ef(42)],
           [ef(0), ef(7), ef(14), ef(21),  ef(28), ef(35), ef(42), ef(49)] )$;
FF : trigsimp(F)$;

eg(kn) := cos(c*kn) + %i * sin(c*kn)$;
I : matrix([eg(0), eg(0), eg(0), eg(0),  eg(0), eg(0), eg(0), eg(0)],
           [eg(0), eg(1), eg(2), eg(3),  eg(4), eg(5), eg(6), eg(7)],
           [eg(0), eg(2), eg(4), eg(6),  eg(8), eg(10), eg(12), eg(14)],
           [eg(0), eg(3), eg(6), eg(9),  eg(12), eg(15), eg(18), eg(21)],
           [eg(0), eg(4), eg(8), eg(12),  eg(16), eg(20), eg(24), eg(28)],
           [eg(0), eg(5), eg(10), eg(15),  eg(20), eg(25), eg(30), eg(35)],
           [eg(0), eg(6), eg(12), eg(18),  eg(24), eg(30), eg(36), eg(42)],
           [eg(0), eg(7), eg(14), eg(21),  eg(28), eg(35), eg(42), eg(49)] )$;
II : trigsimp(I)$;

_x0 : transpose(matrix( [ 0,  0,   0,   0,   x0,  x1,  x2,  x3] ));  /* input, a part of x_n */
_x1 : transpose(matrix( [x0, x1,  x2,  x3,   x4,  x5,  x6,  x7] ));
_x2 : transpose(matrix( [x4, x5,  x6,  x7,   x8,  x9, x10, x11] ));
_x3 : transpose(matrix( [x8, x9, x10, x11,  x12, x13, x14, x15] ));

X0 : FF._x0;                                             /* X0 = DFT(_x0) */
X1 : FF._x1;
X2 : FF._x2;
X3 : FF._x3;

/****
x0_ : ratsimp(ratsimp(( II.X0 )/N));
x1_ : ratsimp(ratsimp(( II.X1 )/N));
x2_ : ratsimp(ratsimp(( II.X2 )/N));
x3_ : ratsimp(ratsimp(( II.X3 )/N));
*****/

/* _h : transpose(matrix( [h0, h1,h2, h3, 0, 0, 0, 0] )); */
_h0 : transpose(matrix( [h0, h1,  0,  0,  0, 0, 0, 0] ));  /* filter, a part of h_n */
_h1 : transpose(matrix( [h2, h3,  0,  0,  0, 0, 0, 0] ));


H0 : FF._h0;                                                /* H0 = DFT(_h0) */
H1 : FF._h1;


/*****
h0_ : ratsimp(ratsimp((II.H0)/N));
h1_ : ratsimp(ratsimp((II.H1)/N));
*****/


Y00 : H0*X0;                                               /* Y00 = H0 * X0 */
Y01 : H1*X0;

Y10 : H0*X1;
Y11 : H1*X1;

Y20 : H0*X2;
Y21 : H1*X2;

Y30 : H0*X3;
Y31 : H1*X3;

y00_ : ratsimp(ratsimp(( II.Y00 )/N));                      /* y00_ = IEFT(Y00) */ 
y01_ : ratsimp(ratsimp(( II.Y01 )/N));
y10_ : ratsimp(ratsimp(( II.Y10 )/N));
y11_ : ratsimp(ratsimp(( II.Y11 )/N));
y20_ : ratsimp(ratsimp(( II.Y20 )/N));
y21_ : ratsimp(ratsimp(( II.Y21 )/N));
y30_ : ratsimp(ratsimp(( II.Y30 )/N));
y31_ : ratsimp(ratsimp(( II.Y31 )/N));



/* [wxMaxima: input   end   ] */



/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 19.05.7"$

maxima 出力

f:id:shimauma3san:20200715183631j:plain

上記の結果を次のように加算する。　左の数字が $y_n$ の $n$ に対応する。 $y_3$ の場合は、その3行目を足せばよい。
（下の * は、掛け算を表す。)

	y00_           y01_            y10_           y11_         
	=============================================================
0	h0*x0
1	h0*x1 + h1*x0
2	h0*x2 + h1*x1  h2*x0
3	h0*x3 + h1*x2  h2*x1 + h3*x0
4	               h2*x2 + h3*x1   h0*x4 + h1*x3
5	               h2*x3 + h3*x2   h0*x5 + h1*x4
6	                               h0*x6 + h1*x5  h2*x4 + h3*x3
7	                               h0*x7 + h1*x6  h2*x5 + h3*x4
8	                                              h2*x6 + h3*x5  
9	                                              h2*x7 + h3*x6

重畳加算法(overlap add method)

重畳加算法は、 $x_n$ を次のように取る。
　_ $x0_n = (x_0, x_1, x_2, x_3, ~~0, ~0, ~0,~0)$
　_ $x1_n = (x_4, x_5, x_6, x_7, ~~0, ~0, ~0,~0)$
　_ $x2_n = (x_8, x_9, x_{10}, x_{11}, ~~0, ~0, ~0,~0)$
　_ $x4_n = (x_{12}, x_{13}, x_{14}, x_{15}, ~~0, ~0, ~0,~0)$
$h_n$ は、
　 $h_n = (h_0, h_1, h_2, h_3, 0, 0 ,0, 0)$
それぞれの_ $x_n$ のフーリエ変換と、 $h_n$ のフーリエ変換を掛け算する。　その結果の逆フーリエ変換を取り、隣の逆フーリエ変換の結果と加算する。maxima コードと maxima の計算結果と加算の図を次に示す。

maxima コード

/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 19.05.7 ] */
/* [wxMaxima: input   start ] */
N : 8;
c : (2.0 * %pi / N)$;
ef(kn) := cos(c*kn) - %i * sin(c*kn)$;
F : matrix([ef(0), ef(0), ef(0), ef(0),  ef(0), ef(0), ef(0), ef(0)],
           [ef(0), ef(1), ef(2), ef(3),  ef(4), ef(5), ef(6), ef(7)],
           [ef(0), ef(2), ef(4), ef(6),  ef(8), ef(10), ef(12), ef(14)],
           [ef(0), ef(3), ef(6), ef(9),  ef(12), ef(15), ef(18), ef(21)],
           [ef(0), ef(4), ef(8), ef(12),  ef(16), ef(20), ef(24), ef(28)],
           [ef(0), ef(5), ef(10), ef(15),  ef(20), ef(25), ef(30), ef(35)],
           [ef(0), ef(6), ef(12), ef(18),  ef(24), ef(30), ef(36), ef(42)],
           [ef(0), ef(7), ef(14), ef(21),  ef(28), ef(35), ef(42), ef(49)] )$;
FF : trigsimp(F)$;

eg(kn) := cos(c*kn) + %i * sin(c*kn)$;
I : matrix([eg(0), eg(0), eg(0), eg(0),  eg(0), eg(0), eg(0), eg(0)],
           [eg(0), eg(1), eg(2), eg(3),  eg(4), eg(5), eg(6), eg(7)],
           [eg(0), eg(2), eg(4), eg(6),  eg(8), eg(10), eg(12), eg(14)],
           [eg(0), eg(3), eg(6), eg(9),  eg(12), eg(15), eg(18), eg(21)],
           [eg(0), eg(4), eg(8), eg(12),  eg(16), eg(20), eg(24), eg(28)],
           [eg(0), eg(5), eg(10), eg(15),  eg(20), eg(25), eg(30), eg(35)],
           [eg(0), eg(6), eg(12), eg(18),  eg(24), eg(30), eg(36), eg(42)],
           [eg(0), eg(7), eg(14), eg(21),  eg(28), eg(35), eg(42), eg(49)] )$;
II : trigsimp(I)$;

_x0 : transpose(matrix( [x0, x1, x2, x3,  0,0,0,0] ));    /* input, a part of x_n */
_x1 : transpose(matrix( [x4, x5, x6, x7,  0,0,0,0] ));
_x2 : transpose(matrix( [x8, x9, x10, x11, 0,0,0,0] ));
_x3 : transpose(matrix( [x12, x13, x14, x15, 0,0,0,0] ));

X0 : FF._x0;                              /* X0 = DFT(_x0) */
X1 : FF._x1;
X2 : FF._x2;
X3 : FF._x3;

/****
x0_ : ratsimp(ratsimp(( II.X0 )/N));
x1_ : ratsimp(ratsimp(( II.X1 )/N));
x2_ : ratsimp(ratsimp(( II.X2 )/N));
x3_ : ratsimp(ratsimp(( II.X3 )/N));
*****/

_h : transpose(matrix( [h0, h1,h2, h3, 0, 0, 0, 0] )); /* filter h_n */
H : FF._h;                                             /* H = DFT(h_n) */

/****
h_ : ratsimp(ratsimp((II.H)/N));
***/


Y0 : H*X0;                                              /* Y0 = H*X0 */
Y1 : H*X1;
Y2 : H*X2;
Y3 : H*X3;


y0_ : ratsimp(ratsimp(( II.Y0 )/N));                    /* y0_ = IDFT(Y0) */
y1_ : ratsimp(ratsimp(( II.Y1 )/N));
y2_ : ratsimp(ratsimp(( II.Y2 )/N));
y3_ : ratsimp(ratsimp(( II.Y3 )/N));



/* [wxMaxima: input   end   ] */



/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 19.05.7"$

maxima 出力

f:id:shimauma3san:20200716151948j:plain

（下の * は、掛け算を表す。)

	y0_                            y1_                             y2_
	================================================================
0	h0*x0
1	h0*x1 + h1*x0
2	h0*x2 + h1*x1 + h2*x0
3	h0*x3 + h1*x2 + h2*x1 + h3*x0
4	h1*x3 + h2*x2 + h3*x1          h0*x4
5	h2*x3 + h3*x2                  h0*x5 + h1*x4
6	h3*x3                          h0*x6 + h1*x5 + h2*x4
7	0                              h0*x7 + h1*x6 + h2*x5 + h3*x4
8	                               h1*x7 + h2*x6 + h3*x5  
9	                               h2*x7 + h3*x6 
10	                               h3*x7                
11	                               0

重畳加算法でフィルタ系列をブロックに分割する方法
(overlap add method for split block processing)

$x_n$ を次のように分割する。
　_ $x0_n = (x_0, x_1, x_2, x_3, ~~0, ~0, ~0,~0)$
　_ $x1_n = (x_4, x_5, x_6, x_7, ~~0, ~0, ~0,~0)$
　_ $x2_n = (x_8, x_9, x_{10}, x_{11}, ~~0, ~0, ~0,~0)$
　_ $x4_n = (x_{12}, x_{13}, x_{14}, x_{15}, ~~0, ~0, ~0,~0)$
$h_n$ を次のように分割する。
　_ $h0_n = (h_0, h_1, ~0, ~0, ~~0, ~0, ~0, ~0)$
　_ $h1_n = (~~0,~~0, h_2, h_3, ~~0, ~0, ~0, ~0)$
それぞれの_ $x_n$ のフーリエ変換と、 $h_n$ のフーリエ変換を掛け算する。　その結果の逆フーリエ変換を取り、近くの逆フーリエ変換の結果と加算する。maxima コードと、maxima の計算結果と加算の図を次に示す。

maxima コード

/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 19.05.7 ] */
/* [wxMaxima: input   start ] */
N : 8;
c : (2.0 * %pi / N)$;
ef(kn) := cos(c*kn) - %i * sin(c*kn)$;
F : matrix([ef(0), ef(0), ef(0), ef(0),  ef(0), ef(0), ef(0), ef(0)],
           [ef(0), ef(1), ef(2), ef(3),  ef(4), ef(5), ef(6), ef(7)],
           [ef(0), ef(2), ef(4), ef(6),  ef(8), ef(10), ef(12), ef(14)],
           [ef(0), ef(3), ef(6), ef(9),  ef(12), ef(15), ef(18), ef(21)],
           [ef(0), ef(4), ef(8), ef(12),  ef(16), ef(20), ef(24), ef(28)],
           [ef(0), ef(5), ef(10), ef(15),  ef(20), ef(25), ef(30), ef(35)],
           [ef(0), ef(6), ef(12), ef(18),  ef(24), ef(30), ef(36), ef(42)],
           [ef(0), ef(7), ef(14), ef(21),  ef(28), ef(35), ef(42), ef(49)] )$;
FF : trigsimp(F)$;

eg(kn) := cos(c*kn) + %i * sin(c*kn)$;
I : matrix([eg(0), eg(0), eg(0), eg(0),  eg(0), eg(0), eg(0), eg(0)],
           [eg(0), eg(1), eg(2), eg(3),  eg(4), eg(5), eg(6), eg(7)],
           [eg(0), eg(2), eg(4), eg(6),  eg(8), eg(10), eg(12), eg(14)],
           [eg(0), eg(3), eg(6), eg(9),  eg(12), eg(15), eg(18), eg(21)],
           [eg(0), eg(4), eg(8), eg(12),  eg(16), eg(20), eg(24), eg(28)],
           [eg(0), eg(5), eg(10), eg(15),  eg(20), eg(25), eg(30), eg(35)],
           [eg(0), eg(6), eg(12), eg(18),  eg(24), eg(30), eg(36), eg(42)],
           [eg(0), eg(7), eg(14), eg(21),  eg(28), eg(35), eg(42), eg(49)] )$;
II : trigsimp(I)$;

_x0 : transpose(matrix( [x0, x1, x2, x3,  0,0,0,0] ));  /* input, a part of x_n */
_x1 : transpose(matrix( [x4, x5, x6, x7,  0,0,0,0] ));
_x2 : transpose(matrix( [x8, x9, x10, x11, 0,0,0,0] ));
_x3 : transpose(matrix( [x12, x13, x14, x15, 0,0,0,0] ));

X0 : FF._x0;                                            /* X0 = DFT(_x0) */
X1 : FF._x1;
X2 : FF._x2;
X3 : FF._x3;

/****
x0_ : ratsimp(ratsimp(( II.X0 )/N));
x1_ : ratsimp(ratsimp(( II.X1 )/N));
x2_ : ratsimp(ratsimp(( II.X2 )/N));
x3_ : ratsimp(ratsimp(( II.X3 )/N));
*****/

/* _h : transpose(matrix( [h0, h1,h2, h3, 0, 0, 0, 0] )); */
_h0 : transpose(matrix( [h0, h1,  0,  0,  0, 0, 0, 0] ));  /* filter, a part of h_n */
_h1 : transpose(matrix( [ 0,  0, h2, h3,  0, 0, 0, 0] ));


H0 : FF._h0;                                               /* H0 = DFT(_h0) */
H1 : FF._h1;


/*****
h0_ : ratsimp(ratsimp((II.H0)/N));
h1_ : ratsimp(ratsimp((II.H1)/N));
*****/


Y00 : H0*X0;                                               /* Y00 = H0 * X0 */
Y01 : H1*X0;

Y10 : H0*X1;
Y11 : H1*X1;

Y20 : H0*X2;
Y21 : H1*X2;

Y30 : H0*X3;
Y31 : H1*X3;

y00_ : ratsimp(ratsimp(( II.Y00 )/N));                     /* y00_ = IDFT(Y00) */ 
y01_ : ratsimp(ratsimp(( II.Y01 )/N));
y10_ : ratsimp(ratsimp(( II.Y10 )/N));
y11_ : ratsimp(ratsimp(( II.Y11 )/N));
y20_ : ratsimp(ratsimp(( II.Y20 )/N));
y21_ : ratsimp(ratsimp(( II.Y21 )/N));
y30_ : ratsimp(ratsimp(( II.Y30 )/N));
y31_ : ratsimp(ratsimp(( II.Y31 )/N));

maxima 出力

f:id:shimauma3san:20200716115259j:plain

（下の * は、掛け算を表す。)

	y00_           y01_           y10_           y11_            y20_
	=================================================================
0	h0*x0          0
1	h0*x1 + h1*x0  0
2	h0*x2 + h1*x1  h2*x0
3	h0*x3 + h1*x2  h2*x1 + h3*x0
4	h1*x3          h2*x2 + h3*x1  h0*x4          0
5	0              h2*x3 + h3*x2  h0*x5 + h1*x4  0
6	0              h3*x3          h0*x6 + h1*x5  h2*x4
7	0              0              h0*x7 + h1*x6  h2*x5 + h3*x4
8	                              h1*x7          h2*x6 + h3*x5
9	                              0              h2*x7 + h3*x6
10	                              0              h3*x7           ??
11	                              0              0               ??

その他

$x_n$ と $h_n$ の分割方法は他にもある。

shimauma3sanのブログ

オペアンプの動き、言葉で説明

Windows 11 で削除できないファイルを削除できた方法。

重畳保留法(overlap save method)と重畳加算法(overlap add method) とそれぞれの split block processing の帰納法的説明

畳み込み、フーリエ変換の関係式の整理

畳み込みをフーリエ変換、フーリエ逆変換を使って計算する

重畳保留法(overlap save method)

重畳保留法でフィルタ系列をブロックに分割する方法(overlap save method for split block processing)

重畳加算法(overlap add method)

重畳加算法でフィルタ系列をブロックに分割する方法
(overlap add method for split block processing)

その他

定点観測 5月21日今日のコロナ君

定点観測 5月20日今日のコロナ君

定点観測 5月19日今日のコロナ君

定点観測 5月18日今日のコロナ君

畳み込み、フーリエ変換の関係式の整理

畳み込みをフーリエ変換、フーリエ逆変換を使って計算する

重畳保留法(overlap save method)

重畳保留法でフィルタ系列をブロックに分割する方法(overlap save method for split block processing)

重畳加算法(overlap add method)

重畳加算法でフィルタ系列をブロックに分割する方法 (overlap add method for split block processing)

その他

重畳加算法でフィルタ系列をブロックに分割する方法
(overlap add method for split block processing)