Transcript PPT版

デジタル信号処理②
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
 bn 
n  tan   振幅情報
位相情報
 an 
1 2

an   0 f ( x) cos(nx)dx

bn  1 2 f ( x) sin(nx)dx

 0
1
フーリエ級数展開により、
１）振幅の情報 an2  bn2
２）位相の情報  n
を取り出すことが可能
an2  bn2 「振幅スペクトル」
と呼ばれている。
an2  bn2 「パワースペクトル」
フーリエ級数⑩
f (t )  0.3sin(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]が
以下のように表現可能。
1Ω
v
a0 
v(t )    an cos(nt)  bn sin(nt)
2 n1
1 2

a

 0  0 v(t )dt

1 2

a

v(t ) cos(nt)dt, (n  1,2,3,)
 n

0


1 2

b

 n  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

a
2 2
0
a0 2
0
v (t )dt   v(t )dt 
2 0
n
n

2
2
 a

v
(
t
)
cos(
nt
)
dt

b
v
(
t
)
sin(
nt
)
dt

n 0
n 0


n1
a
b
 a02  2 2 
    an  bn 
 2 n1



パワースペクトルの
積分が電力に相当
1 2 2
1  a02  2 2 
v (t )dt     an  bn 
平均電力 P 

2 0
2  2 n1



フーリエ級数⑬
f (t )  0.3sin(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



a

b


n
n 

 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 
cn   f ( x)e jnxdx
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   cne jnx   cne 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
  f ( x)e dx   f ( x)e jnxdx
と変形できる。
2 0
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 
cn   f ( x)e jnxdx
2 
f ( x) sin(nx)dx
an , bn :「実フーリエ係数」

cn
を求めることを、
「関数f(x)のスペクトルを調べる」
「関数f(x)をスペクトル分解する」
（光のスペクトル分光の類似から）
フーリエ変換①
扱える関数を周期関数から、非周期関数へ拡張する
(周期Tが∞であると考える)
下の複素フーリエ級数の式を
f ( x) 

jnx
c
e
n
n
1 
cn   f ( y)e jnydy
2 

1 

f ( x)     f ( y)e jnydy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
 T2
 j 2ft
 j 2 fu
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
k 0
j
2 nk
N
離散フーリエ変換
(Discrete Fourier Transform)
  2nk 
 2nk 
  f (k )cos
  j sin

 N 
k 0
  N 
N 1
2 nk
j
1 N 1
f (k )   F (n)e N 離散フーリエ逆変換
(Discrete Fourier Inverse Transform)
N n 0
  2nk 
1 N 1
 2nk 
  F (n)cos
  j sin

N n 0
 N 
  N 
離散フーリエ変換③
離散フーリエ変換の出力の解釈
計測時間・サンプリング周波数・計測点とフーリエ係数との関係
N 1
F (n)   f (k )e
j
2 nk
N
k 0
  2nk 
 2nk 
  f (k )cos
  j sin

 N 
k 0
  N 
N 1
f(k) (k=0,1,2,….,N-1)が
Fs: サンプリング周波数
T: 計測時間
により得られた時系列データの時、
時間変数や周波数変数を含んで
いないため、
時間、周波数の解釈は、計測条件
から行う必要がある。
1
f 
とおくと、F(n)は、周波数 f  n の成分の強さを
T
示している。
離散フーリエ変換④
スペクトルの計算法
離散フーリエ変換の結果からどのようにスペクトルを計算するか？
N 1
F (n)   f (k )e
j
2 nk
N
k 0
  2nk 
 2nk 
  f (k )cos
  j sin

 N 
k 0
  N 
N 1
N 1
F ( N  n)   f (k )e
j
より、
2 ( N n ) k
N
k 0
 
2nk 
2nk 

  f (k )cos 2k 
  j sin 2k 

N 
N 

k 0
 
 F * (n)
*
N 1
F ( N  n)  F (n)
離散フーリエ変換⑤
スペクトルの計算法
F ( N  n)  F * (n) より
F ( N 1)  F * (1)
F ( N  2)  F * (2)

N
* N
F (  1)  F ( 1)
2
2
F(n)は、N/2で折り返して、大きさが同じ
であることが分かる。
F(N/2)は、サンプリング周波数の1/2の
周波数成分を示すため、F(0)～F(N/2)
までをスペクトルの計算のために使用
する。
具体的に、どのように計算すると
「振幅スペクトル」が計算できるか
導出してみる。
離散フーリエ変換⑥
スペクトルの計算
2nk
j
 1 N 1

N
f (k )  Re  F (n)e

 N n 0

N
N
 
1
1
2

nk
2nk 
2
2
j
j
1
N


N
 Re F (0)  F ( ) cosk   F (n)e
  F ( N  n)e N 
N
2
n1
n1

 

N
N
 
1
1
2nk
2nk 
2
2
j
j
1
N


k
*
N
 Re F (0)  F ( )(1)   F (n)e
  F (n)e N 
N
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
2
2
2
A (n)  B (n)  F (n)
N
N
ＦＦＴを使った実際のフーリエ変換①
データ数が2のべき乗であるとき、高速に離散フーリエ変換
を行うアルゴリズム「高速フーリエ変換」(Fast Fourier Transform: FFT)
が知られている。多くの数値処理ソフトウェア(Matlabなど）で提供
されている。
実際に、Mablabを使って、
１）適当な波形をサンプリングし、
２）振幅スペクトルを計算する
３）簡単なフィルタ処理を行ってみる
不要な箇所のフーリエ係数を０にして、
逆フーリエ変換する。
ＦＦＴを使った実際のフーリエ変換②
スペクトルの計算
サンプリング周波数 8192[Hz]
1[s]計測(8192点計測）
f (t )  0.3sin(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. デジタルフィルタ