2021-10-24

双二次フィルタの係数計算方法

Step1 基準ＬＰＦ設計
Step2 プリワーピング（予補償）
Step3 周波数・型変換
Step4 双一次s-z変換
おまけフィルタ特性の確認

　双二次フィルタの係数計算方法を勉強したので備忘録として残しておきます。下記のサイトが大変参考になりました。他のIIRフィルタに関する記事もぜひご覧ください。

スマホで学ぶ！実践ディジタル信号処理入門

hwswsgps.hatenablog.com

まずは簡単に今回対象とする３つのフィルタの特徴をまとめます。チェビシェフには遮断領域にリプルを持たせるチェビシェフⅡ型がありますが今回は触れません。

バタワース特性：最も代表的なフィルタ特性で、他の型よりも簡単な式で表される。振幅特性が優れており（最大平坦特性）、次数を上げていくと理想フィルタ特性に近づいていく。通過域の特性が良いため音声処理ではよく使用される。

チェビシェフ特性：通過域に一定振幅のリプルを持つのが最大の特徴である。バタワースフィルタが周波数の増加に従ってなだらかに減衰していくのに対して、チェビシェフフィルタはカットオフ周波数付近で急峻な特性を持っている。

ベッセル特性：群遅延特性に優れており、通過域の位相特性が周波数に対して直線的に変化する（位相直線）。単位ステップ応答の立ち上がりおよびリンギングが最も良好である。その代わりカットオフ特性がなだらかとなる。

設計には大きく分けて4つのStepがあります。Step1：フィルタの原型となる基準LPF（低域原型フィルタ）の作成、Step2：双一次変換の際の周波数軸のゆがみに対応するのための周波数変換（予補償、プリワーピング）、Step3：目的のカットオフ周波数とフィルタ形式に変換する周波数・型変換、Step4：アナログ→デジタル空間に変換する双一次s-z変換で成り立っています。

f:id:hsy221:20211024152025p:plain

Step1 基準ＬＰＦ設計

まずはアナログの基準LPFを設計します。説明の単純化のために双二次フィルタ（２次IIRフィルタ）の場合に限った伝達関数を示します。すべての次数に一般化した式は先ほどの参考リンクに載っています。各フィルタ型の基準LPFの伝達関数は下記の通りです。

・バタワース型

f:id:hsy221:20211024222540p:plain

・チェビシェフ型

f:id:hsy221:20211024222530j:plain

A：通過域リプル＝遮断周波数におけるゲイン

f:id:hsy221:20211024222534j:plain

f:id:hsy221:20211024222536j:plain

f:id:hsy221:20211024222532j:plain

f:id:hsy221:20211024222538j:plain

・ベッセル型

f:id:hsy221:20211024222544j:plain

Step2 プリワーピング（予補償）

後段の双一次s-z変換は周波数領域で非線形変換が発生するため、アナログフィルタとデジタルフィルタの間で周波数がイコールになりません。設計しようとしているデジタルフィルタのカットオフ周波数に対して、前もってその逆変換をアナログフィルタに対して実施します。

f:id:hsy221:20211024222542j:plain

T：サンプリング周期

Ωc：カットオフ角周波数（デジタル領域）

Fs：サンプリング周波数

Fc：カットオフ周波数

ωc：カットオフ角周波数（アナログ領域）

Step3 周波数・型変換

作成した基準LPFに対して周波数とフィルタタイプの変換を行います。ここでの角周波数ωcはStep2 プリワーピングで求めた角周波数です。作成したいフィルタタイプに合わせて、Step1で求めた基本LPF伝達関数のsを置き換えます。これを実施する事で基準LPFのカットオフ周波数がプリワーピングで求めた角周波数にシフトします。

・LPF

f:id:hsy221:20211024222526j:plain

・HPF

f:id:hsy221:20211024222524j:plain

・BPF

f:id:hsy221:20211024222519j:plain

・BRF

f:id:hsy221:20211024222521j:plain

ωc1：低域側カットオフ角周波数

ωc2：高域側カットオフ角周波数

BPFとBRFの場合は次数が4次になってしまうため、双二次フィルタとして設計するには因数分解して二つの二次式の積の形で表現する必要があります。真正面から計算すると大変なので、最初の参考リンクの表2-12を参考にするといいと思います。

