Transcript 5月12日分

情報通信システム（3）
http://www10.plala.or.jp/katofmly/chiba-u/
2015年5月12日 火曜日
午後4時10分～5時40分
NTT-IT Corp.
加藤 洋一
千葉大学 3- 2
FFTの例
• フーリエ級数（三角関数版→指数関数版） →
フーリエ変換→標本化定理（ナイキスト定理，サ
ンプリング定理，ともいう） →ディジタルフーリエ
変換→高速ディジタルフーリエ変換（FFT）
WaveSpectra
波形表示部（縦軸は振幅，
横軸は時間）
周波数成分表示部（縦軸は
振幅，横軸は周波数）
千葉大学 3- 3
フーリエ級数まとめ（前回の講義）
周期的に繰り返す信号は、その繰り返しの周波数（基本周波数）の整数倍の
サイン波とコサイン波の全てに重みをかけたものの和に展開できる。
基本周波数=
１
０
繰り返しの周期T
T
-3T/2
t
T
-T/2
T
T/2
3T/2
１周期の平均値a0
重みb1と基本周波数のサイン波をかける
重みa1と基本周波数のコサイン波をかける
重みb2と2・基本周波数のサイン波をかける
全て加
える
重みa2と2・基本周波数のコサイン波をかける
重みb3と3・基本周波数のサイン波をかける
重みa3と3・基本周波数のコサイン波をかける
：
：
千葉大学 3- 4
フーリエ級数（前回の講義）
f (t )が周期Tで繰り返す周期関数であるとき、
a0  
t
t 
f (t )    an cos(n  2 )  bn sin(n  2 )
2 n1 
T
T 
と展開することができ る。ここで、
2 T /2
t
an   f (t ) cos(n  2 )dt,
T T / 2
T
T /2
2
t
bn   f (t ) sin(n  2 )dt,
T T / 2
T
である。
n  0の整数
n  0の整数
千葉大学 3- 5
周期的信号の複素フーリエ級数（前回の講義）
 jn2 Tt 
f (t )   cne

n 


T /2
1
cn   f (t )e
T T / 2
（周期的）信号
強度 f(t)
時刻t
周波数n/Tの
サイン波
t
 jn2
T
T
強度 c(n)
dt
複素フーリエ級数の場合のcnは、その絶
対値が振幅、偏角が初期位相を表す。
n/T
基本周波数のn
倍の周波数
基本周波数=1/T
フーリエ級数では、「元の信号が周期的信号である」という条
件がある。この条件がはずせれば、適用範囲が広がる。
フーリエ変換
千葉大学 3- 6
フーリエ変換
t
 jn2
1 T /2
1
T
cn   f (t )e
dt, で、 T0  T , f 0 
とすると、
T T / 2
T0
cn  f0
T0 / 2
T0 / 2
T0 / 2
T0 / 2
 j 2 n f 0t
f
(
t
)
e
dt となる。 cnを f0と

 j 2 n f 0t
f
(
t
)
e
dt

の積、即ち、面積と考 えて図にしてみると以 下のように
T /2
T /2
書ける。
f (t )dt
f (t )e j2 f t dt
0
0


T 0 / 2
0
T 0/ 2
 j 2 n f 0t
f
(
t
)
e
dt

T 0 / 2
T 0 / 2
c0
n f 0
c1
f0
f0
c2
3f0
2f0
f0
・・・・・・・・・
f0
・・・・・・・・・
T 0/ 2
 j 2 ( n1) f 0t
f
(
t
)
e
dt

cn
cn+1
nf0
f0
・・・・・・・・・
n f 0
(n+1)f0
f0
・・・・・・・・・
T 0 / 2
周波数
千葉大学 3- 7
フーリエ変換
T 0/2

F (n  f0 )
T 0/ 2

f (t )dt
f (t )e j2 f 0t dt
T 0 / 2
T 0 / 2
T 0/2

f (t )e j 2 n f 0t dt
T 0 / 2
c0
n f 0
c1
f0
f0
c2
3f0
2f0
f0
・・・・・・・・・
・・・・・・・・・
f0
cn
T 0/ 2
 j 2 ( n1) f 0t
