テキスト

生産システム工学専攻総合演習 「信号の自己相関とパワースペクト ルの計算」
実施日: 2008/1/18 、1/25 、2/1 担当: 東海林
1
CSV ファイル
csv ファイル ( comma separated value : カンマ区切り値 ) はテーブルデータをテキスト形式で表現するファイル形式
であり、表計算ソフトなどで広く使用されている。テーブルの各行は 1 行のテキストで表され、各行の列 (フィールド と
呼ぶ) は数値、もしくは文字列で表され 、各フィールド の間はカンマ (,) で区切られる。例えば


1.0 2.0 3.0


4.0 5.0 6.0
7.0 8.0 9.0
というテーブルを csv ファイルで表すと
1.0, 2.0, 3.0
4.0, 5.0, 6.0
7.0, 8.0, 9.0
となる。
問題 1 任意のプログラム言語で一次元デジタル信号列 x[i], i = 0,1,2,· · · を csv ファイルに出力する関数を作成せよ。
例えば C 言語では
void print csv( char* filename, double* x, int N );
という関数になる。ここで filename は csv ファイルのファイル名、N は配列の大きさである。
なお csv ファイルの各行には時刻と信号値を入れること。例えば x[0] = 1.5 、x[1] = 2.1 、x[2] = 3.9 のとき、csv ファ
イルの中身は
0, 1.5
1, 2.1
2, 3.9
となる。
問題 2 任意のプログラム言語で時刻 i における値が
x[i] = 0.1 · cos
µ
2πi π
−
10
2
¶
+ 0.2 · cos
µ
2πi π
−
4
4
¶
+ 0.3 · cos
µ
2πi π
+
3
6
¶
(1)
となる一次元デジタル信号列を出力する関数を作成せよ。例えば C 言語では
double func x( int i )
という関数になる。また時刻 0 から 59 までの信号列を csv ファイルに出力して任意のグラフ作成ツールを用いてグラ
フ表現せよ。
2
離散フーリエ変換 (DFT)
長さ N の一次元デジタル信号列 x[i], (i = 0, · · · , N − 1) が与えられたとき、離散フーリエ変換 (DFT) は次式により
定義される。
ck =
µ
¶
N−1
2π
1 X
·
x[i] · exp −j ki , k = 0, · · · , N − 1
N i=0
N
1
(2)
ここで ck , (k = 0, · · · , N − 1) は k 番目の離散フーリエ係数である。(2) はオイラー公式から実数行列演算に置き換える
ことができて、ck を
ck = ak + j · bk
(3)
と実数部と虚数部に分けたとき、ak , bk は次の行列演算により求められる。
A=
1
· Dc · X
N
B=−
(4)
1
· Ds · X
N
(5)
ここで A 、B 、X は N × 1 行列であり、それぞれ
A = [a0 , · · · , aN−1 ]t
(6)
B = [b0 , · · · , bN−1 ]t
(7)
X = [x[0], · · · , x[N − 1]]t
(8)
である。また N × N の DFT 行列 Dc は k を行番号、i を列番号として次のように定義される対称行列である。


¡
¢

Dc =  cos 2π
N ki


(9)
すなわち、Dc の k 行 i 列目の値が cos( 2π
N ki) となる。同様に Ds も次のように定義される。

¢
¡

Ds =  sin 2π
N ki



(10)
また、離散フーリエ係数 ck のノルムの二乗
Pk = |ck |2 = a2k + b2k
(11)
をパワースペクトルと呼ぶ。
問題 3 任意のプログラム言語で行列の積 C = AB をおこなう関数を作成せよ。ただし 、A の行数と列数を L と M 、
B の行数と列数を M 、N 、C の行数と列数を L と N とする。例えば C 言語では
void prod mat( double** A, double** B, double** C, int L, int M, int N );
となる。なお、C 言語で L × M の行列 (2 次元配列) のメモリを確保／解放するには次のような関数を用いる。
double** make_mat( int L, int M )
{
double** A;
A = (double**)malloc( sizeof( double* ) * L );
A[0] = (double*)malloc( sizeof( double ) * L * M );
for( int i = 0; i < L; ++i ) A[i] = A[0] + i * M;
return A;
}
2
void free_mat( double** A )
{
free( A[0] );
free( A );
}
問題 4 任意のプログラム言語で行列の積 (prod mat 関数) を用いて DFT を計算する関数を作成せよ。例えば C 言語
では
void DFT( double **X, double **A, double **B, int N );
となる。また、デジタル信号列 (1) のパワースペクトルを csv に出力してグラフで表現し 、元の信号列のグラフと比較
せよ。ただし N = 60 とし 、問題 1 で作成した csv 出力関数を行列に対応させて用いること。
3
自己相関関数
長さ N の一次元デジタルデジタル信号列 x[i], (i = 0, · · · , N − 1) が周期的であると仮定する。すなわち x[N] =
x[0], x[N + 1] = x[1], · · · とする。これを記号 mod N を用いて x[i] = x[i + n]mod
φxx (n), (n = 0, · · · , N − 1) は次式で定義される。
φxx (n) =
N−1
1 X
x[i + n]mod
N i=0
N
N
と記述するとき、自己相関関数
· x[i]
(12)
問題 5 任意のプログラム言語で、行列の積 (prod mat 関数) を用いて自己相関関数値を計算する関数を作成せよ。例
えば C 言語では
void phi( double** X, double** P, int N );
となる。ここで P = [φxx (0), · · · , φxx (N − 1)]t である。またデジタル信号列 (1) の 0 から 59 番までの自己相関関数の
出力列を csv に出力しグラフで表現せよ。
問題 6 任意のデジタル信号列のパワースペクトルと自己相関関数の DFT の実数成分は一致する。デジタル信号列 (1)
の自己相関関数の DFT の実数成分を csv に出力しグラフで表現してこのことを確認せよ。
問題 7 自己相関関数はどのような性質を持ち、さらにその性質を利用してどのような用途に用いられているのか説明
せよ。
問題 8 自己相関関数の DFT とパワースペクトルが一致することにはどのような利点があるか説明せよ。
問題 9 自己相関関数をあなたが現在おこなっている研究に役立たせるにはど うすれば良いか説明せよ。
3