Step4 双一次s-z変換

最後にアナログフィルタをデジタルフィルタに変換する双一次s-z変換を行います。下記の式をStep3で求めた伝達関数に代入します。式を整理して分母の定数項が1になるように分子分母を割ってやれば（正規化）双二次フィルタの伝達関数の形になるはずです。

f:id:hsy221:20211024222546j:plain

ちなみに過去記事でも紹介していますが、ピーキングフィルタとかシェルビングフィルタとか他のフィルタの伝達関数はこちらのpdfが見やすいです。

https://wsignal.sakura.ne.jp/onsei2007/kousaku/column.pdf

おまけフィルタ特性の確認

求めた伝達関数からフィルタ特性を求めるpythonプログラムを作成しました。設計した所望のフィルタになりましたか？私は計算ミスしまくってなかなか思い通りのフィルタが出来ませんでした。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 48000

a2 = 0.1716
a1 = 0
a0 = 1
b2 = 0.2929
b1 = 0.5858
b0 = 0.2929

a = [a0, a1, a2]
b = [b0, b1, b2]
w, h = signal.freqz(b, a, 65536)
f = w / (2 * np.pi) * fs
amp = 20 * np.log10(abs(h))

plt.figure(figsize=(10, 4))
plt.xscale('log')
plt.xlabel('Freq [Hz]')
plt.ylabel('Amp [dB]')
plt.xlim(20, fs)
plt.ylim(-120, 20)
plt.grid(which='both', axis='both')

plt.plot(f, amp, color="tab:blue")

plt.show()

2021-08-05

おもちゃキーボードを魔改造してFPGAシンセを作る

1. はじめに
2. 仕様
3. 構成

1. はじめに

　更新が一年以上空いてしまった。。。やっぱり続かなかったね。ずっとFPGAでマルチエフェクターを作っていたのですが、その設計資産を使って簡単なガジェットを何か作れないかなと思ってシンセサイザの自作をしてみる事にしました。FPGA エフェクターもそのうち記事にしたい。

　簡単に作るがコンセプトのため、モジュールなどを買い集めて回路設計はしない方針で作ります。また鍵盤はさすがに自作出来ないので、出来るだけ安いキーボードを探して下記のキーボードを購入しました。 FPGAは家で持て余していたamazonで買える激安FPGAモジュールを使用しました。FPGAモジュールは過去に記事にしていますので、そちらの記事もぜひ見て欲しいです。

CASIO(カシオ) 32ミニ鍵盤電子キーボード SA-46 [ミニキーボード]

発売日: 2010/05/20
メディア: エレクトロニクス

WINGONEER ALTERA FPGA Cyslonell EP2C5T144最小システム学習ボード開発

メディア: エレクトロニクス

hwswsgps.hatenablog.com

2. 仕様

・サンプリング周波数：48kHz

・量子化bit：16bit

・モノフォニック

・ステレオ（現状は左右から同じ音を出す、将来的にロータリー的なステレオエフェクトを実装したい）

パラメータ一覧

パラメータ		説明	値
shape		波形	sin, triangle, squear, up saw, down saw
mod		周波数モジュレーション	エンベロープ100% ~ なし ~ LFO100%
FILTER	type	フィルタタイプ	BYPASS, LPF(Q=1), LPF(Q=6), LPF(Q=12), BPF(Q=1), BPF(Q=6), BPF(Q=12),
	freq	カットオフ周波数	240Hz~2.4kHz
	mod	カットオフ周波数モジュレーション	エンベロープ100% ~ なし ~ LFO100%
LFO	shape	LFO波形	sin, triangle, squear, up saw, down saw
	rate	LFO周波数	約23Hz ~ 1.38Hz
	depth	LFO振幅	0 ~ 100%
MAIN	LFO	振幅LFO モジュレーション	0 ~ 100%
	Vol	マスターボリューム	0 ~ 100%
	Key	キー	0~+12半音
	Mode	将来拡張	-