f
(
t
)
e
dt

T 0 / 2
・・・・・・・・・
cn+1
n f 0
周波数
(n+1)f
0
・・・・・・・・・
nf0
f0
f0
周波数n  f 0での、 Cnを面積と見たときの高 さを F (n  f 0 )とすると、
F (n  f 0 ) 
T0 / 2
2n f 0t
f
(
t
)
e
dt

T0 / 2
ここで、 f 0 （即ち
0
T0  ）とすると、 n  f 0は、飛び飛びの値では
なく連続値となる。そ こで、 n  f 0  fとする。すると、
フーリエ変換の公式
F( f ) 

2 f t
f
(
t
)
e
dt


が得られる。
千葉大学 3- 8
T0を無限大にするということは…
フーリエ級数は f (t )が周期Tで繰り返す周期関数で
あるときに適用できた 。 T0を無限大にするという
ことは、周期が無限大 、即ち、 f (t )は「周期関数
である必要はない」と いうことに他ならない 。
０
t
T
-T/2
T
T/2
T
3T/2
5T/2
フーリエ変換では信号は周期的である必要はない
千葉大学 3- 9
フーリエ逆変換
f (t ) 

c e
n
t
T
n
cn
 F (n  f 0 )
f0
なので、
jn2
(複素フーリエ級数の式 )
（前々頁）
f (t ) 

 f F (n  f ) e
n
0
jn2
t
T
0

j 2nf0t
F
(
n

f
)
e
f0

0
1
(  f0 )
T
n
この式で f 0 ０ (即ちT 0  )の時を考える。



n

 は
, n  f 0はf , f 0は df
となるので、

f (t )   F ( f )e j 2ft df が得られる（フーリエ

逆変換）
千葉大学 3- 10
フーリエ変換対
F( f ) 
フーリエ変換


f (t )e j 2ft dt


フーリエ逆変換
f (t )   F ( f )e j 2ft df

フーリエ変換
信号の強度（振幅）
f(t)
フーリエ逆変換
t
時間
周波数成分の大きさ（複素数）
F(f)
f
周波数
千葉大学 3- 11
フーリエ級数とフーリエ変換の違い
フーリエ級数
対象は周期関数(信号)
フーリエ変換
T→∞としたので、対象は
周期関数に限らない
変換結果（フーリエ係数） 変換結果は連続で、周波
Cnは離散的（ Cnは数列で 数の関数となる。
あり、連続量ではない）
さて、フーリエ級数も、フーリエ変換もこのままでは、実際の役には立たない（これは
極論です。理論的な解析には役に立ちます）。というのは、実際に観測する信号 f(t)
は初等関数で表されるようなものではないので、積分が簡単には計算ではできない
からです（「アナログコンピューター」というのもあることはありますが。。。）。
「使える」形のフーリエ変換は、ディジタル計算が行える「離散フーリエ変換」です。離
散フーリエ変換にたどり着くまで、もう少しの辛抱、です。
千葉大学 3- 12
補足：マイナスの周波数？
sin(2 ft  ( 1))
sin(2 ft 1)
-3
-2
1
1
0.5
0.5
-1
1
2
3
t
-0.5
-3
-2
-1
1
3
t
-0.5
-1
2
-1
三角関数では、どちらも同じに見える。
1Hz
e
-1Hz
j ( 2 f t 1)
jy
t=0のポイント
1.5
e
j ( 2 ft ( 1))
1.5
1
1
0.5
-1
1
-0.5
-1
-1.5
x
指数関数では、
マイナスの周
波数は回転が
逆
0.5
-1
1
-0.5
-1
-1.5
千葉大学 3- 13
フーリエ変換の例題（１）
1.5
f(t)
ステップ関数：時刻-0.5から0.5までの間だけ値が１、
その他の期間は値が0
1.25
1
f (t )  1,  0.5  t  0.5
 0, otherwise
0.75
0.5
0.25
-1
-0.5
t
0.5
1


j f
 j f
1
/
2
1
e

e
F ( f )   e j 2ft dt 
e j 2ft 1/ 2 
 j 2f
j 2f
1/ 2
F(f)
sin(f )

