デジタル信号処理② - AIST Digital Human Research

デジタル信号処理②
2002.5.21
前回の内容
1. デジタル信号処理の概論
•
目的：周波数分析、波形合成、デジタルフィルタ、相
関の調査
2. サンプリング定理
•
AD/DA変換の際、連続信号の持つ最大周波数の２
倍以上でサンプリングを行う必要がある
3. フーリエ級数展開
•
•
•
三角関数の直交性
フーリエ級数展開
周波数分析とその応用例の紹介
前回の内容（復習）
• 周期 2 の周期関数f(x)は以下のように
フーリエ級数展開が可能
f ( x  2)  f ( x) 周期関数
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
今回の内容
1. フーリエ級数（つづき）
•
•
フーリエ級数展開により得られる情報
パワースペクトルの物理的意味（パーシバルの関係式）
2. フーリエ級数の拡張
1. 拡張①：複素フーリエ級数
•
実フーリエ級数から複素フーリエ級数
2. 拡張②：フーリエ変換
•
•
フーリエ変換・フーリエ逆変換
離散フーリエ変換 (discrete Fourier transform:DFT)
3. ＦＦＴを使った実際のフーリエ変換
1. FFT処理の留意点
2. Matlab（数値演算処理ソフト）
計測されたデータから必要な情報を取り出す（フィルタリング）
ができるようになる。
フーリエ級数⑨
フーリエ級数展開によって得られる情報
a0 
f ( x)    an cos(nx)  bn sin(nx) 
2 n 1
a0 
   an2  bn2 cos(nx   n )
2 n 1
b 
 n  t an1  n  振幅情報
位相情報
 an 
1 2

an   0 f ( x) cos(nx)dx

 b  1 2  f ( x) sin(nx)dx
 n  0
フーリエ級数展開により、
１）振幅の情報 an2  bn2
２）位相の情報  n
を取り出すことが可能
an2  bn2 「振幅スペクトル」
と呼ばれている。
an2  bn2
「パワースペクトル」
フーリエ級数⑩
f (t )  0.3 sin(2 1000 t ) 
0.4 sin(2  2000 t ) 
0.2 sin(2  3000 t ) 
noise
FFT
振幅スペクトル
振幅スペクトルの物理的意味
2kHz
0.4
0.3
1kHz
0.2
an2  bn2
振幅スペクトルを計算
周波数[Hz]
3kHz
フーリエ級数⑪
パワースペクトルの物理的意味
フーリエ級数展開により電圧ｖ[V]が
以下のように表現可能。
a0 
v(t )    an cos(nt)  bn sin(nt) 
2 n 1
1Ω
v
1 2

a0   0 v(t )dt

1 2

an  0 v(t ) cos(nt)dt, (n  1,2,3, )


1 2

 bn   0 v(t ) sin(nt)dt, (n  1,2,3, )

このとき、
消費される平均電力P[W]を計算する。
1 2
1 2 2
P
i (t )  v(t )dt 
v (t )dt


2 0
2 0
フーリエ級数⑫
パワースペクトルの物理的意味
 a0 

v (t )     an cos(nt)  bn sin(nt)v(t ) より、
 2 n1

2

2
0
a0
2
v (t )dt 
2

2
0
a0
v(t )dt 
an

bn
 a 2  v(t ) cos(nt)dt  b 2  v(t ) sin(nt)dt 

n
n
0
0


n 1
 a02  2
2 
    an  bn 
 2 n 1



パワースペクトルの
積分が電力に相当
1 2 2
1  a02  2
2 
v (t )dt     an  bn 
平均電力 P 

0
2
2  2 n 1



フーリエ級数⑬
f (t )  0.3 sin(2 1000 t ) 
0.4 sin(2  2000 t ) 
0.2 sin(2  3000 t ) 
noise
FFT
パワースペクトル
パワースペクトルの物理的意味
an2  bn2
パワースペクトルを計算
2kHz
面積が
消費電力
に比例する量
1kHz
3kHz
周波数[Hz]
フーリエ級数⑭
振幅スペクトル
パワースペクトル
振幅スペクトルとパワースペクトルの比較
周波数[Hz]
周波数[Hz]
パワースペクトルでは、振幅スペクトルよりも周波数分布
が強調される。
フーリエ級数⑮
パーシバル（Parseval）の等式
パーシバル（Parseval）の等式





