デジタル信号処理③ 2002.5.28

前回までの内容 サンプリング定理 フーリエ変換 実フーリエ級数 ① 拡張①：複素フーリエ級数 拡張②：フーリエ変換 フーリエ級数展開 三角関数の直交性 フーリエ級数展開により得られる情報 振幅スペクトル・パワースペクトルの物理的意味 スペクトル分析とその応用例 拡張①：複素フーリエ級数 実フーリエ級数から複素フーリエ級数 拡張②：フーリエ変換 フーリエ変換・フーリエ逆変換（非周期関数へ拡張） 離散フーリエ変換　(DFT, FFT) DFTを使った実際のフーリエ変換 MATLAB（数値演算処理ソフト）を使ったデジタルフィルタ ① ②

今回の内容 フーリエ変換のまとめ 線形システム 実フーリエ級数・複素フーリエ級数・フーリエ変換・離散フーリエ変換 実際のFFTの補足説明 窓関数・使用の手順 線形システム ラプラス変換（アナログ） Z変換（デジタル） Z変換を使った簡単なデジタルフィルタの分析

前回までの内容（つづき） 各フーリエ級数・変換の関係のまとめ （実）フーリエ級数展開 フーリエ変換・フーリエ逆変換 非周期関数へ拡張 （より簡便な表現・変形 　Laplace変換,Z変換の基礎） 離散化（実際の計測データを処理・ 　　　　　離散化に伴い再び周期関数を仮定） 複素数へ拡張（簡便な表現） 離散フーリエ変換・離散フーリエ逆変換 複素フーリエ級数展開

ＦＦＴを使った実際のフーリエ変換⑦ FFT (DFT)の使い方 A/D変換 サンプリング定理 に留意する データ数 N点 アナログ信号 デジタル信号 FFTにより求めたフーリエ係数の解釈法 FFT 周波数分解能 （データの数は２のべき乗） ナイキスト周波数 （意味のある周波数の最大値） [Hz]のフーリエ係数 複素数であることに注意