1
f
0.75
1/ 2
0.5
0.25
-10
-5
5
10
f
時間的には、有限時間（-0.5 =< t =< 0.5）にしか存在しない（0ではない）f(t)が、無限に
続く無数の周波数の波から成り立っている、ということになる。
フーリエ級数のときの結果と比較 fourier1.gcd: Step Shape
千葉大学 3- 14
フーリエ変換の例題（2）
1.5
F(f)
F ( f )  1,  0.5  f  0.5
 0, otherwise
1.25
1
0.75
0.5
0.25
-1
-0.5
周波数が0.5より上では値が0、ということは、「帯
域制限」されている、ということである。
f
0.5
1


j t
 j t
1
/
2
1
e

e
f (t )   e j 2ft df 
e j 2ft 1/ 2 
j 2t
j 2t
1/ 2
sin(t )
f(t)

t
1/ 2
1
0.75
0.5
0.25
-10
-5
5
10
t
周波数帯域は制限されているが、信号は無限の時間存在する
（「信号が無限の時間存在する」。。なんだか「無理」がありそうな感じですね。。。）
千葉大学 3- 15
時間領域と周波数領域
信号の周波数領域での表現
信号の時間領域での表現
信号強度
1.5
f(t)
周波数成分の大きさ F(f)
1.25
1
1
0.75
0.75
0.5
0.5
時間
t
0.25
-1
-0.5
0.5
1
周波数
0.25
-10
-5
5
• フーリエ級数やフーリエ変換は、信号の時間領域で
の表現を、周波数領域での表現に変える変換
• 両者は同じ現象を異なる視点で見ているだけ（逆変
換可能な変換である）
信号強度と周波数成分の大きさは複素数
10
f
千葉大学 3- 16
ちょっと息抜き
• 人間の耳は「位相」を聞き分けることができる
か？
• 近接した２つの周波数のサイン波を同時に鳴
らすと？
WaveGenを使って実験します。
千葉大学 3- 17
標本化定理
標本化間隔
ディジタル信号
アナログ信号
標本化
時間
時間
• 標本点の間の値は保存しなくても大丈夫なのだろう
か？上記の例では、標本点を「滑らかに」つないでいけ
ば元の波形になるように見える。
時間
しかし、上記の例ではうまくいきそうにない。
時間
千葉大学 3- 18
標本化定理（サンプリング定理）の導出
F(f)
1
周波数
0.5
-3
-2
-1
1
-W
W
2
3
f
周波数制限された信号を考える。つまり、－W＜ｆ＜W以外の区
間では、 F(f)は0。本来はF(f)は複素数だが、上図は便宜上実数
のように書いてある。
さらに、便宜的に、 F(f)は、周期２Wで繰り返すとする（下記）。
F(f)
便宜的に繰り返した部分
便宜的に繰り返した部分
1
0.5
-3
-3W
-2
-1
-W
1
W
2
3
f
3W
この関数をフーリエ級数で展開する（普通、フーリエ級数は、時
間領域の信号を対象にする。周波数領域の関数であるF(f)をさ
らにフーリエ級数で展開する、というのはちょっと奇妙な感じだが、
もちろん、数学的には問題ない）。
千葉大学 3- 19
標本化定理（サンプリング定理）
f (t ) 

c e
jn2
n
t
T
フーリエ級数展開の式 。 tを f , f (t )を F ( f ), Tを 2Wとすると、
n
F( f ) 

c e
n
n
jn2
f
2W


c
e
n
 jn2
f
2W
n  nとおいた。
n
t
 jn2
1 T /2
T
cn   f (t )e
dt
T T / 2
フーリエ係数を求める 式。
同じように、 tを f , f (t )を F ( f ), Tを 2W , nを  nとすると、
f
jn2
1 W
2W
cn 
F
(
f
)
e
df

