Transcript Z変換
デジタル信号処理③
2002.5.28
前回までの内容
1. サンプリング定理
2. フーリエ変換
1.
実フーリエ級数
•
フーリエ級数展開
•
•
•
•
2.
三角関数の直交性
フーリエ級数展開により得られる情報
振幅スペクトル・パワースペクトルの物理的意味
スペクトル分析とその応用例
拡張①:複素フーリエ級数
•
3.
①
実フーリエ級数から複素フーリエ級数
拡張②:フーリエ変換
•
フーリエ変換・フーリエ逆変換(非周期関数へ拡張)
4. 離散フーリエ変換 (DFT, FFT)
5. 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 ( x)
c e
n
jnx
n
1
jnx
cn
f
(
x
)
e
dx
2
f (t )e
j 2 ft
dt
f (t ) F ( f )e j 2 ft df
離散化(実際の計測データを処理・
離散化に伴い再び周期関数を仮定)
Laplace変換,Z変換の基礎)
複素数へ拡張(簡便な表現)
離散フーリエ変換・離散フーリエ逆変換
N 1
F ( n ) f ( k )e
j
2 nk
N
k 0
2nk
2nk
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
2nk
2nk
j sin
N
N
F (n)cos
FFTを使った実際のフーリエ変換⑦
FFT (DFT)の使い方
Ts
f (t )
f (n)
A/D変換
アナログ信号
サンプリング定理
に留意する
データ数
N点
デジタル信号
FFTにより求めたフーリエ係数の解釈法
1
f
周波数分解能
(データの数は2のべき乗)
Ts
FFT
N
ナイキスト周波数
f N f
(意味のある周波数の最大値)
2
F (n) f n [Hz]のフーリエ係数
F (n)
複素数であることに注意
FFTを使った実際のフーリエ変換⑧
実際のフーリエ係数の計算結果
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番目。
FFTを使った実際のフーリエ変換⑨
離散フーリエ変換を用いた簡単なデジタルフィルタ
時間領域
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)
FFTを使った実際のフーリエ変換⑩
窓関数
6Hz
5.5Hz
整数ではない
周波数の時
実際に存在しない周波数
成分が出てくる。
FFTを使った実際のフーリエ変換⑪
窓関数
両端が合っていない
ことが不要な周波数
成分が出現する原因
両端を合わせる
方法に窓関数を
利用した方法がある
不要な周波数成分
周波数は変化させない
振幅を変化させる
FFTを使った実際のフーリエ変換⑫
窓関数
不要な周波数成分
が少なくなる。
FFTを使った実際のフーリエ変換⑬
いろいろな窓関数
矩形窓
w(n) 1
ハニング窓
2n
w(n) 0.5 0.5 cos
M
この窓関数は、
何もしない。
ハミング窓
2n
w(n) 0.54 0.46 cos
M
ブラックマン窓
2n
w(n) 0.42 0.5 cos
M
4n
0.08 cos
M
FFTを使った実際のフーリエ変換⑭
具体的な手順・使用上の注意点
1. 入力データにローパスフィルタを通して帯域を処理
•
対象データの上限周波数以上の周波数成分をカット
2. A/D変換を行う
•
サンプリング定理に注意
(対象データの上限周波数の2倍以上でサンプリングする)
3. 平均値を0にする
•
直流成分を0にするため
4. 窓関数をかける
5. FFT処理を行う
6. スペクトルを求めたり、フィルタ処理を行う。
•
使っているソフトウェアの離散フーリエ変換の定義に注意s
(離散フーリエ変換の定義は統一されていない。)
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 )
が成立するシステムのこと。
• それぞれの入力信号に対する応答が分かれば、それらの入力
を重み付き加算した入力に対応する出力は、それぞれの対応
する出力の同じ重み付け加算になるシステム。
(例)入力が2倍になると、出力が2倍になるシステムは、線形システム。
入力が2倍になると、出力が4倍になるシステムは、線形システムではない。
線形システム②
時不変システム
時不変システムとは、
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
C
R2
インパルス応答は、以下。
R1 R2
1
h(t )
exp
t
CR1
CR1 R2
y(t)
0 .7 ( 0 )
時不変性
線形性
x
x(t)
x(3) (3 3)
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 iT ) x(iT )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 ()
x(t )e
j t
dt
1
j t
x(t )
X
(
)
e
d
2
L{x(t )} X ( s )
jωの部分を
σ+jωに拡張
st
x
(
t
)
e
dt
1
st
x(t )
X
(
s
)
e
ds
2j
ラプラス変換②
ラプラス変換の性質
代数演算が簡単になる
例えば、
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 ddt
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
ddt
j t
x
(
)
y
(
t
)
d
e
dt
j ( T )
x
(
)
y
(
T
)
e
dTd
j
x
(
)
e
j T
y
(
T
)
e
dTd
Y ( s ) x ( ) e s d
Y () x()e j d
Y ( s) X (s)
Y () X ()
ほとんど同じ
Z変換①
Z変換の目的
ラプラス変換:アナログ(連続)信号の時系列信号を周波数
解析する場合、s(=σ+jω)の多項式で表現することにより、
時間領域と周波数領域との相互交換における代数演算を
簡単にするため。
Z変換:同じ目的で、デジタル信号の場合はZ変換を用いる。
(z=exp(sT) インパルス応答のz変換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
2j
離散データを扱う式にする。
離散化し、サンプリング周期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.0 z 1.0 z 1.2 z 0.8z
7
Z変換⑤
Z変換をデジタル信号処理で使う利点
• ラプラス変換と同様、代数操作が簡便である。
• Z変換は、Z 1の多項式によって表されており、 Z 1
のべき乗と時間の指標nが一致している。
• 離散的な時間系列とz変換との間に非常に簡単
な写像関係が存在する。
Z 1が時間的な1単位の遅れを表すパラメータ(遅延演算
子)である点
Z変換⑥
フーリエ変換・Laplace変換・Z変換
Z変換
フーリエ変換
X ( z ) x ( n) z n
n 0
zを周波数に変換するときは、
z e jTで変換する。
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
2j
Z変換⑦
Z変換の性質
線形性
推移定理
Z[ax1 (n) bx2 (n)] aX1[ z] bX 2[ 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変換の性質(たたみ込み)
ラプラス変換・フーリエ変換と同様、たたみ込み積分
の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変換
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;
% 移動平均をとる。最初は、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変換を使ったデジタルフィルタ②
時系列データの移動平均
どんな周波数応答を持つフィルタかを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
H ( z)
(1 z 1 )
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
H ( z)
(1 z 1 )
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( jT )
H () は、周波数応答
関数と呼ばれる
1
(1 exp( jT ))
2
1
(1 cos T j sin T )
2
z exp( jT )
を代入して、
周波数応答
関数を出す。
1
(1 cos T ) 2 sin 2 T
2
T
f
cos
cos
T:サンプリング周期
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]);
% 表示座標を調整
前回までの内容の
補足・訂正箇所
FFTを使った実際のフーリエ変換⑥
Matlabプログラム例
逆フーリエ変換する
• y4=real(ifft(fft_y2));
• sound(y4,Fs);
%逆FFTを行い実数部分だけを採用す
今回のプログラムでは、
なぜ、実数だけのデータをFFTし、
逆FFTすると、虚数項が出てくる
か?
本来、実数だけからなるデータをFFTし、逆FFTすると復元されたf(x)は、すべて
実数だけになるが、以下の理由により、虚数成分が出てくることがある。
1)計算上の丸め誤差によって、大きさがゼロに近い虚数が出てくることがある。
→したがって、実数部分だけを取り出す、という操作をしておくと安全。
2)今回は、途中で、半分のフーリエ係数を捨てているため、虚数が出てきてしま
う。(実数に関しては、完全な復元が保証される。)
→この場合は、かならず、実数部分だけを取り出す必要がある。
今回の内容
1. フーリエ変換のまとめ
1. 実フーリエ級数・複素フーリエ級数・フーリエ変換・
離散フーリエ変換
2. 実際のFFTの補足説明
•
窓関数・使用の手順
2. 線形システム
1. ラプラス変換(アナログ)
2. Z変換(デジタル)
3. Z変換を使った簡単なデジタルフィルタの分析
次回の内容
1. デジタルフィルタ
•
•
フィルタの種類
FIRフィルタ