a
2
2
2 
f ( x)dx     an  bn 
 2 n1

2
0


さきほど、電力とパワースペクトルの関係のところで、導いた
上の式は、パーシバルの等式と呼ばれている。
複素フーリエ級数①
• フーリエ級数を、複素関数を使って書き表すと、
式の表現や変形が簡単になることが多い。
• ラプラス変換・Ｚ変換への拡張の基礎

a
f ( x)  0   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 ( x) 
jnx
c
e
n
n  
j 

1
1 
 jnx
cn 
f
(
x
)
e
dx

2  
f ( x) 

 ( jn)  c e
n  
jnx
n
例えば、「微分」は、「jnをかける」
という簡単な操作に置き変わる。
複素フーリエ級数②
e jx  cos x  j sin x （オイラーの公式）
e jnx  e  jnx
cos nx 
（ド・モアブルの公式）
2
e jnx  e  jnx
sin nx 
2j
ド・モアブルの公式をフーリエ級数の式に代入する
a0 
f ( x)    an cos nx  bn sin nx
2 n 1
a0   1
1

jnx
    an  jbn e  an  jbn e  jnx 
2 n 1  2
2

複素フーリエ級数③
a0
an  jbn
an  jbn
c0  , cn 
, c n 
, (n  1,2,3,, )
2
2
2
とおくと、

n 1
n 1
f ( x)  c0   cn e jnx   cn e  jnx

一方、


jnx
c
e
n
と変形できる。
n  
an  jbn
1 2
cn 

f ( x)cos nx  j sin nxdx

2
2 0
1 2
1 
 jnx
 jnx

f
(
x
)
e
dx

f
(
x
)
e
dx


0


と変形できる。
2
2
複素フーリエ級数④
ーまとめー
フーリエ級数
フーリエ級数展開

複素フーリエ級数
複素フーリエ級数展開
a
f ( x)  0   an cos(nx)  bn sin(nx) 
2 n 1
1 2
an   f ( x) cos(nx)dx
bn 

0
1
2

0
f ( x) 
cn :「複素フーリエ係数」
「スペクトル」
と呼ばれる
jnx
c
e
n
n  
1 
 jnx
cn 
f
(
x
)
e
dx

2  
f ( x) sin(nx)dx
an , bn :「実フーリエ係数」

cn
を求めることを、
「関数f(x)のスペクトルを調べる」
「関数f(x)をスペクトル分解する」
（光のスペクトル分光の類似から）
フーリエ変換①
扱える関数を周期関数から、非周期関数へ拡張する
(周期Tが∞であると考える)
下の複素フーリエ級数の式を
f ( x) 

jnx
c
e
n
n  
1 
 jny
cn 
f
(
y
)
e
dy



2

1 

f ( x)     f ( y )e  jny dy e jnx


n    2
周期を変数Tを使って表現する。
f(x)を周期Tの関数f(t)と考える。
2
x
t
T
で変数をｔに変える。
2 n
2 n
j
u
j
t
 1 T2

T
T
f (t )     T f (u )e
due

n    T
2


1
f 
T
とおくと
n
f  f  n 
T
  T
 j 2 ft
 j 2 fu
2
f (t )     T f (u )e
due f

n    2

となる。
フーリエ変換②
 T2
 j 2 ft
 j 2 fu
f (t )     T f (u )e
due f

n    2


非周期関数（Ｔが∞）を考えるために、f(x)の極限をとる。
T  f  0
  T
 j 2 ft
 j 2 fu
2
f (t )  lim    T f (u )e
due f
f 0

n    2




   f (u )e  j 2 fu due j 2 ft df

 
 
周波数fの関数F(f)とおくと
フーリエ変換③
フーリエ変換とフーリエ逆変換

F ( f )   f (t )e
 j 2 ft


f (t )   F ( f )e
j 2 ft

F ()  


f (t )e
 j t
dt
df
dt
フーリエ変換
(Fourier Transform)
フーリエ逆変換
(Fourier Inverse Transform)
フーリエ変換
(Fourier Transform)
1 
j t
フーリエ逆変換
f (t ) 
F
(

)
e
d


