上にカーソルを合わせ，マウスの左ボタンをクリックする（以下，クリックと略す）と，ターミナルが起動する。医学データ処理概論
２． 時系列信号処理の基礎
(1) Octave によるプログラミング
今回の実習は，プログラミング言語octaveを用いて行ないます。
最近，科学技術計算にプログラミング言語MATLABがよく用いられます。Octaveはそのクローンです。
MATLABやOctaveの特徴は，ベクトルや行列の演算を計算式どおりに書けることです。
（c，PASCAL，Javaなどの言語では，行列の演算はやっかいです。）
MATLABは有料ですが，同様に行列をそのまま扱えるフリーの言語にOctave，Scilabな
どがあります。
ここではOctaveを用います。Octave
のマニュアルはここ。
Octaveのプログラムの多くは，そのままか，あるいは少し修正すれば，MATLAB上で実行できます。
MATLABは，本郷の医学図書館の計算機システムで利用可能です。
MRIなどの医療機器の制御用計算機ではunixがよく使われており，unixを知っておくと便利ですので，本実習はunix上で行ないます。
１．unixの立ち上げ
Dock（下または横にあるアイコンが並んでいるバー）内にあるターミナル・アイコン
上にカーソルを合わせ，マウスの左ボタンをクリックする（以下，クリックと略す）と，ターミナルが起動する。
２．Octaveの起動と終了
キーボードを英数モードにして、ターミナル窓に，次のようにタイプする。（赤字のみ）（g12345の部分は，ユー
ザにより異なる）
（行の最後に改行(Return)を入力すること。）
g12345$>
octave
GNU Octave, version 3.0.0
Copyright (C) 2007 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTIBILITY or
FITNESS FOR A PARTICULAR PURPOSE. For details, type `warranty'.
Octave was configured for "i386-apple-darwin9.2.2".
Additional information about Octave is available at
http://www.octave.org.
Please contribute if you find this software useful.
For more information, visit http://www.octave.org/help-wanted.html
Report bugs to <bug@octave.org> (but first, please read
http://www.octave.org/bugs.html to learn how to write a helpful report).
For information about changes from previous versions, type `news'.
octave:#n> ・・・・・ %
以降の演習問題を入力する（#nは
計算機から出力される行番号の数字）
・ ・・・・・
・ ・・・・・
・ ・・・・・
octave:#n> exit
% Octaveの終了
g12345$>
以下を単に読むだけではなく，実際にキーボードを叩いて実行してみましょう。
３．Octaveのプログラム要素
赤字の部分だけタイプ入力すること。（行
の最後に改行(Return)を入力すること。）
以下の例にある % の右側の文字列は説明文（コメント）なので，入力する必要はない。
Octaveのプロンプト（入力催促）： octave:#n> が計算機から出力さ
れるが，
#n の部分は，場合により
異なるので気にしないこと。
なお，行数の多い時に，次のような行が現れることがあるが，f キーを叩くと次の数行を見ることができる。
q キーを叩くと表示が終わって次のコマンド入力ができる。
-- less -- (f)orward, (b)ack,
(q)uit
１）代入文
演算記号
| 加算 | 減算 | 乗算 | 除算 | べき乗 |
| + | - | * | / | .^ |
スカラーの演算
octave:1> x=10*3+5
% 変数 x に 10*3+5 の計算結果を代入する。（＝の左辺の変数に，右辺の計算結果を代入する。）
x =
35 %
計算結果が自動的に表示される。
octave:2> y=9*5+3;
% ； を付けると，計算結果を表示しない。（変数 y には計算結果が代入される。）
octave:3>
y % ｙ の計算結果を表示させたい時は，ｙ とタイプする。
y = 48
octave:4>
z=x+y % 変数 x と y の内容を加算して z に代入する。
z = 83
ベクトルの演算
octave:5> v=1:10
% 変数 v に，（初期値）：（最終値）を代入する。
v =
1 2 3 4 5 6 7 8 9 10
octave:6> u=2:3:20
% 変数 u に，（初期値）：（ステップ）：（最終値）を代入する。
u
=
初期値は０や負数でもよい。（ステップ）を省略すると１となる。
2 5 8 11 14 17 20
octave:7> v=v*0.2
% 現在の v の各要素に 0.2 を掛けて，結果を
v=
元の v の各要素と入れ替える。
Columns 1 through 8:
0.20000 0.40000 0.60000 0.80000 1.00000 1.20000 1.40000 1.60000
Columns 9 and 10:
1.80000 2.00000
octave:8> w=u*2+1
w =
5 11 17 23 29 35 41
octave:9> u(3)
% u の３番目の要素。u(i) i 番目の要素。
ans = 8
行列の入力と演算
ctave:10>
A=[1,2,3;4,5,6;7,8,9] %
行列の入力。要素間は，カンマ( , ) またはスペースで区切る。
A = % 行の区切りは，セミコロン( ; ) または改行(Return)。
1 2 3
4 5 6
7 8 9
octave:11> B=[ 1 2 3
% ３行に分けて入力してもよい。
> 4 5 6
> 7 8 9 ]
B =
1 2 3
4 5 6
7 8 9
octave:12> A(2,3)
% Aの２行目，第３列の要素。A(i, j) i行，j列の要素を表示する。
ans = 6
octave:13> C=A*0.1
C =
0.10000 0.20000 0.30000
0.40000 0.50000 0.60000
0.70000 0.80000 0.90000
octave:14> C(2,:)=0
% 行列Cの２行目の全ての要素（：で指定）に０を代入する。
C =
0.10000 0.20000 0.30000
0.00000 0.00000 0.00000
0.70000 0.80000 0.90000
octave:15> Z=zeros(3,4) % ゼロ行列 zeros は組込み関数（Octaveで用意されている関数）。
Z =
0 0 0 0
0 0 0 0
0 0 0 0
行列の加算と減算
octave:16> A2=[2 0 0;0 2 0;0 0 2]
A2 =
2 0 0
0 2 0
0 0 2
octave:17> B3=[0 0 3;0 3 0;3 0 0]
B3 =
0 0 3
0 3 0
3 0 0
octave:18> AplusB=A2+B3
AplusB =
2 0 3
0 5 0
3 0 2
octave:19> AsubB=A2-B3
AsubB =
2 0 -3
0 -1 0
-3 0 2
行列の乗算と除算
octave:20> AmultB=A2*B3
AmultB =
0 0 6
0 6 0
6 0 0
octave:21> AdivB=A2/B3
AdivB =
0.00000 0.00000 0.66667
0.00000 0.66667 0.00000
0.66667 0.00000 0.00000
行列の要素ごとの乗算と除算（ドット . に注意，後で使うので覚えておくこと）
octave:22> AemultB=A2 .* B3
AemultB =
0 0 0
0 6 0
0 0 0
octave:23> A2edivB3=A2 ./ B3
A2edivB3 =
Inf NaN 0.00000
NaN 0.66667 NaN % NaN : Not-a-Number
定義できない結果
0.00000 NaN Inf % Inf : 無限大
２）繰り返し
octave:24> for i=1:10 % x の各要素に i の２乗を代入する例題。
> x(i)=i*i;
% x の要素の指定は１番目からでないといけない。
>
end
% for文とend文に挟まれた演算式を10回繰り返す。
octave:25> x
x =
1
4
9
16
25
36
49
64
81
100
octave:26> for
j=1:10;y(j)=j^2;end; %
一行に書く事もできる。
octave:27> y
y =
1
4
9
16
25
36
49
64
81
100
行列の転置
octave:28> b=(1:10)'
% ' は，転置を表す。
b =
1
2
3
4
5
6
7
8
9
10
３）端末からの入出力
x=input('Please input a value:
') % プロンプト（催促）の文字列を ’
と ’ で囲む。
%
” と ” でも
よい。変数 x にはタイプ入力した値が代入される。
a=input('Please input string: ', 's') % 変数 a にはタイプ入力した文字列が代入される。
例
octave:29> a=input('input string: ','s')
input string: abc
a = abc
disp('The value is : ') % The value is と端末に表示する。
disp(x) % 変数 x の値を表示する。
４）条件判断（if 文）
octave:30> a=input('x:
') % 端末からの入力を要求し，値を a
に代入する。
x: 7
a = 7
octave:31> if rem(a,2)==0 % a を 2 で割った余りが 0 ならば次行を実行する。
> disp('a is even') % rem() は剰余を求める組込み関数。
>
else
% そうでないならば，次行を実行する。
> disp('a is odd')
> end
a is odd
if 文の構文
if 条件式
条件を満たす場合に実行する文
else
条件を満たさない場合に実行する文
end
if 文の条件に用いられる論理演算記号
| > | より大きい | < | より小さい |
| >= | 以上 | <= | 以下 |
| == | 等しい | ~= | 等しくない |
| & | および (and) | | | または (or) |
| ! | でない |
５）グラフの表示
octave:32> x=1:10;
octave:33> y=x .^ 2;
octave:34> plot(x, y)
% (x, y) グラフをプロットする。
６) 変数の表示と消去
変数を使用していくうちに，使用済みの不要な変数がでてくる。使われた変数の一覧は，whos
-variables
で見ることができる。不要な変数を消去するには，clear を使う。
octave:35> whos -variables % ユーザーが作成した変数のリストを表示する
octave:36> clear c
% 変数 c を消去
octave:37> clear
% ユーザーが作成した変数をすべて消去
４．テキスト・エディタ ( mi ) の使用
Octaveは，入力されたコマンドを即時に解釈しながら実行するインタプリタ方式
の言語です。
インタプリタ方式は，プログラムを修正しながら即実行できるので便利ですが，次に述べるコンパイラ方式の言語に比べて実行時間が長いと
いう弱点があります。
これに対し，コンパイラ方式の言語であるｃやpascalでは，プログ
ラムを計算機に最適な実行形式に予め翻訳（コンパイルという）しておくので実行時間が短く，大規模な計算を実行するプログラムや，固定された作業を実行す
る プログラムに適します。
Octaveのプログラムで何回も同じ仕事をさせたい時，同じコマンドを始めからタイプし直すのは時間の無駄です。
テキスト・エディタを使って， コマンド列のテキスト・ファイルを作っておき，Octaveの実行時にそれを呼び出して使うことができます。
ここでは，テキスト・エディタとして，mi 
を使う ことにします。
（ mi は初めから入っているので上のm i をクリックする必要は無い）
最初アイコンが入っているDockは、画面の下側に大きく表示されていますが、邪魔なので画面の
右端に小さく表示されるようにします。この操作は，最初に１回だけ行えばよい。
Dockを小さくして画面の右側に配置する。
アップル・マーク→Dock→Dock環境設定
だいたい下のように設定する。
を押して
Dock環境設定窓を閉じる。
５．作業フォルダの作成
１）Finderアイコン
（通常一番上にある）をクリックしてFinderウインドウを表示させる。
各自のユーザ名（ca12345な
ど）がFinderウインドウの最上部に表示されるようにする。
２）Finderウインドウの上端をクリックすると，画面一番上のバーがFinder用になる。
３）ファイル→新規フォルダ
４）「名称未設定フォルダ」の文字
をクリックして適当な名前に変更する（この例ではSIGNAL_PROC）。
６．テキスト・エディタ mi
を使ったoctaveプログラムの実行
Dock内にある mi アイ
コン をクリックする。上の条件判断のところで試したプログラムを
mi に入力して，handan.m
というテキスト・ファイルを作る。
テキストエンコーディングが UTF-8, 改行コードが LF(UNIX) であることを確認して、
ファイル →別名で保存を選
ぶ。
新規フォルダを作成後１回目は，場所に新規フォルダ名を表示できないので，
を押して場所窓
(Finder窓と同様の窓)を開き，新規フオルダ名を選択する。
ここで一旦octaveを終了して，作業フォルダ(unixの場合はディレクトリdirectoryという）に移る。
ターミナル窓内で以下の赤字のコマンドをタイプする。（右側の緑の
文字は説明）
ca12345$ cd
SIGNAL_PROC (change directory)
SIGNAL_PROCは自分で作った作業フォルダ名
ca12345$ pwd (print working directory) 現在のdirectoryを表示
ca12345$ octave 再びoctaveを実行
handan.m という名前で保存された Octave のプログラムを実行する。Octave を実行して，
octave:1>
handan % [注
意] .m
の部分はタイプしないこと。
x: 13
a = 13
a is odd
octave:2> handan
x: 16
a = 16
a is even
octave:3>
以降，Octaveが動作しているターミナルと，テキスト・エディ
タ mi
とを同時に開いておき， テキストを修正したら別名保存（同じファイルを修正した場合は保存）して，即Octave
で実行すれば，効率よく実習が行なえる。
ファイルを修正して保存し再実行する場合，ターミナル窓内にカーソルを置いてクリックし，キー・ボードの
[上向き矢印・キー↑] を押すと，前にタイプしたコマンドの文字列が現れるので，[リターン・キー] を押すと再
実行される。
(2) 信号処理入門
１．アナログとディジタル (Analog and Digital)
医学で用いられるデータには，時間的に変化するものがある。脳波，筋電図，心電図，音声波形，など
いろいろある。
これらのデータは，通常，電気信号に変換されてデータ処理される。
物理現象の場合には，物理量を電気信号に変換して用いる。この変換器を一般に，「トラン
スデューサ」という。
例えば，音声波形では，音波を，「マイクロホン」という「トランスデューサ」を通して，電気信号に変換する。
これらの電気信号は，アナログ・データと呼ばれ，時間的に連続して変化する量である。
これを電子計算機で処理する場合，連続量のアナログ・データから離散量のディジタル・データに変換する必要がある（ＡＤ変換）。
これには，A/D変換器（アナログ-ディジタル変換器，A/D Converter) が用いられる。
A/D変換器を用いて，アナログ・データから一定時間ごと（周期的）に信号値を読み取り，その値を時間的に順に並べた時系列 (Time series) データに変換す
る。この操作を，標本化 (Sampling) という。
また，アナログ・データは連続量であるので，標本化と同時に，信号値を計算機で扱える最小単位（の値）の整数倍の数値に変換しなければならない。この変換
を量子化 (Quantization) という。
１）標本化（サンプリング Sampling）
標本化は，通常，一定時間間隔（周期）ごとに行なうが，その周期を，標本化周期 (Sampling period)，その
逆数を標本化周波数 (Sampling frequency) という。
標本化周期は，アナログ・データの変動周期との時間関係を考えて選ぶ必要がある。例えば，１時間周期で変動する現象を１秒ごとに標本化すると，データが
多過ぎて現象の変動傾向を捕えにくいし，同じ現象を１日周期で標本化すると，変動は捕えられない。したがって，標本化周期の選び方には，何らかの適当な基
準があるはずである。この時間関係を規定するサンプリング定理（標本化定理）が知られている。
サンプリング定理（標本化定理）
アナログ・データは，その信号中に含まれる成分の最高周波数の，２倍以上高い標本化周波数
で標本化しなければならない。
 |
 |
○印は標本点 （サンプル・ポイント）を示す。 標本化周波数が低すぎる場合には，実際とは異なる周波数のディジタル・データ（--印）として読込まれているのがわかる。 |
|
十分な標本化周波数で標本化 |
標本化周波数が低すぎる場合 |
２）量子化(Quantization)
 |
 |
量子化では，入力信号を何段階の値に分割するかが問題となる。計算機内部では通常２進法でデータが取扱われるので，分割数は普通２のべき乗の値が とられる。この何乗かを量子化ビット(bit)数という。左図に示すように，量子化は非線形の変換である。例えば，右図のように，入力信号の最大値が± １Vの場合，この信号を５bit（32段階）で量子化すると，量子化ステップ0.0625V以下の信号の変動は捕らえられない。このように， もとのアナログ源波形と出力波形との間には，中図のように常に量子化誤差が発生する。量子化誤差を小さくするには，量子化ビット数を増やせばよい。通常， 物理現象の観測データの有効数字は３桁（10進）程度であるので，量子化ビット数は10bit以上必要である（２10= 1024）。パソコンでは 16bitが用いられる場合が多い （最近のパソコンのオーディオ入出力では24bitが用いられることもある）。
3）複数byteデータとエンディアン (endian)
計算機のメモリはbyte単位でaddressが振られている（要するに０から何番目のbyteであるかということ）。数値を2byte以上で表
現する場合、上位の桁(byte)が若い（０に近い）addressになるようにして順に並べる方法と、下位の桁(byte)が若いaddressになる
ようにして順に並べる方法とがあり得る。このbyteを並べる順序をバイト・オーダーと言い、計算機の機種（アーキテクチャー）により異なる。それぞれ
Big Endian, Little Endian という。Endian
の異なる機種間でデータをやり取りする場合には注意を要する。（バイト・オーダーを変換する計算機プログラムを用いて変換する必要がある。）例えば数値が
4byteで表現されている場合は、下図のようになる。
| Big Endian |
|
Mac(PowerPC, 680x0, unix(Sun, SGI, HP) | ||||
| Little Endian |
|
Windows, Mac(Intel) |
２．正弦波の重畳と周波数スペ クトル
どのような信号波形でも，種々の正弦波（sin, cos 関数）の足し合わせで作り出せること，また，どのような信号波形でも，種々の正弦波に分解できることが知られている。ここでは，正弦波の足し合わせについ て Octave を使って演習する。自分でタイプして試してみること。
１） ３つの正弦波の加算
行列演算を使ったMATLAB/Octave方式のプログラム
|
Y(i, :) 中の ： は，行列Ｙの i
行目の「すべての要素」を意味する記号である。（ワイルドカードという。） plot(x, Y(1,:),x,Y(2,:),x,Y(3,:), x, Sum) は，(x, Y) と (x, Sum) のプロットを重ねて表示する。（ x, y のペアを追加していけば，いくつも重ねて表示できる。） pi は予約定数であり，自動 的に3.141592・・・に置き換えられる。 |
位相
 |
３つの正弦波の位相を変えてみる。 sin関数の引数（関数に与えるパラメタ）中に位相項 p(i) * pi を追加する。 ｐ:３つの正弦波の位相の係数が入るベクトル。0π，π，π/2と順に位相が変わる。 ３つの正弦波の周波数は同じでも，位相が変わると，加算した波形は異なる。 |
 |
方形波 (Square wave) の合成
 |
上の「行列演算を使ったMATLAB/Octave方式のプログラム」で，適当な振幅係数を掛け合せた周
波数の異なる正弦波の加算数を増やしていくと，方形波を合成することができる。
方形波は，基本波の奇数倍の周波数をもつ正弦波を加算することにより合成できる。個々の振幅係数は， |
 |
演習課題
[1]
上に述べた，１５個の正弦波成分を加算して方形波のグラフをプロットするプログラムを作ること。
ヒント：方形波 の合成のプログラムと３つの正弦波の加算のプログラムの一部を組み合わせる必要があります。
[2] 三角波も，奇数
倍の周波数の正弦波の加算により合成できる。個々の振幅係数は，
2n-1倍の周波数成分の振幅＝
8・(-1)(n-1)/（π2・(2n-1)2）
で求まる。１５個の正弦波成分を加算して，下の三角波のグラフをプロットするプログラムを作ること。
（π2は，プログラムでは，べき乗を使うより pi*pi と書いた方が簡単である。）
２）加算波形の各正弦波成分の大きさを求める
加算された信号波形は，種々の正弦波成分に分解でき，個々の
成分の大きさを求めることができる。これは，離散フーリエ変換という方法を用いて行なうことができる。通常，その高速計算法である FFT(Fast Fourier
Transform) を用いて計算される。FFTは，Octave の組込み関数である。
個々の周波数成分の大きさを周波数スペクトル(Frequency Spectrum）という。
 |
動作しない場合は、改行コードを
LF(UNIX) にしてみること。
３つの正弦波（1Hz, 3Hz, 5Hz）を加算した波形の周波数スペクトルを求める。サンプリング周波数は51.2Hz，FFT点数は512点である。FFT点数には２のべき乗の値が用 いられることが多い。[周波数分解能]＝[サンプリング周波数]÷[FFT点数]。後で説明するが，FFTの結果には，サンプリング周波数を周期とするス ペクトルの折返しが含まれているので，周波数軸の下半分（25.6Hzまで）の結果だけを用いる。 clear は，Ovtaveを初期状態にもどす。 |
 |
演習課題
[3] ３つの位相の異なる正弦波を加算した波形の周波数スペクトルをプロットすること。
[4] １５個の正弦波成分を加算して作った方形波の周波数スペクトルをプロットすること。
[5] １５個の正弦波成分を加算し て作った三角波の周波数スペクトルをプロットすること。
３）標本化による折返し歪 (aliasing)
演習課題[4]の結果のプロットは下図のようになる。しかし，この図を見るとおかし い所がある。それは，下から１２番目以上の成分が，周波数軸（横軸）上で等間隔になっていないことである。どうしてであろうか？
周波数軸（横軸）を，サンプリング周波数（51.2Hz）まで拡げてプロットしてみると下図のようになる。１５個の正弦波の各成分（元の成 分）が，周波数0Hzを原点として周波数の高い方に順に並んで表示されているのと同時に，サンプリング周波数51.2Hzを原点として，それらを折返した 成分（折返し成分）が，周波数の低いほうに順に並んで重なって表示されているいるのが分かる。上図で不等間隔に表われている周波数成分は，周波数の高い方 から折返してきた，14番目と15番目の成分である。14番目の成分の周波数は，(2×14-1)=27Hz，15番目の成分の周波数は，(2×15- 1)=29Hzであり，サンプリング周波数（51.2Hz）の半分の25.6Hz以上の周波数成分である。このように，入力データにサンプリング周波数の 1/2以上の周波数成分が含まれていると，正しい標本化が行なわれない。これを，標本化による折返し歪 (aliasing エイリアジング) という。
４）折返し歪の防止
折返し歪を防ぐ方法は，2種類ある。
1番目は，入力データにサンプリング周波数の1/2以上の周波数成分が含まれないように，サンプリング周波数を高くすることである。下の図
は，上図のサンプリング周波数を102.4Hzにしたものである。元の周波数成分と折返し成分とは，全く重なっておらず，折返し歪のないことが分かる。
2番目は，入力データをＡＤ変換する前に，サンプリング周波数の1/2以上の周波数成分を，ローパス・フィルタ (Low pass Filter) を用いて取り除いておく方法である。この目的のローパス・フィルタは，折返し防止フィルタ (anti-aliasing filter) と呼ばれ，アナログ信号を扱うアナログ・フィルタである。折返し防止フィルタについては後で述べる。
３．時系列信号処理
１）ディジタル・フィルタ
標本化されたディジタル信号を扱うディジタル・フィルタというものもある。その入力信号の標本化には，やはり，折返し防止フィルタが必要である。 アナログ・フィルタ，ディジタル・フィルタを問わず，一般に，信号処理においては，フィルタは特定の周波数成分を通過または除去（阻止）するために用いら れる。このようなフィルタには，下図のようなものがある。名称には、通過する信号に着目した場合と、除去する信号に着目する場合のものがある。（図中の阻 止 域における櫛型カーブは一部のディジタル・フィルタに特有の特性であり，アナログ・フィルタでは現れない。）
 |
低域通過フィルタ ローパス・フィルタ LPF (Low Pass Filter) 高域除去フィルタ ハイカット・フィルタ |
 |
高域通過フィルタ ハイパス・フィルタ HPF (High Pass Filter) 低域除去フィルタ ローカット・フィルタ |
 |
帯域通過フィルタ バンドパス・フィルタ BPF (Band Pass Filter) |
 |
帯域除去フィルタ 帯域阻止フィルタ バンドエリミネーション・フィルタ BEF (Band Elimination Filter) |
デシベル（dB）
上図の縦軸の目盛りは，dB（デシベル）である。フィルタの特性や周波数スペクトルの表示には，一般にデシベル(dB)表示が使われる。ある入力波形を
フィルタに通した周波数スペクトルは，入力波形の周波数スペクトルに，フィルタ特性を掛け合わせた（正確には，たたみ込み convolution
という）ものとなる。デシベル(dB)表示を用いた場合，各周波数の成分の大きさは，その周波数における周波数スペクトルとフィルタ特性のデシベル値の和
で計算できる。この場合，フィルタの特性は，その入力と出力の関係を表した関数と考えることができ，工学的には伝達関数と呼ばれる。
臨床や研究の場においても周波数スペクトルはデシベル表示されるのが一般的である。例えば純音聴覚検査では、種々の周波数の正弦波（純音）の音圧がどの
程度大きければ聞こえるかを検査する。健常な聴覚を持つ若年者が聞き取れる音圧の範囲は，振幅で約20μPaから約20,000,000μPa（Paは音
圧の単位パスカル=N/m2、μは10-6）もある。振幅をそのまま表示すると８桁もの数値を書くことにな
る。これに対してデシベル（dB）表示では、 20log10(P/P0)
というデシベル値を使う。これを使うと0〜120dBまでの範囲で表現できる。物理学的な音圧表示（dB SPL）（sound pressure
level）では，P0=20μPa（健聴者の標準的な最小可聴値）を用いる。(1Pa = 94 dB SPL)
（注意：音圧や電圧のデシベルは20log10であるが，パワーの場合は10log10で計算される。）
|
dB (20log) |
0 | 2 | 3 | 4 | 5 | 6 | 7 | 10 | 20 | 30 | 40 |
|
比率 |
1.00 | 1.26 | 1.41 | 1.59 | 1.78 | 2.00 | 2.24 | 3.16 | 10.0 | 31.6 | 100 |
デシベルの計算（Octave による計算）
Octaveでは，function 文を使って関数を定義できる。（end が定義の終わり。)
octave:1> function y=dB(x)
> y=20*log10(x);
> end
octave:2> dB(10)
ans = 20
octave:3>
ディジタル・フィルタによる信号処理の例
 |
この課題では，雑音波形を作り，その周波数スペクトルを求めてプロットする。その雑音波形をディジタル・ローパス・フィルタに通し
た後，周波数スペクトルを再度求めてプロットし，元の雑音波形の周波数スペクトルと比べてみる。 プロットの縦軸は20・log10のデシベル表示とする。 Octave にはフィルタを設計する機能は標準装備されていないので（MATLABの信号処理ツールボックスやScilabにはある），ローパス・フィルタ（LPF） を予め設計して，フィルタ係数を LPF10_15_coef.mat というファイルにしまってある。それを load して用いる。 （青い部分を クリックすると LPF10_15_coef.mat というファイルが Dock 内（下の方）のダウンロードに入る。ダウンロードをクリック するとLPF10_15_coef.matのアイコンが表示される。Finderウインドウを表示し、作業フォルダSIGNAL_PROC内にドラッグする。） （うまくいかない場合は，上の図をクリックしてみる。） プログラムを実行できない場合は，unix窓内で ls -l *.mat を実行して， -rw-r--r-- 1 ca12345 student 437 11 16 13:17 LPF10_15_coef.mat ファイルのサイズが 437 になっていることを確認する。 |
 |
雑音波形（縦軸は振幅，横軸は時間（サンプル点）） |
 |
ローパス・フィルタの特性（縦軸は減衰量dB，横軸は周波数Hz） |
 |
雑音波形の周波数スペクトル |
 |
ローパス・フィルタを通した後の雑音波形の周波数スペクトル（青）と，ローパス・フィルタの特性（緑） |
２）折返し防止フィルタ (anti-aliasing filter)
折返し防止フィルタ (anti-aliasing filter) は，先に述べたようにアナログ・ローパス・フィルタである。
AD変換の前や，DAの後に用いる。アナログ・フィルタの特性例を下図に示す。
 |
折返し防止のために用いるアナログ・ローパス・フィルタの選定にあたっては，遮断周波数と遮断特性が重要となる。遮断周波数はフィルタによる入力信号の減衰が始まる周波数，遮断特性は減衰の傾斜を表す。遮断特性は，２倍の周波数で入力信号が何dB減衰するかの指標 dB/oct で表示される。 アナログ・フィルタでは，上記のディジタル・フィルタ特性に見られるような阻止域における櫛型カーブは見られない。しかし，アナログ・フィルタを構成す る電子回路が雑音を発生し，それに信号が埋もれてしまうので，-120dB程度以下のフィルタの特性は実際には意味をもたない。電子回路が発生する雑音の 大きさは，カタログに記載されている信号対雑音比（S/N比Signal to Noise Ratio）または最大減衰量から判る。 通常，折返し防止に用いられるローパス・フィルタの遮断特性は，-48dB/oct〜-135dB/oct 程度である。図は，典型的な，遮断周波数が4000Hz，遮断特性が-96dB/octのフィルタ特性であるが，6000Hzにおいても-50dB程度の 減衰しかない。このことよ り，アナログ信号の標本化にあたっては，実際には，サンプリング定理で与えられる標本化周波数よりも十分高い周波数で標本化する必要のあることが分かる。 |
３）FFTによるフィルタ
離散フーリエ変換（および，その高速計算法であるFFT）を用いると，信号波形が種々の正弦波成分 に分解できることを先に示した。これは，離散フーリエ変換により，信号波形の時間領域の表現を，周波数領域の表現に変換したと考えることができる。逆に， 信号波形の周波数領域の表現を，信号波形の時間領域の表現に変換する方法として，離散逆フーリエ変換（および，その高速計算法であるIFFT inverse fast fourier transform）がある。FFTにより求めた信号波形の周波数成分のうち，除去したい周波数成分を零にしてIFFTしてやれば，フィルタが実現でき る。その例は，下の ５）ＥＣＧ分析 で示す。
４）自己相関関数
信号波形が周期性をもつ場合，その周期を計測したいことがある。例えば，下の ５）ＥＣＧ分析 で示す心電図における心拍の周期がある。このような信号波形の周期の計算法として，自己相関関数がある。N個の標本値からなる信号波形x(i)と，それを 時間遅れ(lag) t だけずらした x(i+t) とを掛け合せ，累積したものを平均化する操作である。
 |
自己相関関数は，N個の標本値からなる信号波形x(i)と，それを時間遅れ(lag)
t だけずらした x(i+t) とを掛け合せ，累積したものを平均化する計算である。 (式)
3-4-1自己相関関数は， x(t) と x(i+t) とが類似しているとき，その値が大きくなるという性質がある。周期的な信号波形は，周期毎に波形が類似しているので，その周期に相当する時間遅れ (lag)において自己相関関数の値が極大値を示す。なお，時間遅れがない t=0 のときには x(t) と x(i+t) は同じなので，自己相関関数の値が大きくなる。 |
５）ＥＣＧ分析
ここでは，ECG(心電図)波形を用
いて，FFTによるフィルタと自己相関関数の計算例を示す。ECG波形は，500Hz
でサンプリングしたものであり，体の動きにより，計測に不必要な波形のゆらぎがある。このような不必要な波形成分をアーティファクト(artifact)
という。また，商用電源から誘導した50Hzの雑音（日本では，静岡県の富士川，新潟県の糸魚川あたりを境にして東側は50Hz，西側は60Hz）が重畳
している。商用電源から誘導した成分をハム(hum)ということがある。
本例題の目的は，このような不要な波形成分をFFTによるフィルタで除去した後，自己相関関数を計算して，心拍の周期を求めることである。
|
プログラム・リスト { プログラム部分をセレクト（白黒反転）してコピー，mi 画面にペースト可能 } サンプリング周波数は500Hz，FFT点数は4096点。 ECG波形は ecg.mat というファイルにしまってある。それを load して用いる。（青い部分を クリックすると ecg.mat というファイルが Dock 内（下の方）のダ ウンロードに入る。ダウンロードをクリックするとecg.matのアイコンが表 示される。Finderウインドウを表示し、作業フォルダSIGNAL_PROC内にドラッ グする。） プログラムを実行できない場合は，unix窓内で ls -l *.mat を実行して， -rw-r--r-- 1 ca12345 student 52501 11 16 13:53 ecg.mat ファイルのサイズが 52501 になっていることを確認する。 （ ファイルのサイズが 20570位の時は，ecg.mat ファイルが圧縮されたままになっている。この場合は，アイコンを再度クリックすると，52501に展開されることがある。） 自己相関関数の計算には数分以上かかる。 ifft
逆FFTを計算する組込み関数。 |
||
 |
図3-5-1 元波形 横軸は時間（秒），縦軸は電圧(mV)。 約5秒周期の波形のうねりは，体の動きによるアーティファクトである。 |
||
 |
図3-5-2 元波形の周波数スペクトル 横軸は周波数(Hz)，縦軸は振幅。 |
||
 |
図3-5-3 元波形の周波数スペクトルの周波数軸を拡大したもの。 横軸は周波数(Hz)，縦軸は振幅。 50Hzに商用電源からの誘導による雑音成分が見える。 |
||
 |
図3-5-4 IFFTにより，周波数成分から再現されたECG波形。 図3-5-1と比べてみること。 |
||
 |
図3-5-5 50Hz近辺の周波数成分の大きさを棒グラフ表示したもの。横軸は周波数(Hz)。 FFTの点数が4096，サンプリング周波数が500Hzなので，50Hz成分は410番目近辺にあるはず。408〜413番目の周波数成分を示した。 |
||
 |
図3-5-6 同様に，1Hz以下の周波数成分の大きさを棒グラフ表示したもの。横軸は周波数(Hz)。 |
||
 |
図3-5-7 50Hz近辺と，1Hz以下の周波数成分を除去した周波数スペクトルの，周波数軸を拡大したもの。横軸は周波数(Hz)，縦軸は振幅。 図3-5-3と比べてみること。 |
||
 |
図3-5-8 不要な周波数成分を除去してＩＦＦＴにより求めた再現ECG波形（緑） と，元のECG波形（青）。ＩＦＦＴによる波形は，プ ロットが重ならないように0.2だけ上方にずらして表示してある。アー ティファクトが取り除かれているのが分かる。 |
||
 |
図3-5-9 ＩＦＦＴによる再現ECG波形より自己相関関数を計算した結果。 横軸は時間遅れ(lag)(秒)。 lag 0.9秒付近に波形のピークが見られる。心拍周期は約0.9秒(1.1Hz)であることが分かる。 |
演習課題
[6] 上のプログラムの自己相関関数の計算では，(式)3-4-1 に示した
x(t)*x(i+t)の累積結果を，(N-t)で割り平均化することが省略されている。この計算をプログラムに追加して，自己相関関数の計算結果をプ
ロットすること。
[7] Octaveには，ベクトル x の自己相関関数を，遅れ 0 から h
まで計算する組込み関数 autocor(x,h) が用意されている。
これを R=autocor( x, h ); のように使って自己相関関数の計算を行ない，結果をプロットすること。x はデータのベクトル，h
は遅れ(lag)の最大値(= m-1)。
６）窓関数
実際の時系列データのサンプル点数は，今回用いたFFT点数(4096)よりはるかに多い。実際の時系列データの一部を（例えば
4096点）切り出すと，切り出した両端の部分では波形が不連続になる（両端の外側の値が0になる）ので，FFTの結果得られる周波数スペクトルは，連続
した波形の周波数スペクトルとは異なったものとなる。波形の切り出しによる周波数スペクトルへの影響を低減するため，実際の分析では窓関数が用いられる。窓関数は，両端が0に近く，中央が1となる時間関数で，これを切り出した時系列
データに掛け合せてからFFT分析を行なう。
窓関数にはいろいろあるが，ハニング(Hanning)ウインドウ，ハミング(Hamming)ウインドウ，ブラックマン(Blackman)ウインドウ
などがよく用いられる。
| Hanning窓 | 0.5-0.5cos(2πn/N) | |
| Hamming窓 | 0.54-0.46cos(2πn/N) | |
| Blackman窓 | 0.423-0.498cos(2πn/N)+0.0792cos(4πn/N) |
[8] Hanning窓，Hamming窓，Blackman窓の形を，N=
4096として，octaveを使ってプロットすること。（ECG分析を行なう必要はない。）
[9] Octaveには，組込み関数 hamming(N),
hanning(N), blackman(N) が用意されているので，これらを使って[8]と
同様にプロットすること。
Octave には他にも多くの信
号処理用の組込み関数が用意されている。
7) 雑音を伴った時系列信号の観測
信号と雑音の性質が分かっているときは，雑音を伴った信号から信号だけを強調して観測できる場
合が多い。
信号の性質としては，その周波数スペクトルや周期性，再現性など，雑音の性質としては，周期性の有無や周波数スペクトル分布，定常的か非定常的（不規則
的）かなどを知る必要がある。雑音が特定の周波数成分を含む場合の例については，既に5) ECG分析で述べた。
i) 同期加算法（積算平均化，平均応答法など種々な名前で呼ばれている）
信号が何らかの反復する事象に対する応答であり，かつ雑音が信号と相関がない場合（無相関），信号を事象に同期させて加算すると信号を見やすくすること
ができる。正弦波信号に平均0，分散1の雑音を重畳した場合の波形と，同期加算した場合の波形をoctaveプログラムを書いて比べてみる。
 |