ENV	atk	アタックタイム	5ms ~ 1.36s
	decay	ディケイタイム	5ms ~ 1.36s
	sustain	サスティンレベル	0 ~ 100%
	release	リリースタイム	5ms ~ 1.36s

3. 構成

　構成は図の通り、FPGAを中心にSA-46の鍵盤、パラメータ設定用の可変抵抗器を乗せた基板、スピーカアンプモジュールを接続します。

　可変抵抗器は12bit/8chのAD変換ICであるMCP3208を使用します。SPI接続です。スピーカアンプはDACとアンプが一つになったICがないかなーと思って探していた所見つけたMAX98357モジュールを使用します。鍵盤とスピーカはSA-46のものを流用します。鍵盤は32鍵ありますが、回路を調べてみた所8in4outのキーマトリクス方式でした。

f:id:hsy221:20210404222215p:plain

　実物の写真はこちらです。メイン基板を取っ払って代わりにFPGAモジュール等を内部配線しています。可変抵抗はユニバーサル基板に実装して白いプラ板をナットで固定、MCP3208からのSPIは元々LCD用の穴だった所に通しています。

f:id:hsy221:20210805225544j:plain

外観

f:id:hsy221:20210805225602j:plain

内部配線

　配線がめちゃめちゃ汚いですが、案の定動作がたまに安定しません。最初はスピーカの近くにMAX98357モジュールを配置していたのですが、SPIとI2Sが干渉しているのか、パラメータはバタバタするし、スピーカはノイズ祭りでうまくいきませんでした。MAX98357をFPGAボードのピンに直挿しする事で解決しました、その分スピーカ配線は長くなっていますが。

　一番苦労したのがGNDです。筐体がプラスチックで個々のモジュールのGNDをとる方法がFPGAボードと線でつなぐしかなく、まともに面でGNDをとれないのでこれも動作不安定の原因のようです。

内部のFPGA設計については次回以降で。

2020-05-06

lattice MachXO2に書き込もうとしたらパスワードを求められた？

mouserで購入したlattce MXO2(LCMXO2-640HC-4TG100C)に書き込もうとしたら、こんなエラーが出た。

-------------------------------------------------------------
Starting: "pgr_program run"
INFO - Check configuration setup: Start.
INFO - JTAG Chain Verification. No Errors.
INFO - Check configuration setup: Successful.
INFO - Device1 LCMXO2-640HC: FLASH Program
"A Password Key is required. Please select a different operation in the Advanced Security Programming Mode."
ERROR - Process Operation Failed.
INFO - Elapsed time: 00 min : 01 sec
ERROR - Operation: unsuccessful.
ERROR: pgr_program failed.
ERROR - Programming failed.
-------------------------------------------------------------

Lattice diamondを色々漁っていると対象デバイス右クリックメニューのDevice Propertiesの中にAcces Modeがあり、Advanced Security Keys ProgrammingやAdvanced Security File Programmingなるモードがあった。Operationを色々変更するとウインドウ下部にはデフォルトで8桁のパスワードが黒丸で隠されて書かれている。

Access modeやOperationを色々いじって試してみたが当然パスワードが合っているはずもなく、パスワードが違いますという旨のエラーが出る。

Latticeはネット上にもほとんど情報がなく、完全にお手上げ状態。以前秋月で買ったデバイスには書き込めた。型番を見ると同じシリーズの上位互換のはずなのに書き込み方法が異なるはずがない。

Lattice Diamondで変な設定をしてしまったのか、最初からパスワード付データが書き込まれた個体をつかまされたのか（mouserでそんなことあるのか？）。

2020-03-03

FPGAセレクタ回路の記述法による違い

process文外にif elseでセレクタ回路を記述した場合と、process文内でcase when文で記述した場合に組まれる回路をRTL Viewerで見てみた記録です。

・同時処理文(process文外)

LED0 <= DIPSW0 when(sel = "00") else
        DIPSW1 when(sel = "01") else
        DIPSW2 when(sel = "10") else
        DIPSW3 when(sel = "11") else
        '0';

f:id:hsy221:20200302232035p:plain

・process文

process(sel) begin
    case sel is
    when "00"   => LED0 <= DIPSW0;
    when "01"   => LED0 <= DIPSW1;
    when "10"   => LED0 <= DIPSW2;
    when "11"   => LED0 <= DIPSW3;
    when others => LED0 <= '0';
    end case;