(Fourier Inverse Transform)
2  
離散フーリエ変換①
フーリエ変換の式を離散化

F ( f )   f (t )e
 j 2 ft


f (t )   F ( f )e

j 2 ft
dt
df
フーリエ変換
(Fourier Transform)
フーリエ逆変換
(Fourier Inverse Transform)
１）実際に使う場合、無限の長さを持った時系列データを使う
ことはない
２）サンプリングされた時系列データは、連続信号ではなく、
離散化された信号
実用上は、
コンピュータで処理できる離散フーリエ変換を用いる。
離散フーリエ変換②
離散フーリエ変換と離散フーリエ逆変換
ある時系列データf(k) (k=0,…,N-1)があったとき
N 1
F ( n )   f ( k )e
j
2 nk
N
k 0
N 1

k 0
1
f (k ) 
N
1

N
離散フーリエ変換
(Discrete Fourier Transform)
  2nk 
 2nk 
f (k )cos
  j sin 

 N 
  N 
N 1
 F ( n )e
n 0
j
2 nk
N
離散フーリエ逆変換
(Discrete Fourier Inverse Transform)
  2nk 
 2nk 
F (n)cos
  j sin 


 N 
n 0
  N 
N 1
離散フーリエ変換③
離散フーリエ変換の出力の解釈
計測時間・サンプリング周波数・計測点とフーリエ係数との関係
N 1
F ( n )   f ( k )e
j
2 nk
N
k 0
N 1

k 0
  2nk 
 2nk 
f (k )cos
  j sin 

 N 
  N 
f(k) (k=0,1,2,….,N-1)が
Fs: サンプリング周波数
T: 計測時間
により得られた時系列データの時、
1
f 
T
時間変数や周波数変数を含んで
いないため、
時間、周波数の解釈は、計測条件
から行う必要がある。
とおくと、F(n)は、周波数 f  n の成分の強さを
示している。
離散フーリエ変換④
スペクトルの計算法
離散フーリエ変換の結果からどのようにスペクトルを計算するか？
N 1
F ( n )   f ( k )e
j
2 nk
N
k 0
N 1

k 0
  2nk 
 2nk 
f (k )cos
  j sin 

 N 
  N 
N 1
F ( N  n)   f ( k )e
j
より、
2 ( N n ) k
N
k 0
N 1

k 0
 
2nk 
2nk 

f (k )cos 2k 
  j sin  2k 

N 
N 

 
 F * ( n)
F ( N  n)  F * (n)
離散フーリエ変換⑤
スペクトルの計算法
F ( N  n)  F * (n) より
F ( N  1)  F * (1)
F ( N  2)  F * (2)

F(
N
N
 1)  F * (  1)
2
2
F(n)は、N/2で折り返して、大きさが同じ
であることが分かる。
F(N/2)は、サンプリング周波数の1/2の
周波数成分を示すため、F(0)～F(N/2)
までをスペクトルの計算のために使用
する。
具体的に、どのように計算すると
「振幅スペクトル」が計算できるか
導出してみる。
離散フーリエ変換⑥
スペクトルの計算
1
f (k )  Re 
N
N 1
 F ( n)e
n 0
j
2nk
N




1
 Re 
N

N
N

1
1

2

nk
2nk 
2
2
j
j
N

N
N 
F
(
0
)

F
(
)
cos

k

F
(
n
)
e

F
(
N

n
)
e




2
n 1
n 1





1
 Re 
N

N
N

1
1

2nk
2nk 
2
2
j
j
N

k
*
N
N 
F
(
0
)

F
(
)(

1
)

F
(
n
)
e

F
(
n
)
e




2
n 1
n 1




F (n)  A(n)  jB(n)
とすると、
N / 2 1
F (0) F ( N / 2)
2nk 2  B(n)
2nk 
 2  A(n)
k


(1)   
cos(
)
sin(
)
N
N
N
N
N
N 
n 1 
離散フーリエ変換⑦
スペクトルの計算
N / 21
F (0) F ( N / 2)
2nk 2  B(n)
2nk 
 2  A(n)
k
f ( x) 

