理想のローパスフィルタと現実のフィルタ

## はじめに

### 対象読者

デジタルフィルタについてなんとなくわかっているけどどう設計したら良いかわからない人

### 書いてる人

オーディオマニアです。
エンジニアではありません。
デジタル信号処理を学んだこもありません。
書いてる内容はすべてネット上の情報の受け売りです。

### 参考サイト

<https://jp.mathworks.com/help/signal/ug/practical-introduction-to-digital-filter-design.html>

### 環境

Linux(Ubuntu 24.04) 上の GNU Octave 8.4.0 で計算

参考サイトは MATLAB なのでサンプルプログラム(Signal Processing Toolbox)はOctaveでは動きません

CDの再生を想定して fs 88200Hz、 フィルタのカットオフ22050Hz で計算
波形を表示するときは必要に応じてアップサンプリングしてます

### お願い

素人がデジタルフィルタを解説するという無謀な記事です。
説明が雑なのは許してください。
間違っていたら指摘をお願いします。
ただし私の能力が足りないため、指摘されても記事が適切に修正できるかわかりません。

## 理想のローパスフィルタ

皆さんご存知、sinc関数です。

`sin(x)/x, sinc(0) = 1`

Octave や Python の sinc関数は正規化sinc関数です。 sin(πx)/πx

非正規化sinc関数として扱いたいならxの代わりに x/π を使います。

sinc関数の波形
![fig01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/3ad31281-9e14-9ab6-6432-b156f2b797e2.png)

前後に無限に続きます。

理想ローパスフィルタの周波数特性

振幅特性(横:周波数 Linear、縦:dB)
![fig02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/515259d8-0bb2-179d-e4f9-64ce9e86c2a9.png)


理想というだけあって以下の特徴があります。

+ 直線位相
+ 通過帯域は完全にフラット(0dB)
+ 阻止帯域はgin 0(-∞dB)
+ 遷移帯域は無い
+ 実現するためには前後に無限の長さの時間が必要

前後に無限の長さが必要なので実現不可能です。

## 現実のフィルタ

### 単純にぶった切る(矩形窓)

前後に無限に続いた信号を単純に切り詰めてみます。

波形: 青が長さ511、赤が16倍の8191です。
![fig03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/ea7d1a2f-117d-bbde-acb0-0c528737e560.png)

![fig04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/42477313-8382-4d3a-bbb7-ca111661a232.png)

青は端のところで打ち切られています。赤も端で同様に打ち切られてます。

理想ローパスフィルタを有限に打ち切ったときの特性

青が長さ511、赤が8191です。

![fig05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/24d6b9a4-3e4d-8f7e-f710-cf27cfd5e937.png)

肩部分の拡大

![fig06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/815ce32c-ba0f-24ba-0f92-c922b823e354.png)

+ 通過帯域がフラットではない(リップルがある)
+ 阻止帯域が−∞ではない
+ 遷移帯域がある

フィルタを長くすると改善しますが、長くするにも限界があります。
sinc関数はxに反比例するので長さを16倍にしても端の大きさは1/16にしかなりません。

### 窓関数

有限に打ち切るときにそのままだと端の部分でいきなりゼロになります。
この急激な変化により特性が悪化します。
このため端の部分で滑らかにゼロに近づくように窓関数をかけます。

#### 窓関数の波形  

青:Hamming、赤:Kaiser(β=8)、マゼンタ:Kaiser(β=20)

![fig07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/209252bc-404b-ad70-1ab5-49b4da42cada.png)

端の部分の波形
青:窓関数なし、赤:Kaiser(β=8)
![fig08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/804cb29a-d9d0-3dca-aa55-769a8677bf70.png)

#### 窓関数を使用したときの周波数特性

青:長さ511窓関数なし、赤:長さ8191窓関数なし、マゼンタ:長さ511窓関数Kaiser(β=8)
![fig09](/img/upsample02/fig09.png)![fig09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/dc4bfde1-50c5-a7d7-c358-fa8eef1b7385.png)

阻止帯域の減衰量は長さ8191で打ち切ったものよりも長さ1/16で窓関数を使ったもののほうが遥かに良いです。

肩部分を拡大します。
![fig10](/img/upsample02/fig10.png)![fig10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/6deb5c12-d99c-6d68-49b8-76d6ce7ea736.png)

Kaiser窓を使ったものはリプルがなくフラットです。

![fig11](/img/upsample02/fig11.png)![fig11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/c0d8a7b1-ee8f-9821-66b2-679105194d3c.png)

Kaiser窓を使用したものは肩がなだらかで遷移帯域は広くなります。

#### 窓関数を変えたときの周波数特性  

長さは511
青:Hamming、赤:Kaiser(β=8)、マゼンタ:Kaiser(β=20)
![fig12](/img/upsample02/fig12.png)![fig12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/935b3a13-99a3-c1c5-c929-8d62ed32cd1f.png)

![fig13](/img/upsample02/fig13.png)![fig13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/5189ded1-b175-6cfc-f9d6-0a1b9ba0d5c2.png)

Kaiser(β=20)が一番阻止帯域の減衰量が大きいのですが、肩はなだらかで遷移帯域が一番広くなっています。
フィルタの肩特性と阻止帯域の減衰量は相反します。
両方改善しようとするとフィルタ長を長くする必要があります。

#### 同じ窓関数で長さを変える

窓関数はKaiser(β=8)
青:長さ511、赤:長さ1023、マゼンタ:長さ8191
![fig14](/img/upsample02/fig14.png)![fig14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/22500c07-06b2-0634-5481-00698a423655.png)

![fig15](/img/upsample02/fig15.png)![fig15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/bae5e1d9-0620-ff16-d129-98292fec282c.png)


長くすると阻止帯域の減衰量も改善し、遷移帯域も狭くなります。

#### カットオフ周波数を変えてみる

長さ511、窓関数Kaiser(β=8)
青:22050Hz、赤:22000Hz、マゼンタ:21950Hz
![fig16](/img/upsample02/fig16.png)
![fig16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/a93557e2-d169-fa1e-1db9-a74faa647c0f.png)

ちゃんと50Hzずつずれています。
ちなみにカットオフ周波数でのゲインは0.5、約-6dBとなります。

## 雑に作ってみる

MATLABなんかでは必要な特性を入力すれば最小の長さのフィルタを自動で作ってくれます。
今回は手動で適当にパラメータを決めていきます。

+ サンプリング周波数 88200Hz
+ カットオフ周波数 20940Hz
+ 長さ511
+ 窓関数:Kaiser(β=20)

作ったフィルタの特性
![fig17](/img/upsample02/fig17.png)![fig17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/547684/97b3b385-4daa-7d87-6e6e-a9a0a504f41e.png)


阻止帯域は22050Hz で約-190dB、通過帯域は-0.5dBが20552Hz(22050の93%)、カットオフ(-6dB)が20940Hz(22050Hzの94.7%)と充分な特性です。

このフィルタを作るOctaveのコードを載せておきます。
Octaveのsignalパッケージにはfir1関数があってローパスフィルタを簡単に作れますが、学習目的なのでこの関数は使ってません。

~~~octave
pkg load signal; # signal パッケージをロード
fs = 88200; # サンプリング周波数
fc = 20940; # カットオフ周波数
x = -1/fs*255:1/fs:1/fs*255; # 長さ511 (255*2+1)
low = sinc(x*fc*2)*fc/fs*2; # 正規化sinc関数
low = low .* kaiser(length(low), 20)'; # 窓関数
~~~