end process;

f:id:hsy221:20200302232126p:plain

やっぱり両方とも記述そのままの回路って感じですね。ソースの可読性もRTLの見やすさとしてもセレクタ信号が2bit以上であればprocess文で書いたほうが良さそうです。ちなみにロジック規模は両方とも2LEで同じでした。

脱初心者のために、どういう記述をしたらどういう回路が合成されるかまで意識して設計していきたいです。

2020-02-24

FPGAでIIRフィルタ（双二次フィルタ）を設計する

1. はじめに
2. 設計
3. ソースコード

1. はじめに

FPGAにディジタルフィルタを実装したい時、FIRフィルタであればQuartus IPが用意されていますが、IIRは見つけられませんでした。そこで学びなおしも含めてVHDLでIIRフィルタを設計してみました。FPGAでディジタル信号処理を実装する時は、少数の扱いやビットあふれ（オーバーフロー）など、ソフトウェアで実装する時には無い少し面倒な点があります。

今回実装するのは双二次フィルタと呼ばれる二次のIIRフィルタです。ブロック図と伝達関数は以下の通り。a0が1となるように全体を正規化します。

f:id:hsy221:20200223174606p:plain

f:id:hsy221:20200223174410p:plain

2. 設計

FPGA上に設計する時は遅延器（Dフロップフロップ）、乗算器、加算器と後述のクリップ回路を使用します。実際の設計としてのブロック図は以下のようになります。

図中の数字はbit数を表し、ドットがついているのは固定少数形式であり、「整数部bit数 . 少数部bit数」となります。

f:id:hsy221:20200224153454p:plain

①上位4bit付加

加算器によるオーバーフロー対策です。例えば4bit同士の加算を行うとき、

1001b(4bit) + 1110b(4bit) = 10111b(5bit)

のように1bit増える可能性があるため、あらかじめ1bit付加して加算を行う事で、オーバーフローによるエラーを防ぎます。

加算器は図中には3つですが、中段の加算器は3つの信号を加算しているため、実際には2つの加算器で実現されます。よって計4つの加算器に対するオーバーフロー対策をしなければなりません。また符号付き信号にパディングを行う時は単純に0パディングをすると負数だった場合に符号が変わってしまうため、最上位bitをコピーします。例えば4bit信号に対して2bit付加する場合、

0011 → 000011

1010 → 111010

とします。こうする事で、符号と数値を維持したまま、bit拡張を行う事ができます。

②下位12bit切捨て

小数bitの切り捨てです。本当は少数bitの最上位を見て四捨五入した方が正確なのですが、入出力データの16bitに対して切り捨てだろうと四捨五入だろうと大した差はないので設計簡易化のために切り捨てにしました。

・乗算器

符号拡張された20bitデータと16bit(整数部4bit、少数部12bit)フィルタ係数の乗算を行います。普通に記述すればFPGAのエンベデッド乗算器が割り当てられ、自動で符号拡張が行われます。乗算器は二つの信号のbit長の合計が出力のbit長となります。今回は固定小数点演算なので整数部は20 + 4 = 24bit、少数部は0 + 12 = 12bitが出力のbit数となります。

・加算器

24bit信号同士の加算を行います。上述の通り、入力の段階でオーバーフロー対策はしているので、ここでは普通に加算するだけです。

・FF（フリップフロップ）

1クロック分信号を遅延させます。

・CLIP回路

オーバーフロー対策や乗算器で拡張していたbitを減らします。この時、上位bitをただ無視するだけだと、bitエラーを起こすので出力のbit数で表せる最大値より大きい場合や最小値より小さい値はそれぞれ最大値、最小値に丸めます。このように一定値より大きい（小さい）値をその一定値に丸める処理をCLIPと呼びます。例として17bit符号付き整数を16bit符号付き整数にCLIP処理をかけると下図のようになります。

f:id:hsy221:20200224144625p:plain

3. ソースコード

library IEEE;
use IEEE.std_logic_1164.all;
use IEEE.std_logic_signed.all;