2W W
W
さて、元の信号を f (t )とすると、 f (t )   F ( f )e
W
W
f
jn2
n
f ( )   F ( f )e 2W df
2W W
j 2ft
n
df t 
のとき
2W
上の式と比較すると、 c  n 
1
n
f( )
2W 2W
千葉大学 3- 20
標本化定理（サンプリング定理）
n
1 0 1 2
)は連続信号f (t )の t     
,
,
,
,   
2W
2W 2W 2W 2W
n
の時の値を示している 。つまり、全ての nに対するf ( )がきまれば、
2W
c  nというフーリエ係数が 求まり、 c  nからフーリエ級数で F ( f )が求まる。
F ( f )が求まれば、フーリエ 逆変換で、元のf (t )が求まる。
n  0.5
n
n 1
つまり、例えば、 f (
)のように f ( )と f (
)の間の値は知る必要が
2W
2W
2W
n
1
なく、 f ( )という
秒毎の値だけわかれば、元の f (t )が完全に再現で
2W
2W
きる。
cn 
1
n
f( )
2W 2W
f(
f
t
n
2W
n  0.5
2W
n 1
2W
不要
千葉大学 3- 21
標本化定理（サンプリング定理）
最大の周波数成分が WHz に制限されている信号 fw(t) を時間
間隔 1/2W でサンプルした場合、そのサンプル系列（つまりディ
ジタル信号）から元の 連続信号 fw(t) を完全に再現できる。
例：人間が聞くことができる音の周波数は、20Hzから20KHzとい
われている。音楽CDでは、44.1KHzで音楽をサンプルし、ディジ
タルデータとして保存している。即ち、2W=44.1KHzなので、
W=22.05KHzとなる。理論的には、音楽CDは22.05KHzまでの
周波数を含む信号を完全に再生できる（もちろん実際には理論
どおりにはいかないので、大体20KHz程度までの音を再生でき
るようである）。
ところで、 1/44.1KHz＝２２．６７６マイクロ秒。即ち、音楽CDは
連続信号である音楽を２２．６７６マイクロ秒ごとにサンプルし、そ
のサンプル値のみを保存している。
千葉大学 3- 22
標本化定理（完全再現の方法）
1
n
cn 
f ( ),
2W 2W

F( f ) 
c
e
n
 jn2
f
2W
から、
n

1
n  jn2 2Wf
F( f )  
f ( )e
2W
n 2W

f (t )   F ( f )e

j 2ft


1
n  jn2 2Wf j 2ft
df   
f ( )e
e df
2W
n 2W
W
 j 2 (t  n ) f 
2W

1 
n
1 
n  e

f ( ) e
df 
f ( )



2W n 2W W
2W n 2W  j 2 (t  n ) 
2W  W


n  e j ( 2Wt n)  e j ( 2Wt n)  
n sin( (2Wt  n))
  f ( )
  f( )

2W 
j 2 (2Wt  n)  n 2W  (2Wt  n)
n
W
n
j 2 (t 
)f
2W
千葉大学 3- 23
標本化定理（完全再現の方法）
n sin( (2Wt  n))
f (t )   f ( )
2W  (2Wt  n)
n
n
この式で、サンプル値 列f ( )から、元の信号 f (t )を再生できる。
2W

Gcalcで確認
標本化定理が成立するには、元の信号が周波数Wに制限さ
れている必要がある。この条件は、一般的な音声信号、音楽
信号、映像信号では普通に達成できる。たとえば、音楽の場
合、人間の耳に聞こえる周波数の最大値は決まっているなど。
標本化定理は、アナログ信号とディジタル信号の関係を規定
する最も重要な法則である。
サンプリング定理，ナイキストの定理，とも呼ばれています
千葉大学 3- 24
標本化定理（完全再現の確認）
W=100Hz
2W=200Hz
1/ 2W=0.005sec=5m sec
f (t ) = 3sin ( 2  15 t + 1 ) +
(f=15Hz)
 2 sin ( 2  70 t + 2 ) +
(f=70Hz)
 cos( 2  80 t + 3 )
(f=80Hz)
0 <= t < 0.2 (second)
実際の数値をプログラ ムで計算し、グラフソ フトで表示してみる。
信号のグラフと、サンプル値から再生された信号のグラフは、両端を除き、
よく一致している。両端のさらに外側では信号値は０と考えられるが、そこ
で、「急に」値が変化するため、W以上の周波数成分が存在する。このため、
両端ではうまく再生できないこととなる。
（Sampling.demをgnuplotで表示）
千葉大学 3- 25
標本化定理（完全再現ができない場合）
W=100Hz
2W=200Hz
1/ 2W=0.005sec=5m sec
f (t ) = 3sin ( 2  15 t + 1 ) +
(f=15Hz)
 2 sin ( 2  70 t + 2 ) +
(f=70Hz)
 cos( 2  120 t + 3 )
(f=120Hz)
0 <= t < 0.2 (second)
実際の数値をプログラ ムで計算し、グラフソ フトで表示してみる。
W=100Hzにもかかわらず、上記信号は120Hzの成分を含んでいる。
この場合は、再生信号とオリジナル信号との誤差が大きくなることが
わかる。
千葉大学 3- 26
予習のお願い
• 今回の講義でも使いましたが、教材として
Pythonという言語で書いたプログラムを用い
ます。
• Pythonは、最初に学習するプログラミング言
語として最適といわれています。
• Pythonに関する説明をしようと考えています。
そこで、皆さんもPythonの処理系を自分の
PCにインストールし、できれば、予めいろいろ
試してください。
• PythonのWindows版のダウンロードURLは、
本講義のホームページにあります。
千葉大学 3- 27
標本化した信号の数学的取り扱いについて
 f (nT ), t  nT , nは整数
f * (t )  
otherwise
 0,
f (t )
標本化
t
t
連続信号なので積分ができる
そこで、
このままでは、積分ができない。
t 0
 ,
 (t )  
0, otherwise

 (t )dt  T
（Diracのデルタ関数という）
 は小さな実数

a 
 (t  a) f (t )dt  Tf (a)
a

f (t )   (t  nT ) f (nT )
*
n 0
として、f*(t)を標本化された信号列を表す関数とする。f(nT)は数値の列（ディジタル）で
あるが、 f*(t)は（一応）連続関数であることに注意。各サンプルにそれぞれデルタ関
数をかけたものの総和になっている。
千葉大学 3- 28
標本化した信号の表現方法について
積分してみると、
f (t )

f * (t )   (t  nT ) f (nT )
n0

 

n0
*
f
 (t)dt 
  (t  nT ) f (nT )dt


n0

  f (nT )   (t  nT )dt

  f (nT )  T
t
f(t)の積分（斜線の面積）
f (t )  f (nT )
間隔=T
n0
t
f*(t)の積分（長方形の面積の合計）
両方とも大体同じ値になる。
千葉大学 3- 29
離散フーリエ変換（DFT)
周期的（と仮定した）なディジタル信号を扱います
f (nT )
これは含まない（次の
周期の最初のサンプ
ル）
n=N-1
t
時刻
N個のサンプル
n=0
n=1
繰り返しの周期 N T
T の意味が変わっ
た。。以前は周期、い
まはサンプル間隔
NT
0
サンプリングレート 2W
サンプリング間隔 = T = 1 / 2W