ＦＦＴを使った実際のフーリエ変換⑧ 実際のフーリエ係数の計算結果 ＦＦＴを使った実際のフーリエ変換⑧ 実際のフーリエ係数の計算結果 計測時間が3.2[s] のとき、 -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.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 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] FFT この範囲(1-16)が スペクトルとして 意味のある範囲 f(n) F(n) 32点 32点 (注意）ちなみに、f(n)が実数の時、 1番目と(N/2+1)番目は、必ず実数になる ことが知られている。この例では、1番目と17番目。

ＦＦＴを使った実際のフーリエ変換⑨ 離散フーリエ変換を用いた簡単なデジタルフィルタ ＦＦＴを使った実際のフーリエ変換⑨ 離散フーリエ変換を用いた簡単なデジタルフィルタ 時間領域 周波数領域 時間領域 FFT IFFT 振幅（パワー） スペクトル計算 振幅（パワー） スペクトル計算 ノイズを含んだ データ ノイズがなくなった データ 厳密には、

ＦＦＴを使った実際のフーリエ変換⑩ 窓関数 ＦＦＴを使った実際のフーリエ変換⑩ 窓関数 6Hz 5.5Hz 実際に存在しない周波数 成分が出てくる。 整数ではない 周波数の時

ＦＦＴを使った実際のフーリエ変換⑪ 窓関数 ＦＦＴを使った実際のフーリエ変換⑪ 窓関数 両端を合わせる 方法に窓関数を 利用した方法がある 周波数は変化させない 振幅を変化させる 両端が合っていない ことが不要な周波数 成分が出現する原因 不要な周波数成分

ＦＦＴを使った実際のフーリエ変換⑫ 窓関数 不要な周波数成分 が少なくなる。

ＦＦＴを使った実際のフーリエ変換⑬ いろいろな窓関数 矩形窓 ハニング窓 ハミング窓 ブラックマン窓 この窓関数は、 何もしない。

ＦＦＴを使った実際のフーリエ変換⑭ 具体的な手順・使用上の注意点 ＦＦＴを使った実際のフーリエ変換⑭ 具体的な手順・使用上の注意点 入力データにローパスフィルタを通して帯域を処理 対象データの上限周波数以上の周波数成分をカット A/D変換を行う サンプリング定理に注意 （対象データの上限周波数の２倍以上でサンプリングする） 平均値を0にする 直流成分を0にするため 窓関数をかける FFT処理を行う スペクトルを求めたり、フィルタ処理を行う。 使っているソフトウェアの離散フーリエ変換の定義に注意ｓ（離散フーリエ変換の定義は統一されていない。）

Z変換と線形システム

線形システム① 線形システム L x(t) y(t) それぞれの入力信号に対する応答が分かれば、それらの入力を重み付き加算した入力に対応する出力は、それぞれの対応する出力の同じ重み付け加算になるシステム。 （例）入力が２倍になると、出力が２倍になるシステムは、線形システム。 　　　入力が２倍になると、出力が４倍になるシステムは、線形システムではない。

線形システム② 時不変システム 同じ入力に対しては、時刻によらず（今日も明日も）、同じ結果が得られるということ。 線形で、時不変の性質をもつシステムを 「線形時不変システム(linear time invariant)」 と呼ぶ。

たたみ込み① δ関数とインパルス応答 線形システムを分析する上で、もっとも重要な概念である 「たたみ込み（積分）(Convolution)」を理解する。 離散化されている信号を表現するのにδ(n)関数を用いる。 T 2T 3T

たたみ込み② δ関数とインパルス応答 入力として時刻0でインパルス(関数δ(t))を加えた時の 出力をインパルス応答と呼ぶ。 δ(t) たたみ込み② δ関数とインパルス応答 入力として時刻0でインパルス(関数δ(t))を加えた時の 出力をインパルス応答と呼ぶ。 δ(t) h(t) 線形システム インパルス インパルス応答

たたみ込み③ インパルス応答 x x(t) y(t) y インパルス応答は、以下。 1 2 3 4 重ね合わせ 時不変性 線形性 1.0 0.8 0.7 x(t) y(t) 0.5 0.3 y インパルス応答は、以下。 1 2 3 4 重ね合わせ

たたみ込み④ インパルス応答とたたみ込み 線形時不変システムにおいて、インパルス信号が、t=0で入力されたときの出力がh(t)の時、 と表現できる。 xが連続な信号であっても、 δ関数集まりとして表現可能 この数学的操作を 畳み込み積分(Convolution) と呼ぶ。

たたみ込み⑤ 離散データのたたみ込み 線形時不変システムにおいて、インパルスが、 t=0で入力されたときの出力パルス列がh(n)の時、 任意の入力パルス列x(n)に対する、出力パルス列y(n)は、 離散データでも連続データと同じようにたたみ込みで表される 連続データ 離散データ （通常は、入力されてから 　出力されるため）

たたみ込み⑥ たたみ込みの重要性 線形時不変システムにおいては、 インパルス応答h(t)さえ分かれば、 いかなる入力x(t)に対しても、 この式が意味することは、以下の点で重要 線形時不変システムにおいては、 インパルス応答h(t)さえ分かれば、 いかなる入力x(t)に対しても、 出力y(t)がどうなるか計算できる

ラプラス変換① ラプラス変換の定義 ラプラス変換の目的：アナログ（連続）信号の時系列信号を周波数解析する場合、s(=σ+jω)の多項式で表現することにより、時間領域と周波数領域との相互交換における代数演算を簡単にするため。 ラプラス変換の定義 jωの部分を σ+jωに拡張

ラプラス変換② ラプラス変換の性質 代数演算が簡単になる 例えば、 s t(時間領域) 1 h(t)が時不変線形システムのインパルス応答、 たたみ込み h(t)が時不変線形システムのインパルス応答、 x(t)が入力信号であるとき、H(s)は「伝達関数」と呼ばれる。

ラプラス変換③ たたみ込み ちなみに、フーリエ変換の場合 ラプラス変換 ほとんど同じ

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 X(z)×Y(z) -m

Z変換② フーリエ変換・Laplace変換・Z変換 離散データを扱う式にする。 離散化し、サンプリング周期T で正規化する。 へ拡張 ラプラス変換

Z変換③ Z変換の定義 Z変換とは、 で置き換え、有理多項式 に変換したもの (Tは、サンプリング周期) nを0,1,2,3からなる自然数 とすると、

Z変換④ Z変換の具体例 1.0 1.2 0.8 x(n)

Z変換⑤ Z変換をデジタル信号処理で使う利点 ラプラス変換と同様、代数操作が簡便である。 Z変換は、 の多項式によって表されており、　　のべき乗と時間の指標ｎが一致している。 離散的な時間系列とｚ変換との間に非常に簡単な写像関係が存在する。 　　　が時間的な１単位の遅れを表すパラメータ（遅延演算子）である点

Z変換⑥ フーリエ変換・Laplace変換・Z変換 ラプラス変換 zを周波数に変換するときは、 　　　　　で変換する。 sを周波数に変換するときは、 　　　　　で変換する。

Z変換⑦ Z変換の性質 線形性 推移定理

Z変換⑧ Z変換の性質（たたみ込み） ラプラス変換・フーリエ変換と同様、たたみ込み積分 のＺ変換は、掛け算になる。 ラプラス変換

Z変換を使ったデジタルフィルタ① 時系列データの移動平均 時系列信号の高周波成分をもっとも簡単に取り除く方法として、移動平均がある。 2点の平均値をとる場合

MATLABによるプログラム例 Fs=10000; % 10[kHz]のサンプリング周期 t=linspace(0,1,Fs);　　　　　　　% 10[kHz]で1[s]の時間データを作成 y1=0.7*sin(2*pi*100*t); % 100[Hz]の信号 y2=0.2*sin(2*pi*4000*t);　　　% 4[kHz]の信号 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]);　　　　　　　　　　　 % 表示する座標軸を調整する。

