# 信号処理工学

###### tags: `eeic3S`

## 計算全般

### 伝達関数

$H_a(s) = \frac{\frac{1}{sC}}{r+\frac{1}{sC}}$

### オイラーの公式

定番の公式

$$
e^{j\theta} = cos\theta + jsin\theta
$$

### Z変換

wikiに載っているので随時参照

<https://en.wikipedia.org/wiki/Z-transform>

- zとsの関係式: $z=exp(sT)$
- $z^{-1}=exp(-sT)$
- sの一次式で近似すると$z^{-1}=1-sT$
- $s=\frac{1}{T}(1-z^{-1})$...(1)
- $s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$という式を双一次変換では使う...(2)

近似的に$s$となるのは,ラプラス変換では$s$をかけることは微分に対応しているから

### フーリエ変換

$x(t)$のフーリエ変換を$X(\omega)$として,

$$
X(\omega) = \int^{\infty}_{-\infty}x(t)e^{-j\omega t}dt
$$

## 連続時間回路から離散時間回路への変換

### 伝達関数とは

> 伝達関数 (transfer function) とはシステムへの入力を出力に変換する関数のことをいう
>

### 伝達関数の求め方

|電気回路の要素|直列回路(Z:インピーダンス)|並列回路(Y:アドミタンス)|
|-|-|-|
|抵抗|$R$|$\frac{1}{R}$|
|インダクタ|$sL$|$\frac{1}{sL}$|
|コンデンサ|$\frac{1}{sC}$|$sC$|

電気回路の要素から出力と入力を求めます

s領域(s空間)における入力を$X(s)$,出力を$Y(s)$とすると,伝達関数は

$$
H(s) = \frac{Y(s)}{X(s)}
$$

### 電気回路の伝達

- 決まりごと1
- 直列回路は $V = Z_{total} \times I$
- 並列回路は $I = Y_{total} \times V$
- 決まりごと2
- 直列回路の合成インピーダンスは$Z_{total} = Z1+Z2+Z3+...$
- 並列回路の合成アドミタンスは$Y_{total} = Y1+Y2+Y3+...$
- 決まりごと3
- インピーダンスとアドミタンスの変換は$Z=\frac{1}{Y}$

### (インパルス不変法|インパルス不変法変換|標準z変換)

デジタルのインパルス応答がアナログのインパルス応答をサンプルした例になるように変換する手法

#### 離散回路の求め方

標準z変換で離散回路の伝達関数が求める方法は以下の通り

伝達関数を
$$
H_a(s) = \sum_{n} \frac{R_n}{s+\alpha_n}
$$
とすると,その離散回路の伝達関数は
$$
H_d(z) = \sum_{n} \frac{R_n}{1-Z^{-1}-e^{-\alpha_n T}}
$$

とします

覚える部分は分母の変換です.以下にも書いておきます

$$
s+\alpha_n \rightarrow 1-Z^{-1}-e^{-\alpha_n T}
$$

#### 回路構成の求め方

- 入力から来るのが掛ける方
- 出力側が割る方

具体例は下の通り(多分間違いを含んでいます)


**注意**:出力側では`+`が`-`になる部分がある(つまり`-`がかけられている)

#### ラプラス変換からz変換を求める際の欠点

- 周波数応答は繰り返し応答になり,入力信号に対して十分高い標本化周波数を採用しないと,$\omega = \pi / T$($T$は標本化周波数)を対称軸としてエイリアシングが生じることがある
- HPFやBEF(band elimination filter)で高周波数の入力に対して必ずエイリアシングが起きるので,原理的にこれらのフィルタは作成できない

### (双一次変換法|双一次変換)

アナログ信号の周波数領域を$tanh$関数で$|\omega|<\frac{\pi}{T}$に押し込めてからz変換を行うことで折り返し誤差を塞ぐ手法

### 離散回路の求め方

双一次変換で離散回路の伝達関数を求める方法は以下の通り

$$
s = \frac{2}{T}\frac{1-Z^{-1}}{1+Z^{-1}}
$$

とおいて,これを元の式$H(s)$に代入するだけです

### 両手法の長所と短所とその回避法

|手法|長所|短所|回避法|
|-|-|-|-|
|インパルス不変法|インパルス応答を変えずに変換可能|周波数領域でエイリアシングが起きることがある|エイリアシングはTを小さくすると抑えられる|
|双一次変換|簡単な設計である･周波数応答でのエイリアシングが起きない|インパルス応答が不変でない|予め目標とする周波数特性を高周波側にシフト|

