PowerPoint プレゼンテーション - AIST Digital Human

デジタル信号処理③
2002.5.28
前回までの内容
1. サンプリング定理
2. フーリエ変換
1.
実フーリエ級数
•
フーリエ級数展開
•
•
•
•
2.
三角関数の直交性
フーリエ級数展開により得られる情報
振幅スペクトル・パワースペクトルの物理的意味
スペクトル分析とその応用例
拡張①：複素フーリエ級数
•
実フーリエ級数から複素フーリエ級数
3. 拡張②：フーリエ変換
•
4.
5.
①
フーリエ変換・フーリエ逆変換（非周期関数へ拡張）
離散フーリエ変換 (DFT, FFT)
DFTを使った実際のフーリエ変換
•
MATLAB（数値演算処理ソフト）を使ったデジタルフィルタ
②
今回の内容
1. フーリエ変換のまとめ
1. 実フーリエ級数・複素フーリエ級数・フーリエ変換・
離散フーリエ変換
2. 実際のFFTの補足説明
•
窓関数・使用の手順
2. 線形システム
1. ラプラス変換（アナログ）
2. Z変換（デジタル）
3. Z変換を使った簡単なデジタルフィルタの分析
前回までの内容（つづき）
各フーリエ級数・変換の関係のまとめ
（実）フーリエ級数展開
フーリエ変換・フーリエ逆変換
a0 
f ( x)    an cos(nx)  bn sin(nx) 
2 n 1
1 2
an   f ( x) cos(nx)dx
 0
1 2
bn   f ( x) sin(nx)dx
非周期関数へ拡張
 0
（より簡便な表現・変形

F ( f )   f (t )e

複素フーリエ級数展開
f ( x) 

c e
n  
jnx
n
1 
 jnx
cn 
f
(
x
)
e
dx



2
dt

f (t )   F ( f )e j 2 ft df

離散化（実際の計測データを処理・
離散化に伴い再び周期関数を仮定）
Laplace変換,Z変換の基礎）
複素数へ拡張（簡便な表現）
 j 2 ft
離散フーリエ変換・離散フーリエ逆変換
N 1
F ( n )   f ( k )e
j
2 nk
N
k 0
  2nk 
 2nk 
  f (k )cos
  j sin 

N
N



k 0
 
N 1
1
f (k ) 
N
1

N
N 1
 F ( n )e
j
2 nk
N
n 0
N 1

n 0

 2nk 
 2nk 
  j sin 

N 
 N 
 F (n)cos
ＦＦＴを使った実際のフーリエ変換⑦
FFT (DFT)の使い方
Ts
f (t )
f (n)
A/D変換
アナログ信号
サンプリング定理
に留意する
データ数
N点
デジタル信号
FFTにより求めたフーリエ係数の解釈法
1
f 
周波数分解能
（データの数は２のべき乗）
Ts
FFT
N
ナイキスト周波数
f N  f 
（意味のある周波数の最大値）
2
F ( n)  f  n [Hz]のフーリエ係数
F (n)
複素数であることに注意
ＦＦＴを使った実際のフーリエ変換⑧
実際のフーリエ係数の計算結果
f(n)
32点
-0.0433
0.3991
0.0010
0.0836
-0.1682
-0.1237
-0.3928
0.6362
0.1327
0.1031
-0.1061
0.0594
-0.7124
0.5028
0.3991
0.0227
0.1067
-0.0054
-0.4223
-0.3677
0.6830
-0.1205
0.1589
0.0767
-0.1692
-0.5542
0.6371
0.0834
-0.0905
0.0023
-0.0285
-0.4966
FFT
F(n)
32点
0.2865
-0.4039 - 0.2695i
-0.3387 + 0.0922i
-0.2183 - 0.3678i
0.3589 + 0.5210i
0.0980 - 4.2816i
-0.0829 + 0.0915i
-0.2634 + 0.2688i
-0.5076 - 0.1058i
-0.4996 + 0.0068i
0.7744 - 4.7456i
0.2837 + 0.9688i
0.2713 + 0.4669i
-0.1144 - 0.4914i
0.0467 + 0.5260i
-0.0816 - 3.4903i
-0.3160
-0.0816 + 3.4903i
0.0467 - 0.5260i
-0.1144 + 0.4914i
0.2713 - 0.4669i
0.2837 - 0.9688i
0.7744 + 4.7456i
-0.4996 - 0.0068i
-0.5076 + 0.1058i
-0.2634 - 0.2688i
-0.0829 - 0.0915i
0.0980 + 4.2816i
0.3589 - 0.5210i
-0.2183 + 0.3678i
-0.3387 - 0.0922i
-0.4039 + 0.2695i
0 [Hz]
0.3125[Hz]
0.625[Hz]
0.9375[Hz]
1.25[Hz]
1.5625[Hz]
1.875[Hz]
2.1875[Hz]
4.375[Hz]
4.6875[Hz]
5[Hz]
計測時間が3.2[s]
のとき、
1
f 
 0.3125[ Hz ]