Z変換を使ったデジタルフィルタ② 時系列データの移動平均 どんな周波数応答を持つフィルタかをＺ変換を利用して 調べてみる。 Y(z) X(z) -1 Z変換 + Z 1/2 + ブロック線図

Z変換を使ったデジタルフィルタ② 時系列データの移動平均 たたみ込みの観点から見る Z変換 となり、一致する。

Z変換を使ったデジタルフィルタ③ 時系列データの移動平均 は、周波数応答 関数と呼ばれる を代入して、 周波数応答 関数を出す。 T:サンプリング周期 今,サンプリング周波数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));　　%逆FFTを行い実数部分だけを採用す sound(y4,Fs); 今回のプログラムでは、 なぜ、実数だけのデータをFFTし、 逆FFTすると、虚数項が出てくるか？ 本来、実数だけからなるデータをFFTし、逆FFTすると復元されたf(x)は、すべて実数だけになるが、以下の理由により、虚数成分が出てくることがある。 １）計算上の丸め誤差によって、大きさがゼロに近い虚数が出てくることがある。 　　→したがって、実数部分だけを取り出す、という操作をしておくと安全。 ２）今回は、途中で、半分のフーリエ係数を捨てているため、虚数が出てきてしまう。（実数に関しては、完全な復元が保証される。） 　　→この場合は、かならず、実数部分だけを取り出す必要がある。

今回の内容 フーリエ変換のまとめ 線形システム 実フーリエ級数・複素フーリエ級数・フーリエ変換・離散フーリエ変換 実際のFFTの補足説明 窓関数・使用の手順 線形システム ラプラス変換（アナログ） Z変換（デジタル） Z変換を使った簡単なデジタルフィルタの分析

次回の内容 デジタルフィルタ フィルタの種類 FIRフィルタ