加算回数が１回の場合と，400回の場合とを比べる。N回加算した場合，信号の振幅はN倍になるのに対
し，雑音は 倍になることが知られている。その結果，信号対雑音比(S/N比)はN回加算により倍改善される。 |
 |
|
 |
信号と雑音の周波数成分が異なる場合には，移動平均法を用いる
ことができる。移動平均の処理は入力波形にローパス・フィルタをかけることと等しい。 プログラムでは，randnで得られた雑音波形の隣接する値を引き算することにより，高域(高周波数域)を強調した雑音波形を生成している。 下の図では，移動平均によりnoiseが減っているのが分かる。 移動平均法の別の利用法として，体動など非常に低い周波数成分をもつ雑音を低減させるために，もとの観測波形から，その移動平均波形を引き算することが ある。 信号の性質に応じて重み付け加算を行う方法もあるが，ここでは述べない。詳しくは， 「南茂夫 編著，科学計測のための波形データ処理，CQ出版社 1986 (現在絶版)」などを参考のこと。 |
 |
|
 |
|
が現れたら，作業用フォルダ内にドラッグする。
 |
周波数880Hz，長さ1秒の正弦波と，基本周波数が440Hzの３つの正弦波の加算波形をつくる。サンプリング周波数は16kHz
である。 saveaudioは，波形をファイル化する関数である。'raw'，16 は波形を加工しないで16ビットでしまうことを表す。 |
4) wavesurfer の実行
ア
イコンをダブル・クリックする。
File→Open より
testsound3.rawを選択する。
Raw File の特性を指定する。
Demonstration モードで起動する。
再生ボタンを長めに押すとイヤホンから音が出る。（再生音量に十分注意すること）
ある時刻における spectrum (spectrum section と呼ぶ) に３つの周波数成分が示されている。
ここで，Windowの種類やFFT点数を変えてみること。
control キーを押しながらマウス・ボタンを押すとメニューが現れるので，いろんな機能を試してみると面白い。
5) 音声分析
wavesurfer はもともと音声分析用に作られたアプリケーション・ソフトェアである。試してみよう。
File→New で Choose Configuration が現れたら Speech analysis を選択する。
パソコンにマイクが内蔵されているので，パソコン本体(丸い部分)に顔を近づけて録音ボタン(赤丸)をクリックすると録音が開始される。発声し終わったら
すぐに停止ボタン(黒四角)をクリックする。下図は「あいうえお」と速く発声した例である。
図は，上から順に，音声波形，サウンド・スペクトログラム，ピッチ・カーブである。横軸は時間。
音声波形は，マイクに入力された音圧の時間変化を示す。
サウンド・スペクトログラムは，周波数スペクトルの時間変化を示し，縦軸は周波数，絵の濃淡はそれぞれの周波数における(相対的な)振幅の大きさを示
す。声を発するとき，喉(声帯)から口までの音波の通り道（声道という）は共鳴管として働くので，共鳴のある周波数成分は強く表われ，図では黒く太く表示
されている。赤，緑，青，黄色の線は，共鳴周波数の時間変化を示す。この共鳴周波数はホルマント周波数( Formant Frequency )
と呼ばれ，周波数の低いものから順に第１ホルマント，第２ホルマント，・・・，として区別される。ホルマント周波数は，「あ, い, う, え,
お」などの母音の違いにより異なり，個人に共通した特徴を示す。
ピッチ・カーブは，声の高さ(声帯の振動周波数)の時間変化を示す。縦軸は周波数。基本周波数と呼ばれている。アクセントやイントネーションなどの特徴
を示す。
control キーを押しなが
らマウス・ボタンを押すとポップアップ・メニューが現れるので Spectrogram Controls... を選び，Image
Controls の小さい四角をマウスで動かして絵の濃淡を調節する。
 |
 |
画面上でマウスを左右に動かすと赤い縦線が移動する。ある時点のスペクトルの詳細を見るには，control
キーとマウス・ボタンを押す。ホルマント周波数における共鳴的特徴が現われている。
 |
 |