f (t )   (t  nT ) f (nT )
*
アナログ的な扱いのための記法
n0
千葉大学 3- 30
離散フーリエ変換（DFT: Discrete Fourier Transform）
N 1
f (t )   (t  nT ) f (nT ) を周期N  Tで繰り返す周期信号と する。
*
n 0
N 1
一周期分は、 f (t )   (t  nT ) f (nT ),    t  NT   と表せる。
*
n 0
ここで、  を小さな正の実数とす る。 f * (t ) のフーリエ係数を求め る。
t
 jk2
1 NT  *  jk2 NTt
1 NT  N 1
NT
ck 
f (t )e
dt 
 (t  nT ) f (nT )e
dt



NT 
NT  n0
NT 
t
nT
N 1
 jk2
 jk2
1 N 1
1
NT
NT

f (nT )   (t  nT )e
dt 
f (nT )Te


NT n0
NT n0

n
 jk2
1 N 1
N
  f (nT )e
N n 0
1 N 1  jk2 Nn
ckを Fk、 f (nT )を f nと書き直すと、 Fk   f ne
( DFTの計算式)
N n 0
信号f nはN個である。 e
 j ( k  N )2
n
N
e
 jk2
n
N
なので、 kの範囲は0,1,...,N 1。
千葉大学 3- 31
離散フーリエ逆変換（IDFT: Inverse DFT）
1 N 1  jk2 Nn
Fk   fn e
( DFTの計算式)
N n 0
離散フーリエ逆変換は 、
N 1
fn   Fk e
jk2
n
N
である。これを証明す る。
k 0
N 1
N 1
1
fme