#### 双一次変換について深堀り

- 周波数軸において,ナイキスト周波数以下の領域へと非線形写像を施しているため変換後にエイリアシングが起きない
- $\frac{\omega_{A}T}{2} = tan\frac{\omega_{D}T}{2}$によって写像をしているため,ナイキスト周波数付近で元々の応答から大きく歪む(インパルス応答も)

## 回路･フィルタの位相

### 単位円の問題

- $\theta_1 = Arg(z-\beta_1)$
- $\theta_2 = Arg(z-\beta_2)$
- $\phi_1 = Arg(z-\alpha_1)$
- $\phi_2 = Arg(z-\alpha_2)$
- $|H| = \frac{|z-\beta_1|-|z-\beta_2|}{|z-\alpha_1|-|z-\alpha_2|}=\frac{Q_1Q_2}{P_1P_2}$
- $\angle{H}=\theta_1+\theta_2-\varphi_1-\varphi_2$

$H(z)$が安定である条件は極がすべて単位円の内側にある(つまり絶対値が1未満であること)

### 零点･極配置のフィルタ(最小位相フィルタ)

最小位相フィルタ:
- 極･零点がすべて単位園内にある
- アナログ: 零点がすべて左反面
- ディジタル: 零点がすべて単位園内
- 入出力特性: 等しい振幅特性を持つシステムの中で遅延(群遅延)が最小

#### ディジタルフィルタ

- 位相推移がより小さいことを望まれることが多い理由: 位相推移は応答出力の時間遅延に相当し,遅延を小さくしたい場合が多いから

### 群遅延

- 定義式: $d_g(\omega)=-\frac{d\angle H(\omega)}{d\omega}$
- 物理的な意味: 群遅延の$+t_0$は入力の時間遅れ$-t_0$に対応している.すなわち,角周波数$\omega$の信号の時間遅れを表している

### ナイキスト周波数

- 定義式:$f= \frac{1}{2T}$
- 各周波数での定義式は$\omega = \frac{\pi}{T}$

### ディジタルフィルタ

- ディジタルフィルタを設計する場合,位相推移はより小さなことが望まれる理由
- **理由**: 位相推移を小さくすることで出力の時間遅延を小さくすることが可能でディジタルフィルタとして良い性能を持つから
- ディジタルフィルタが安定動作するための条件･その理由
- **安定動作するための条件**: ディジタルフィルタの伝達関数$H(z)$の極が$z$平面上で単位円内にあれば良い
- **その理由**: 単位円にあるということは$s$平面上で左反面にあることと同じ意味であり,極の実部が負になるのでインパルス応答が発散せず安定するから
- フィルタ最小位相推移フィルタであるための追加条件
- **追加条件**: $H(z)$の零点も$z$平面上で単位円内であればよい

### LPFをFIRで実現する

- $D(\omega)$の特性を代入するときは,積分の範囲を$[\frac{rad}{sec}, -\frac{rad}{sec}]$(閉区間)にすれば良い
- **FIRフィルタの特性**(LPFであることを除く): $h_k=h_{-k}$が成り立つので,$0$に対して$h_k$は対称的である.群遅延が$\omega$によらず一定になるため出力に歪みがなくなる
- **式**: $h_{-k}=\frac{1}{-k\pi}sin(-k\omega_{c})=\frac{1}{k\pi}sink\omega_{c}=h_{k}$

## 信号推定フィルタ･適応フィルタ

### LPF･HPF･FIR･

- **LPF(low-pass filter)**
- > フィルタの一種で、なんらかの信号のうち、遮断周波数より低い周波数の成分はほとんど減衰させず、遮断周波数より高い周波数の成分を逓減させるフィルタ
- **HPF(high-pass filter)**
- > フィルタの一種で、なんらかの信号のうち、遮断周波数より高い周波数の成分はほとんど減衰させず、遮断周波数より低い周波数の成分を逓減させる
- **FIR(finite impulse response)**
- ディジタルフィルタ
- > クロネッカーのデルタ入力に対するフィルタの応答特性であるインパルス応答が「有限」であるとは、有限個の標本でゼロに安定することを意味する
- **IIR(infinite impulse response)**
- 信号処理システムの属性
- > 無限長の時間においてゼロでない値を返すインパルス応答関数を持つ
- 歪まないとは: どの各周波数成分の信号も等しい時間だけ遅延するということ.線形位相システムである条件は$H(z)=\sum_{n=0}^{N-1}h(n)z^{-n}と表したとき以下の2条件を満たすこと.1は条件より自明.2は$n=2$で満たすため,$H(z)$は線形移送システムであり,群遅延が一定で歪まない
1. すべてのnについてh(n)が実数
2. ある$n=n_0$に対して対象ないし反対象に分布していると満たすこと