3. 2
この範囲(1-16)が
スペクトルとして
意味のある範囲
(注意）ちなみに、f(n)が実数の時、
1番目と(N/2+1)番目は、必ず実数になる
ことが知られている。この例では、1番目と17番目。
ＦＦＴを使った実際のフーリエ変換⑨
離散フーリエ変換を用いた簡単なデジタルフィルタ
時間領域
f1 (n)
周波数領域
FFT
F (n)
F (n)  G(n)
振幅（パワー）
スペクトル計算
F (n)
ノイズを含んだ
データ
時間領域
IFFT
f 2 (n)
振幅（パワー）
スペクトル計算
G(n)
ノイズがなくなった
データ
厳密には、
2
F (n)
N
F (n)  G(n)
ＦＦＴを使った実際のフーリエ変換⑩
窓関数
6Hz
5.5Hz
整数ではない
周波数の時
実際に存在しない周波数
成分が出てくる。
ＦＦＴを使った実際のフーリエ変換⑪
窓関数
両端が合っていない
ことが不要な周波数
成分が出現する原因
両端を合わせる
方法に窓関数を
利用した方法がある
不要な周波数成分
周波数は変化させない
振幅を変化させる
ＦＦＴを使った実際のフーリエ変換⑫
窓関数
不要な周波数成分
が少なくなる。
ＦＦＴを使った実際のフーリエ変換⑬
いろいろな窓関数
矩形窓
w(n)  1
ハニング窓
 2n 
w(n)  0.5  0.5 cos

 M 
この窓関数は、
何もしない。
ハミング窓
 2n 
w(n)  0.54  0.46cos

 M 
ブラックマン窓
 2n 
w(n)  0.42  0.5 cos

 M 
 4n 
 0.08cos

 M 
ＦＦＴを使った実際のフーリエ変換⑭
具体的な手順・使用上の注意点
1. 入力データにローパスフィルタを通して帯域を処理
•
対象データの上限周波数以上の周波数成分をカット
2. A/D変換を行う
•
サンプリング定理に注意
（対象データの上限周波数の２倍以上でサンプリングする）
3. 平均値を0にする
•
直流成分を0にするため
4. 窓関数をかける
5. FFT処理を行う
6. スペクトルを求めたり、フィルタ処理を行う。
•
使っているソフトウェアの離散フーリエ変換の定義に注意ｓ
（離散フーリエ変換の定義は統一されていない。）
Z変換と線形システム
線形システム①
線形システム
L
x(t)
線形システムとは、
y(t)
L[ x1 (t )]  y1 (t ), L[ x2 (t )]  y2 (t )の時
L[ax1 (t )  bx2 (t )]  ay1 (t )  by2 (t )
が成立するシステムのこと。
• それぞれの入力信号に対する応答が分かれば、それらの入力
を重み付き加算した入力に対応する出力は、それぞれの対応
する出力の同じ重み付け加算になるシステム。
（例）入力が２倍になると、出力が２倍になるシステムは、線形システム。
入力が２倍になると、出力が４倍になるシステムは、線形システムではない。
線形システム②
時不変システム
時不変システムとは、
L[ x(t )]  y (t )の時
L[ x(t  T )]  y (t  T )
が成立するシステムのこと。
• 同じ入力に対しては、時刻によらず（今日も明日
も）、同じ結果が得られるということ。
• 線形で、時不変の性質をもつシステムを
「線形時不変システム(linear time invariant)」
と呼ぶ。
たたみ込み①
δ関数とインパルス応答
線形システムを分析する上で、もっとも重要な概念である
「たたみ込み（積分）(Convolution)」を理解する。
離散化されている信号を表現するのにδ(n)関数を用いる。
(n)
(n  T )
(n  2T )
(n  3T )
1, t  0
(t )  
0, t  0

 (t )dt  1

0
T
2T
3T
たたみ込み②
δ関数とインパルス応答
入力として時刻0でインパルス(関数δ(t))を加えた時の
出力をインパルス応答と呼ぶ。
δ(t)
インパルス
線形システム
h(t)
インパルス応答
たたみ込み③
インパルス応答
R1
時不変性
線形性
x
x(t)
x(3)  (3  3)
 0.7  (0)
C
R2
インパルス応答は、以下。
 R1  R2 
1
h(t ) 
exp 
t 
CR1
 CR1R2 
y(t)
1.0
0.5
0.7
0.8
0.3
x(3)  h(3  3)
 0.7  h(0)
y
01 2 3 4
重ね合わせ
y(3)  x(0)h(3)  x(1)h(2)  x(2)h(1)  x(3)h(0)
 0.5h(3  0)  1.0h(3  1)  0.3h(3  2)  0.7h(3  3)
たたみ込み④
インパルス応答とたたみ込み
線形時不変システムにおいて、
インパルス信号が、t=0で入力
されたときの出力がh(t)の時、

y(t )   h(t  ) x()d
x(t )  lim
T 0

 (t  iT ) x(iT )T
i  

 h(t )  x(t )

  (t  ) x()d
と表現できる。

xが連続な信号であっても、
δ関数集まりとして表現可能
この数学的操作を
畳み込み積分(Convolution)
と呼ぶ。
たたみ込み⑤
離散データのたたみ込み
線形時不変システムにおいて、インパルスが、
t=0で入力されたときの出力パルス列がh(n)の時、
任意の入力パルス列x(n)に対する、出力パルス列y(n)は、
離散データでも連続データと同じようにたたみ込みで表される
連続データ
y(t )   h(t  ) x()d

 h(t )  x(t )
離散データ
y ( n) 

 h( n  k ) x ( k )
k  

  h( n  k ) x ( k )
k 0
（通常は、入力されてから
出力されるため）
たたみ込み⑥
たたみ込みの重要性

y(t )   h(t  ) x()d


y ( n)   h( n  k ) x ( k )
k 0
 h(t )  x(t )
この式が意味することは、以下の点で重要
線形時不変システムにおいては、
インパルス応答h(t)さえ分かれば、
いかなる入力x(t)に対しても、
出力y(t)がどうなるか計算できる
ラプラス変換①
ラプラス変換の定義
ラプラス変換の目的：アナログ（連続）信号の時系列信号を
周波数解析する場合、s(=σ+jω)の多項式で表現することに
より、時間領域と周波数領域との相互交換における代数演
算を簡単にするため。
ラプラス変換の定義


F{x(t )}  X () 

 jt
x
(
t
)
e
dt


1
jt
x(t ) 
X
(

)
e
d

2  
jωの部分を
σ+jωに拡張
L{x(t )}  X ( s) 

 x(t )e
 st

1
st
x(t ) 
X
(
s
)
e
ds

2j 
dt
ラプラス変換②
ラプラス変換の性質
代数演算が簡単になる
例えば、
t(時間領域)
(t )
t

f ( ) d 
0
s
1
F (s)
s



 F ( s )   f (t )e  st dt 






たたみ込み

y(t ) 
 x()h(t  )d

Y ( s)  X ( s)  H ( s)
 x(t )  h(t )
h(t)が時不変線形システムのインパルス応答、
x(t)が入力信号であるとき、H(s)は「伝達関数」と呼ばれる。
ラプラス変換③
たたみ込み
ちなみに、フーリエ変換の場合
ラプラス変換
 
 
L[ x(t ) * y (t )] 

 st
x
(

)
y
(
t


)
d


e
dt

F [ x(t ) * y (t )] 
  
 


 st
 x() y(t  )e ddt


 x()  y(t  )e

 st
dtd




 x()  y(T )e



 s (T   )
dTd






 s
 sT
x
(

)
e
y
(
T
)
e
dTd



  









 jt
x
(

)
y
(
t


)
e
dtd
 
t    Tとすると、
t    Tとすると、


 jt
x
(

)
y
(
t


)
e
ddt

  
 
  


 jt
x
(

)
y
(
t


)
d


e
dt


 j ( T   )
x
(

)
y
(
T
)
e
dTd
 

 j
 jT
x
(

)
e
y
(
T
)
e
dTd





 Y ( s )  x()e  s d
 Y ()  x()e  j d
 Y (s)  X (s)
 Y ()  X ()


ほとんど同じ
Z変換①
Ｚ変換の目的
ラプラス変換：アナログ（連続）信号の時系列信号を周波数
解析する場合、s(=σ+jω)の多項式で表現することにより、
時間領域と周波数領域との相互交換における代数演算を
簡単にするため。
Z変換:同じ目的で、デジタル信号の場合はZ変換を用いる。
（z=exp(sT) インパルス応答のｚ変換H(z)は、「伝達関数」
と呼ばれる。
時間領域
x(n)
ax(n)+by(n)
x(n-m)
x(n)*y(n)
z領域(周波数領域)
X(z)
aX(z)+bY(z)
X(z)・z -m
X(z)×Y(z)
Z変換②
フーリエ変換・Laplace変換・Z変換
フーリエ変換

F{x(t )}  X () 
 jt
x
(
t
)
e
dt



1
jt
x(t ) 
X
(

)
e
d

2  
j  s    j
へ拡張
ラプラス変換

L{x(t )}  X ( s) 
 st
x
(
t
)
e
dt



1
st
x(t ) 
X
(
s
)
e
ds

2j 
離散データを扱う式にする。
離散化し、サンプリング周期T
で正規化する。

 st
x
(
t
)
e
dt




 sTn
x
(
n

T
)
e
T

n  


 
 x ( n) e
n  
sT  n


n
x
(
n
)
z

n  
Z変換③
Z変換の定義
Z変換とは、
z  e sT で置き換え、有理多項式
に変換したもの

X ( z )   x ( n) z  n

(Tは、サンプリング周期)
nを0,1,2,3からなる自然数
とすると、 
X ( z )   x ( n) z  n
n 0
Z変換④
Z変換の具体例
x(n)
1.0
1.0
1
1.2
0.8
5
6
X ( z)  1.0z  1.0z  1.2z  0.8z
7
Z変換⑤
Z変換をデジタル信号処理で使う利点
• ラプラス変換と同様、代数操作が簡便である。
• Z変換は、Z 1の多項式によって表されており、 Z 1
のべき乗と時間の指標ｎが一致している。
• 離散的な時間系列とｚ変換との間に非常に簡単
な写像関係が存在する。
Z 1が時間的な１単位の遅れを表すパラメータ（遅延演算
子）である点
Z変換⑥
フーリエ変換・Laplace変換・Z変換
Z変換
フーリエ変換

X ( z )   x ( n) z  n
n 0
zを周波数に変換するときは、
z  e jTで変換する。
sを周波数に変換するときは、
s  j で変換する。

F{x(t )}  X () 
 jt
x
(
t
)
e
dt



1
jt
x(t ) 
X
(

)
e
d

2  
ラプラス変換

L{x(t )}  X ( s) 

 st
x
(
t
)
e
dt


1
st
x(t ) 
X
(
s
)
e
ds

2j 
Z変換⑦
Z変換の性質
線形性
推移定理
Z[ax1 (n)  bx2 (n)]  aX1[ z]  bX2 [ z]

Z [ x(n  m)]   x(n  m) z n
n 0

  x ( n  m) z  ( n  m ) z  m
n 0
n  m  kとすると

  x ( n  m) z  ( n  m ) z  m  X ( z ) z  m
k 0
Z変換⑧
Z変換の性質（たたみ込み）
ラプラス変換・フーリエ変換と同様、たたみ込み積分
のＺ変換は、掛け算になる。
 
 n
Z [ x(n) * y (n)]     x(i) y (i  n)z
n   i  




i
i n 
  x(i) z   y (i  n) z 
i  
i 

 X ( z)  Y ( z)

ラプラス変換
L[ x(n) * y (n)]  X ( s )  Y ( s )
フーリエ変換
F [ x(n) * y (n)]  X ()  Y ()
Ｚ変換
Z [ x(n)  y (n)]  X ( z )  Y ( z )
Z変換を使ったデジタルフィルタ①
時系列データの移動平均
• 時系列信号の高周波成分をもっとも簡単
に取り除く方法として、移動平均がある。
2点の平均値をとる場合
1
y (n)  x(n)  x(n  1)
2
MATLABによるプログラム例
•
•
•
•
•
Fs=10000;
t=linspace(0,1,Fs);
y1=0.7*sin(2*pi*100*t);
y2=0.2*sin(2*pi*4000*t);
y=y1+y2;
•
•
•
•
y_mave(1)=0;
% 移動平均をとる。最初は、０にしておく。
for n=2:Fs,
% 移動平均をとる。
y_mave(n)=(y(n)+y(n-1))/2; % 移動平均をとる。
end
% 移動平均をとる。
•
•
•
•
•
hold off;
plot(t,y,‘LineWidth’,2.0);
%合成波形をプロットする。
hold on;
% 次のグラフと重ねて描く。
plot(t,y_mave,‘r’,‘LineWidth’,2.0); % 移動平均後のデータを描く。
axis([0 0.05 -1 1]);
% 表示する座標軸を調整する。
% 10[kHz]のサンプリング周期
% 10[kHz]で1[s]の時間データを作成
% 100[Hz]の信号
% 4[kHz]の信号
% 信号を合成する
Z変換を使ったデジタルフィルタ②
時系列データの移動平均
どんな周波数応答を持つフィルタかをＺ変換を利用して
調べてみる。
1
y (n)  x(n)  x(n  1)
2
X(z)
Z変換

1
1
Y ( z)  X ( z)  X ( z) z
2
Y ( z) 1
1
H ( z) 
 (1  z )
X ( z) 2

Z -1
+
+
Y(z)
1/2
ブロック線図
Z変換を使ったデジタルフィルタ②
時系列データの移動平均
たたみ込みの観点から見る
1
y (n)  x(n)  x(n  1)
2
Z変換

1
1
Y ( z)  X ( z)  X ( z) z
2
Y ( z) 1
1
H ( z) 
 (1  z )
X ( z) 2

1
1
h(0)  , h(1)  ,
2
2
h(n)  0 (n  1)
n
y ( n)   h( n  k ) x ( k )
k 0
1
1
 x(n  1)  x(n)
2
2
となり、一致する。
Z変換を使ったデジタルフィルタ③
時系列データの移動平均
H ()  H ( z ) | z exp( jT )
H () は、周波数応答
1
 (1  exp( jT ))
2
1
 (1  cosT  j sin T )
2
関数と呼ばれる
z exp( jT )
を代入して、
周波数応答
関数を出す。
1
(1  cosT ) 2  sin 2 T
2
T
f
T:サンプリング周期
 cos
 cos
2
fs
| H ( w) |
今,サンプリング周波数fs=1/T=10[kHz]とし、f=ω/2πを
ナイキスト周波数である5[kHz]まで変化させてグラフ
作成する。
Z変換を使ったデジタルフィルタ④
時系列データの移動平均
100Hz
フィルタの振幅
スペクトル
入力信号の振幅
スペクトル
4kHz
MATLABによるプログラム例
• hold off;
• f=linspace(0,Fs,Fs);
• H=cos(pi*f/Fs);
•
•
•
•
•
% グラフを新しく作成
% 周波数データを作成
% 周波数応答関数
fft_y=fft(y);
% 合成波をFFT
plot(f,abs(fft_y)/Fs*2,‘LineWidth’,2.0); %プロット
hold on;
% グラフを重ねて描く
plot(f,H,‘r’,‘LineWidth’,2.0); % 応答関数を描く
axis([0 5000 0 1]);
% 表示座標を調整
前回までの内容の
補足・訂正箇所
ＦＦＴを使った実際のフーリエ変換⑥
Matlabプログラム例
逆フーリエ変換する
• y4=real(ifft(fft_y2));
• sound(y4,Fs);
%逆FFTを行い実数部分だけを採用す
今回のプログラムでは、
なぜ、実数だけのデータをFFTし、
逆FFTすると、虚数項が出てくる
か？
本来、実数だけからなるデータをFFTし、逆FFTすると復元されたf(x)は、すべて
実数だけになるが、以下の理由により、虚数成分が出てくることがある。
１）計算上の丸め誤差によって、大きさがゼロに近い虚数が出てくることがある。
→したがって、実数部分だけを取り出す、という操作をしておくと安全。
２）今回は、途中で、半分のフーリエ係数を捨てているため、虚数が出てきてしま
う。（実数に関しては、完全な復元が保証される。）
→この場合は、かならず、実数部分だけを取り出す必要がある。
今回の内容
1. フーリエ変換のまとめ
1. 実フーリエ級数・複素フーリエ級数・フーリエ変換・
離散フーリエ変換
2. 実際のFFTの補足説明
•
窓関数・使用の手順
2. 線形システム
1. ラプラス変換（アナログ）
2. Z変換（デジタル）
3. Z変換を使った簡単なデジタルフィルタの分析
次回の内容
1. デジタルフィルタ
•
•
フィルタの種類
FIRフィルタ

Download Report