k 0 N m 0
 jk2
m
N
e
jk 2
n
N

N 1 N 1
1
fme

N k 0 m 0
1 N 1 N 1 jk2 nNm
  fm e
,
N m 0 k 0
 fn （証明終わり）
jk2
N 1
n
m
 jk 2
N
N
e
k 0
jk2
n m
N
N, n  m

0, n  m
千葉大学 3- 32
離散フーリエ逆変換（IDFT: Inverse DFT）
N 1
e
n m
jk 2
N
k 0
N 1
e
jk2
N 1
j 2
n m
N
k 0
e
N 1
N 1

0
N , n  m, これは  e  1  N だから分かるけど、

k 0
k 0
0, n  m これは？？？？
N 1
 e
j 2
k
( n  m)
N
k 0
k
N
N 1
 ( e
j 2
k
N n m
)
k 0
の部分は、複素平面での半径１の円を N個の円弧に
k 0
分割したときの、各円 弧の始点の総和である。従って、 0になるのは明らか。
1
1
-0.5
1
0.5
0.5
-1
jy
jy
jy
0.5
1
x
-1
-0.5
0.5
0.5
1
x
-1
-0.5
0.5
-0.5
-0.5
-0.5
-1
-1
-1
N=8
これらの値の総和
N=16
N=7
1
x
千葉大学 3- 33
DFTの行列表現
1 N 1  jk2 Nn
Fk   fn e
( DFTの計算式)
N n 0
 j 2
1
N
 j 2
W  1, W  e
,W e
1 N 1
Fk   fn W kn
N n 0
W 0 W 0
W0
 F0 
 0
 F1 
1
2
W
W
W

 1
 F 2   W 0 W 2
W4

 N


 
  
W 0 W N 1 W 2( N 1)
FN  1

0
1
2
2
N
,  W
N 1
e
 j 2
N 1
N
と定義すると、
 W 0  f 0 

N 1  
f
1
 W


2( N 1) 
 f2 
 W



   
( N 1) 2  
 W
  fN  1
ところで、明らかに、 W m Nn  W n なので、上記の Wの行列の要素は、
N 個の異なる値しかない 。 F 0 は、全信号列の平均値となる。
千葉大学 3- 34
IDFTの行列表現
N 1
fn   Fk e
jk2
n
N
( IDFTの計算式)
k 0
W 0  1, W 1  e
j 2
1
N
, W 2  e
j 2
2
N
,     W ( N 1)  e
j 2
N 1
N
N 1
fn   Fk W kn
k 0
W0
W0
 f 0  W 0
 f1   0
1
2
W
W
W

 
 f 2   W 0 W 2
W 4

 


    
 fN  1 W 0 W ( N 1) W 2( N 1)

W 0  F 0 


 W ( N 1)   F1 
 W 2( N 1)   F 2 



   
( N 1) 2  
 W
 FN  1
と定義すると、
千葉大学 3- 35
DFTとIDFTの実験
DFTとIDFTならディジタル計算ができる。いろいろな信号のDFTを実際に計算してみる。
| Fk | Fk  F *k
（dft.demをgnuplotで表示）
信号が実数のときは、 Fkは以下のような関係に あることが分かる。
Fkの実部  F ( N  k )の実部
Fkの虚部  F ( N  k )の虚部
即ち、
（複素共役）
Fk  F * ( N  k )
*
*
n
 j 2kn
N
1


1
F ( N  k )    fn e
   fn e
 N n 0
 N n 0 
n
 j 2k
1 N 1
N
  fn e
 Fk
(証明終わり)
N n 0
N 1
*
 j 2 ( N k )
n
N
N 1
j 2k



以上から、実数の信号に対するDFTは、0=<k<N/2だけ計算すればよいことが分かる。