### 適応フィルタ

> 最適化アルゴリズムに従ってその伝達関数を自己適応させるフィルタ

#### 離散ウィーナフィルタ(式変形はwikiに丁寧に書いてあります...)

インパルス応答を$h(n)$,観測信号を$y(k)$として,フィルタに通した推定信号$z(k)$は

$$
z(t) = \sum_{n=0}^{\infty}h(n)y_{k-n}
$$

と書けます

この際に

誤差信号を$e_k$,所望信号を$d_k$として

$$
J = \mathbb{E}[e^2_k] = \mathbb{E}[(d_k - z_k)^2] \\
= \mathbb{E}[(d_k - \sum_{n=0}^\infty h(n)y_{k-n})^2]
$$

このJが最小になるようなとき,推定信号が求まります

$\frac{\partial J}{\partial h(l)} = 0$となる条件を求めればよいです

計算をすると

$$
\frac{\partial J}{\partial h(l)} = \mathbb{E}[2(d_{k}-\sum_{n=0}^{\infty}h(n)y_{k-n})(-y_{k-l})]\\
=2\sum_{n=0}^{\infty}h(n)\mathbb{E}[y_{k-n}y_{k-l}] - 2\mathbb{E}[d_{k}y_{k-l}]\\
=0
$$

つまり
$$
\sum_{n=0}^{\infty}h(n)\mathbb{E}[y_{k-n}y_{k-l}] = \mathbb{E}[d_{k}y_{k-l}]
$$

となればよいわけです

ここで,$y_k$の自己相関関数を$\varphi_{yy}(n-l)$,$d_k$と$y_k$の相互相関関数を$\varphi_{dk}(-l)$として

$$
\varphi_{yy}(n-l) = \mathbb{E}[y_{k-n}y_{k-l}] = \varphi_{yy}(l-n)
$$
$$
\varphi_{dk}(-l) = \mathbb{E}[d_{k}y_{k-l}] = \varphi_{dk}(l)
$$

となり, $l=0,1,2,...$として以下のWiener-Hopf方程式が成り立つ

$$
\sum_{n=0}^{\infty}h(n)\varphi_{yy}(l-n)=\varphi_{dy}(l)
$$

つまり,$\varphi_{yy}$と$\varphi_{dy}$が分かればW-H方程式によりインパルス応答$h(n)$が分かるということ

これを行列で表現すると以下のようになる

$$
\left(
\begin{array}{ccccc}
\varphi_{yy}(0) & \varphi_{yy}(1) & \varphi_{yy}(2) & \ldots & \varphi_yy(N-1) \\
\varphi_{yy}(1) & \varphi_{yy}(0) & \varphi_{yy}(1) & \ldots & \varphi_yy(N-2) \\
\varphi_{yy}(0) & \varphi_{yy}(1) & \varphi_{yy}(2) & \ldots & \varphi_yy(N-3) \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\varphi_{yy}(N-1) & \varphi_{yy}(N-2) & \varphi_{yy}(N-3) & \ldots & \varphi_{yy}(N-1)
\end{array}
\right)
\left(
\begin{array}{cccc}
h(0)\\
h(1)\\
h(2)\\
\ldots\\
h(N-1)
\end{array}
\right)=
\left(
\begin{array}{cccc}
\varphi_{yy}(0)\\
\varphi_{yy}(1)\\
\varphi_{yy}(2)\\
\vdots\\
\varphi_{yy}(N-1)
\end{array}
\right)
$$

#### 連続ウィーナフィルタ

連続ウィーナフィルタは以下の通りです

$$
t \geq 0で\\
\int_0^{\infty}h(\tau)\varphi_{yy}(t-\tau)d\tau = \varphi_{dy}(t)
$$

$h(\tau)$が因果性で$0 \leq t \leq \infty$と定義されているため,下限が$0$になっているが,ここで$-\infty \leq t \leq \infty$とすると,

$$
\int_{-\infty}^{\infty}h(\tau)\varphi_{yy}(t-\tau)d\tau = \varphi_{dy}(t)
$$