entity iir_filter is
    port(
    RST         : in std_logic;                         -- Low Active
    CLK         : in std_logic;                         -- Clock
    a1          : in std_logic_vector(15 downto 0);     -- Filter Factor(signed 整数:4bit 少数:12bit）
    a2          : in std_logic_vector(15 downto 0);     -- Filter Factor(signed 整数:4bit 少数:12bit）
    b0          : in std_logic_vector(15 downto 0);     -- Filter Factor(signed 整数:4bit 少数:12bit）
    b1          : in std_logic_vector(15 downto 0);     -- Filter Factor(signed 整数:4bit 少数:12bit）
    b2          : in std_logic_vector(15 downto 0);     -- Filter Factor(signed 整数:4bit 少数:12bit）
    DATA_IN     : in std_logic_vector(15 downto 0);     -- Input Data(signed 16bit)
    DATA_OUT    : out std_logic_vector(15 downto 0);    -- Output Data(signed 16bit)
    OVF_FLAG    : out std_logic                         -- Overflow Flag(H pulse)
    );
end iir_filter;

architecture RTL of iir_filter is

    ------------------------------
    -- signal                   --
    ------------------------------
    signal s_data_in_d1     : std_logic_vector(19 downto 0);
    signal s_data_in_d2     : std_logic_vector(19 downto 0);
    signal s_data_out_d1    : std_logic_vector(19 downto 0);
    signal s_data_out_d2    : std_logic_vector(19 downto 0);
    signal s_mult_b2        : std_logic_vector(35 downto 0);
    signal s_mult_a2        : std_logic_vector(35 downto 0);
    signal s_mult_b1        : std_logic_vector(35 downto 0);
    signal s_mult_a1        : std_logic_vector(35 downto 0);
    signal s_mult_b0        : std_logic_vector(35 downto 0);
    signal s_sum2           : std_logic_vector(23 downto 0);
    signal s_sum1           : std_logic_vector(23 downto 0);
    signal s_sum0           : std_logic_vector(23 downto 0);
    signal s_data_in        : std_logic_vector(19 downto 0);
    signal s_data_out       : std_logic_vector(19 downto 0);

    begin
    
    -- bit 拡張
    s_data_in <= DATA_IN(15) & DATA_IN(15) & DATA_IN(15) & DATA_IN(15) & DATA_IN;

    -- 遅延(FF)
    process(RST, CLK) begin
        if(RST = '0') then
            s_data_in_d1  <= (others => '0');
            s_data_in_d2  <= (others => '0');
            s_data_out_d1 <= (others => '0');
            s_data_out_d2 <= (others => '0');
        elsif(CLK'event and CLK = '1') then
            s_data_in_d1  <= s_data_in;
            s_data_in_d2  <= s_data_in_d1;
            s_data_out_d1 <= s_data_out;
            s_data_out_d2 <= s_data_out_d1;
        end if;
    end process;

    -- 乗算
    s_mult_b2 <= b2 * s_data_in_d2;
    s_mult_a2 <= a2 * s_data_out_d2;
    s_mult_b1 <= b1 * s_data_in_d1;
    s_mult_a1 <= a1 * s_data_out_d1;
    s_mult_b0 <= b0 * s_data_in;

    -- 加算
    s_sum2 <= s_mult_b2(35 downto 12) - s_mult_a2(35 downto 12);
    s_sum1 <= s_mult_b1(35 downto 12) - s_mult_a1(35 downto 12) + s_sum2;
    s_sum0 <= s_mult_b0(35 downto 12) + s_sum1;
    
    -- CLIP処理
    s_data_out <= X"7FFFF" when(s_sum0 > X"07FFFF") else
                  X"80000" when(s_sum0 < X"F80000") else
                  s_sum0(19 downto 0);
    DATA_OUT <= X"7FFF" when(s_data_out > X"07FFF") else
                X"8000" when(s_data_out < X"F8000") else
                s_data_out(15 downto 0);

    -- Overflow flag
    OVF_FLAG <= '1' when(    s_sum0 > X"07FFFF" or     s_sum0 < X"F80000") else
                '1' when(s_data_out > X"07FFF"  or s_data_out < X"F8000" ) else
                '0';

end RTL;