(1)   
cos(
)
sin(
)
N
N
N
N
N
N 
n 1 
f(x)が実数であるときは、半分のフーリエ係数から復元可能
n=0,N/2以外では、
2
振幅スペクトル=
N
2
A ( n)  B ( n) 
F ( n)
N
2
2
ＦＦＴを使った実際のフーリエ変換①
データ数が2のべき乗であるとき、高速に離散フーリエ変換
を行うアルゴリズム「高速フーリエ変換」(Fast Fourier Transform: FFT)
が知られている。多くの数値処理ソフトウェア(Matlabなど）で提供
されている。
実際に、Mablabを使って、
１）適当な波形をサンプリングし、
２）振幅スペクトルを計算する
３）簡単なフィルタ処理を行ってみる
不要な箇所のフーリエ係数を０にして、
逆フーリエ変換する。
ＦＦＴを使った実際のフーリエ変換②
スペクトルの計算
サンプリング周波数 8192[Hz]
1[s]計測(8192点計測）
f (t )  0.3 sin(2 1000 t ) 
0.4 sin(2  2000 t ) 
0.2 sin(2  3000 t )
の信号を計測する。
FFTをかけ振幅スペクトルを計算してみる。
ＦＦＴを使った実際のフーリエ変換③
Matlabプログラム例
•
•
•
•
•
•
•
•
•
•
•
•
•
•
Ts=5;
% 5秒
Fs=8192;
% サンプリング周波数 8192[Hz]
t=linspace(0,Ts,Ts*Fs); %サンプリングする時間のデータ 0[s]～5[s]
hz1=1000;
音を聞くために5[s]分のデータ
hz2=2000;
hz3=3000;
を作成するが、FFTに使用する
y1=0.3*sin(2 *pi *hz1*t); % 1k[Hz]
のは、1[s](8192点)であることに
y2=0.4*sin(2 *pi hz2*t); % 2k[Hz]
注意。
y3=0.2*sin(2 *pi *hz3*t); % 3k[Hz]
all=y1+y2+y3;
%合成波形を作る。
sound(all,Fs);
% 合成波形を鳴らす
f  1[ Hz]
f=linspace(0,Fs,8192);
%グラフを書くときに使う周波数のデータ
fft_y=fft(all,8192);
% 最初から8192点目までのデータをＦＦＴにかける
plot(f,abs(fft_y)/4096,'LineWidth',2.0);axis([0 4096 0 1.0]);
振幅スペクトル
周波数[Hz]
ＦＦＴを使った実際のフーリエ変換④
Matlabプログラム例
filterを定義する。
•
•
•
•
•
•
hz_l = floor(500/Fs*8192);
hz_h = floor(1500/Fs*8192);
filter=zeros(1,8192);
filter(hz_l:hz_h);
filter(hz_l:hz_h)=1;
plot(f,filter,'LineWidth',2.0);axis([0 4000 0 2.0]);
周波数 [Hz]
この部分だけを残す
ＦＦＴを使った実際のフーリエ変換⑤
Matlabプログラム例
不要なフーリエ係数を０にする。
•
•
fft_y2=fft_y.*filter;
plot(f,abs(fft_y2)/4096,'LineWidth',2.0);axis([0 4000 0 2.0]);
1k[Hz]の部分だけが
残る。2k[Hz],3k[Hz]の
周波数成分は、０になる。
ＦＦＴを使った実際のフーリエ変換⑥
Matlabプログラム例
逆フーリエ変換する
• fft_y2=fft_y.*filter;
• plot(f,abs(fft_y2)/4096,'LineWidth',2.0);axis([0 4000 0 2.0]);
• y4=real(ifft(fft_y2));
• sound(y4,Fs);
今回の内容
1. フーリエ級数（つづき）
•
•
フーリエ級数展開により得られる情報
パワースペクトルの物理的意味（パーシバルの関係式）
2. フーリエ級数の拡張
1. 拡張①：複素フーリエ級数
•
実フーリエ級数から複素フーリエ級数
2. 拡張②：フーリエ変換
•
•
フーリエ変換・フーリエ逆変換
離散フーリエ変換 (discrete Fourier transform:DFT)
3. ＦＦＴを使った実際のフーリエ変換
1. Matlab（数値演算処理ソフト）を使ったノイズ除去
次回の内容
1. 線形システム
2. Z変換, ラプラス変換
3. デジタルフィルタ

Download Report