両辺をフーリエ変換すると,

$$
H(jw)\Phi_{yy}(jw)=\Phi_{dy}(jw)
$$

となり,すなわち

$$
H(jw)=\frac{\Phi_{dy}(jw)}{\Phi_{yy}(jw)}=\frac{\text{観測信号と所望信号の相互スペクトル密度}}{観測信号の電力スペクトル密度}
$$

### 谷とか山とか

$$
H(z)=\frac{信号成分}{信号成分+雑音成分}
$$
で表されるとき,

1. 信号成分が雑音成分よりも十分大きいとき
- $|H(z)|\approx 1$となり,急峻なピークとなる
2. 雑音成分が信号成分よりも十分大きいとき
- $|H(x)|\approx 0$となり,谷が目立たなくなる

### フィルタ

これを計算すれば良い

$$
\varphi_{yy}(\tau) = \int_{-\infty}^{\infty}y(t - \tau)y(t)dt\\
=\int_{-\infty}^{\infty}\{x(t-\tau)+v(t-\tau)\}\{x(t)+v(t)\}dt\\
=\int_{-\infty}^{\infty}x(t-\tau)x(t)dt+\int_{-\infty}^{\infty}v(t-\tau)v(t)\\
= \varphi_{xx}(\tau)+\varphi_{vv}(\tau)
$$

よって$\Phi_{yy}(j\omega)=\Phi_{xx}(j\omega)+\Phi_{vv}(j\omega)$

$$
\varphi_{yd}(\tau) = \int_{-\infty}^{\infty}y(t - \tau)d(t)dt\\
= \int_{-\infty}^{\infty}\{x(t-\tau)+v(t-\tau)\}+x(t-\delta)dt\\
= \int_{-\infty}^{\infty}x(t-\tau)x(t-\tau)dt + = \int_{-\infty}^{\infty}v(t-\tau)x(t-\Delta)dt\\
= \varphi(\tau - \Delta)
$$

よって$\Phi_{yd}(j\omega)=\Phi_{xx}(j\omega)e^{-j\omega \Delta}$

以上より
$$
H(j\omega)=\frac{\Phi_{yd}(j\omega)}{\Phi_{yy}(j\omega)}=\frac{\Phi_{xx}(j\omega)e^{-j\omega \Delta}}{\Phi_{xx}(j\omega)+\Phi_{vv}(j\omega)}
$$

#### 逆フーリエ変換した際の副作用

逆フーリエ変換して得られる$h(t)$は,一般に$t<0$において0でない値を持つが,実際に因果律の成り立つフィルタを設計すると$t<0$において,$h(t)=0$であるものしか作ることができない

すなわち,インパルス応答をW-H方程式の解$h(t)$に完全に一致することはできない

## 短時間フーリエ変換･Wavelet変換

- 一般に窓長が長いほど正確なスペクトルが求められる
- 音声信号の周期が短いと大きな信号のある時刻とない時刻でスペクトルに差が生じたことで縦縞が発生する
- 音声時刻の時間領域における周期性(周波数領域における離散性)により,横縞は発生する

### フーリエ変換

定常的な振動波形$e^{jwt}$を基本波として,その線形わとして任意の信号を表現するのがフーリエ変換

### 短時間フーリエ変換

窓関数を使って信号が十分定常である時間幅を設定し,その窓をずらしながら区分的にフーリエ解析を行うのを短時間フーリエ変換と言います

$$
F(\omega, b) = \int_{-\infty}^{\infty}W(t-b)f(t)e^{-j\omega t}dt
$$

これはフーリエ変換に$W(t-b)$をかけただけで,$b$によって窓関数をシフトできます

### 窓関数

入力となる離散信号に窓関数をかけた結果を離散フーリエ変換します

この際に要求される条件は以下のとおりです

1. 主成分(メインローブ)の幅が狭いこと
- この幅が狭いほど,主成分の周波数分解能が高くなります
2. サイドローブの振幅が小さいこと
- この値が小さいほど,省電力のスペクトルを検出する能力が高まります

しかし,この2つの特徴はトレードオフの関係にあります

#### 短時間フーリエ変換で窓長の違いによって生じる結果の差

- 窓関数と信号の関数の**時間領域の乗算**は**周波数領域の畳み込み**で表される
- **窓関数のスペクトル**において,時間幅と周波数幅を同時に小さくすることはできない
- 時間幅を小さくするとスペクトルのメインローブが横に広がり**スペクトルもれ**が生じる
- **時間分解能**と**周波数分解能**の間には一種の**不確定性**がある

|窓長|周波数分解能|時間分解能|スペクトル|
|-|-|-|-|
|長い|高い|低い|スペクトルが細い|
|短い|低い|高い|スペクトルが太い|

### Wavelet変換

#### Wavelet とは

>数学において、局在する波、つまり、有限の長さの波もしくは速やかに減衰する波
>

一般的なマザーウェーブレットの例は以下の通りです

$$
\psi_{a,b}(t)=\frac{1}{\sqrt{\alpha}}\psi\left(\frac{t-b}{a}\right)
$$

$a$, $b$を動かすことで,様々に変形された基本波を作り出し,それらの線形和で任意の波形を表現します(フーリエ変換も基本波$e^{jwt}$を周波数と重み付き重ね合わせでした)

以下はフーリエ変換との対応です

- $a$: **窓幅の伸縮**に相当
- $b$: **時間シフト**に相当

フーリエ変換では窓長が固定された状態で様々な周波数の振動を組み合わせたが,ウェーブレット変換では窓長を可変として窓帳を伸び縮みさせることで周波数を変化させている

#### Wavelet変換と短時間フーリエ変換の違い

- ウェーブレットには平均した値が$0$になるという性質がある
- 短時間フーリエ変換をウェーブレット変換にするためには平均が$0$になるように基本波を変える必要がある

## 参考文献

- 制御工学をCで実現する（積分）<https://qiita.com/hara41/items/861a484e74d95f60c4ff#_reference-30fe7954231f3b9223d4>
- ディジタルフィルタの作り方<https://qiita.com/fukuroder/items/e1cd551b7492020da992>
- 伝達関数の求め方<http://www.kairo-nyumon.com/control_calculation2.html>
- 伝達関数法<https://ja.wikipedia.org/wiki/%E4%BC%9D%E9%81%94%E9%96%A2%E6%95%B0%E6%B3%95>
- 双一次変換<https://ja.wikipedia.org/wiki/%E5%8F%8C%E4%B8%80%E6%AC%A1%E5%A4%89%E6%8F%9B>
- ローパスフィルタ<https://ja.wikipedia.org/wiki/%E3%83%AD%E3%83%BC%E3%83%91%E3%82%B9%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF>
- ハイパスフィルタ<https://ja.wikipedia.org/wiki/%E3%83%8F%E3%82%A4%E3%83%91%E3%82%B9%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF>
- 有限インパルス応答<https://ja.wikipedia.org/wiki/%E6%9C%89%E9%99%90%E3%82%A4%E3%83%B3%E3%83%91%E3%83%AB%E3%82%B9%E5%BF%9C%E7%AD%94>
- 無限インパルス応答<https://ja.wikipedia.org/wiki/%E7%84%A1%E9%99%90%E3%82%A4%E3%83%B3%E3%83%91%E3%83%AB%E3%82%B9%E5%BF%9C%E7%AD%94>
- 信号処理論第二第9回(12/22)<http://www.sp.ipc.i.u-tokyo.ac.jp/~saruwatari/SP17_09.pdf>(Wiener-Hopf方程式に関して)
- DAVにあった"signal_shikepuri.pdf.old.pdf"
- Wiener filter<https://en.wikipedia.org/wiki/Wiener_filter>
- やる夫で学ぶディジタル信号処理<http://www.ic.is.tohoku.ac.jp/~swk/lecture/yaruodsp/main.html>
- 適応フィルタ<https://ja.wikipedia.org/wiki/%E9%81%A9%E5%BF%9C%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF>
- Z-transform<https://en.wikipedia.org/wiki/Z-transform>
- 窓関数<https://ja.wikipedia.org/wiki/%E7%AA%93%E9%96%A2%E6%95%B0>
- 窓関数　　(Window Function)<http://www7b.biglobe.ne.jp/~yizawa/InfSys1/basic/chap9/index.htm>
- ウェーブレット<https://ja.wikipedia.org/wiki/%E3%82%A6%E3%82%A7%E3%83%BC%E3%83%96%E3%83%AC%E3%83%83%E3%83%88>
- Wavelet変換<http://www.spcom.ecei.tohoku.ac.jp/~aito/wavelet/slide.pdf>
- Fourier transform<https://en.wikipedia.org/wiki/Fourier_transform>