任意の振幅と位相フィルター設計

ライブスクリプトを開く

この例では、カスタマイズされた振幅と位相の仕様をもつフィルターを設計する方法を示します。多くのフィルター設計問題では、対称性によって線形の位相応答が得られるという前提で、振幅応答のみに焦点が当てられています。しかし、場合によっては、目的とするフィルターが振幅の制約も位相の制約も満たす必要があります。

たとえば、カスタムの振幅と位相設計仕様は、データ伝送システムにおける振幅と位相の歪みのイコライズ (チャネルのイコライズ) またはオーバーサンプリングされた ADC における振幅と位相の歪みのイコライズ (望ましくないハードウェア特性の補正など) に使用できます。別の用途としては、線形位相フィルターよりも小さい群遅延と特定の次数の最小位相フィルターよりも低い歪みをもつフィルターの設計があります。

周波数応答仕様とフィルター設計

フィルター応答は通常、周波数範囲 (帯域) および各帯域で目的とするゲインによって指定されます。カスタムの振幅と位相フィルター仕様も同様ですが、位相応答も含まれており、通常は複素数値としてゲインと位相応答はともに符号化されます。ほとんどの場合、応答仕様は周波数ベクトル F = [f1, ..., fN] ("N" は増加する周波数) と、フィルターの複素応答値に対応する周波数応答ベクトル H = [h1, ..., hN] から構成されます。DSP System Toolbox™ では、fdesign.arbmagnphase を使用して、目的の周波数応答をもつフィルター仕様オブジェクトを作成することができます。仕様オブジェクトの作成が済むと、関数 design を使用して FIR フィルターまたは IIR フィルターを設計できます。FIR 設計アルゴリズムおよび IIR 設計アルゴリズムの詳細については、[1]を参照してください。

FIR 設計

この最初の例では、いくつかの FIR 設計法を比較して、複素 RF バンドパスフィルターの振幅と位相をモデル化します。まず、周波数をベクトル "F"、複素応答値をベクトル "H" として目的のフィルター仕様を読み込みます。左と右のグラフにゲインと位相周波数応答をそれぞれプロットします。

load('gainAndPhase','F','H') % Load frequency response data
plotResponse(F, H) % A helper plotting function used in this demo

$Figure contains 2 axes objects. Axes object 1 with xlabel Normalized frequency f_n, ylabel Gain response |h_n| contains an object of type scatter. Axes object 2 with xlabel Normalized frequency f_n, ylabel Phase response \angle h_n contains an object of type scatter.$

'N,F,H' の仕様パターンをもつ fdesign.arbmagnphase を使用して仕様オブジェクトを作成します。この仕様は、目的のフィルター次数 "N" と、周波数応答ベクトル "F" および "H" を受け入れます。'N,F,H' のパターンは、ナイキスト範囲全体 (すなわち、緩和 "don't care" 領域のない単一帯域仕様) に対する目的の応答を定義します。この例で、目的の応答データベクトル "F" と "H" の点数は 655 点であり、周波数領域全体で比較的密になっています。

N = 100;
f = fdesign.arbmagnphase('N,F,H', N,F,H);

関数 designmethods を使用して、この仕様オブジェクトに使用できる設計法を判定します。この場合、設計法は equiripple、firls (最小二乗法) および freqsamp (周波数サンプリング法) となります。

designmethods(f,'fir')

FIR Design Methods for class fdesign.arbmagnphase (N,F,H):


equiripple
firls
freqsamp

上のリストから適切な方法を使用し、関数 design でフィルターを設計します。'allfir' を指定して、利用可能な方法をすべて使用して設計することもできます。その場合、関数は System object の cell 配列を返します。

Hd = design(f,'allfir', SystemObject = true);

フィルターの周波数応答および公称応答を破線でプロットします。等リップル設計 Hd(1) は通過帯域で非常に良好な近似であるように見えますが、他の領域ではやや逸脱しています。最小二乗設計 Hd(2) は均一重み付け 2 次ノルムに最適化されており (どの領域にも有利ではない)、周波数サンプリングされた FIR 設計 Hd(3) は 3 通りすべての中で近似が最も良くないように見えます。

hfvt = fvtool(Hd{:});
legend(hfvt,'Equiripple Hd(1)', 'FIR Least-Squares Hd(2)','Frequency Sampling  Hd(3)', ...
    Location = 'NorthEast')
ax = hfvt.CurrentAxes; 
ax.NextPlot = 'add';
plot(ax,F,20*log10(abs(H)),'k--')

Figure Figure 1: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 5 objects of type line. These objects represent Equiripple Hd(1), FIR Least-Squares Hd(2), Frequency Sampling Hd(3).

hfvt(2) = fvtool(Hd{:}, Analysis = 'phase');
legend(hfvt(2),'Equiripple Hd(1)', 'FIR Least-Squares Hd(2)','Frequency Sampling Hd(3)')

ax = hfvt(2).CurrentAxes; 
ax.NextPlot = 'add';
plot(ax,F,unwrap(angle(H))+2*pi,'k--')

Figure Figure 2: Phase Response contains an axes object. The axes object with title Phase Response, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Phase (radians) contains 4 objects of type line. These objects represent Equiripple Hd(1), FIR Least-Squares Hd(2), Frequency Sampling Hd(3).

IIR 設計

次に、IIR フィルターを設計します。目的のフィルターは、通過帯域に線形位相をもつハーフバンドハイパスフィルターです。次の図に示すように、この仕様は周波数領域の 100 個の点で与えられます。

F = [linspace(0,.475,50) linspace(.525,1,50)];
H = [zeros(1,50) exp(-1j*pi*13*F(51:100))];
plotResponse(F, H)

目的の IIR 次数 Na = 10 (分母次数) および Nb = 12 (分子次数) を入力として受け取るシングルバンド設計仕様 'Nb, Na, F, H' を使用して、仕様オブジェクトを作成します。この仕様の場合、利用可能な設計法は 1 つしかありません (最小二乗 IIR 設計 iirls)。

Nb = 12;
Na = 10;
f = fdesign.arbmagnphase('Nb,Na,F,H',Nb,Na,F,H);
designmethods(f)

Design Methods for class fdesign.arbmagnphase (Nb,Na,F,H):


iirls

iirls 設計法では異なる周波数に対して別の重みを指定できるため、各帯域の近似品質をより細かく制御できます。阻止帯域に 1 の重み、通過帯域に 100 の重みを付けるようにフィルターを設計します。通過帯域の重みが大きいと、この領域での近似精度が高くなります。

W = 1*(F<=0.5) + 100*(F>0.5);
Hd = design(f,'iirls', Weights = W, SystemObject = true);

IIR 設計法を使用する場合、フィルターの安定性は保証されません。関数 isstable を使用して、IIR の安定性をチェックします。より徹底的な解析を行うには、極を確認し、それらがどれだけ単位円に近いかを調べます。

isstable(Hd)

ans = logical
   1

IIR 設計応答をプロットします。通過帯域の近似は阻止帯域の近似より優れており、振幅ゲインが小さい (dB が低い) 場合は常に位相応答の重要性が低いことに注意してください。

hfvt = fvtool(Hd);
legend(hfvt,'IIR Least-Squares', Location = 'NorthWest')

Figure Figure 3: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 2 objects of type line. This object represents IIR Least-Squares.

hfvt(2) = fvtool(Hd,Analysis='phase');
legend(hfvt(2),'IIR Least-Squares', Location = 'NorthEast')
ax = hfvt(2).CurrentAxes; 
ax.NextPlot = 'add';
plot(ax,F,unwrap(angle(H))+2*pi,'r--')

Figure Figure 4: Phase Response contains an axes object. The axes object with title Phase Response, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Phase (radians) contains 2 objects of type line. This object represents IIR Least-Squares.

低い群遅延のバンドパス FIR 設計

任意の振幅と位相設計の応用として興味深いのは、小さい方の群遅延を優先して線形位相を放棄する非対称 FIR フィルターの設計です。このようなフィルターを設計すると、通過帯域で線形位相の近似を良好に保つことができます。バンドパスフィルターに、阻止帯域 $F_{1} = [0, 0.3)$ と $F_{3} = (0.6, 1]$ 、および通過帯域 $F_{2} = [0.3, 0.6]$ の 3 つの帯域があるとします。通過帯域で、目的の周波数応答は $H (ω) = e^{- j π ω \cdot gd}$ で、これは群遅延を "gd" とする線形位相応答をもちます。

F1 = linspace(0,.25,30);  % Lower stopband
F2 = linspace(.3,.56,40); % Passband
F3 = linspace(.62,1,30);  % Higher stopband

% Define the desired frequency response on the bands
gd = 12; % Desired Group Delay
H1 = zeros(size(F1));
H2 = exp(-1j*pi*F2*gd);
H3 = zeros(size(F3));
F = [F1 F2 F3];
H = [H1 H2 H3];

目的の周波数応答をプロットします。

plotResponse(F, H)

'N,B,F,H' 仕様パターンを使用してフィルター仕様オブジェクトを作成します。ここで、N = 50 は目的のフィルター次数、B = 3 は帯域の数で、B の後は前と同様 F ベクトルと H ベクトルのペアが続きます。

N = 50;  % Filter Order
B = 3;   % Number of bands
f = fdesign.arbmagnphase('N,B,F,H',N,B, F1,H1, F2,H2, F3,H3); 
Hd_mnp = design(f,'equiripple', SystemObject = true);

関数 islinphase を呼び出すとわかるように、この設計には線形位相がありません。

islinphase(Hd_mnp)

ans = logical
   0

ここで、fdesign.arbmag を使用して、振幅のみのフィルター仕様を作成します。このオブジェクトの 'N,B,F,A' 仕様パターンは、fdesign.argmagnphase オブジェクトの 'N,B,F,H' 仕様に似ています。この 2 つの違いは、'N,B,F,H' の複素フィルター応答 H が、'N,B,F,A' の振幅のみ (非負の実数) 応答 A に置き換えられていることです。

f_magonly = fdesign.arbmag('N,B,F,A',N,3,F1,abs(H1),F2,abs(H2),F3,abs(H3)); 
Hd_mo = design(f_magonly,'equiripple', SystemObject = true);

振幅のみの仕様では、線形位相をもつ対称設計が得られます。

islinphase(Hd_mo)

ans = logical
   1

subplot(1,2,1);
stem(Hd_mnp.Numerator)
title('Magnitude and Phase Design (asymetric)')
subplot(1,2,2);
stem(Hd_mo.Numerator)
title('Magnitude-only Design (symmetric)')

Figure contains 2 axes objects. Axes object 1 with title Magnitude and Phase Design (asymetric) contains an object of type stem. Axes object 2 with title Magnitude-only Design (symmetric) contains an object of type stem.

2 つの設計を比較します。両者のバンドパス振幅応答は非常に似通っていることに注意してください。

hfvt = fvtool(Hd_mnp,Hd_mo);
legend(hfvt,'Magnitude and Phase Design (Low Group Delay)', ...
    'Magnitude Only (Linear Phase, High Group Delay)', Location = 'NorthEast')

Figure Figure 5: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 2 objects of type line. These objects represent Magnitude and Phase Design (Low Group Delay), Magnitude Only (Linear Phase, High Group Delay).

群遅延をプロットします。任意の振幅と位相設計では、群遅延が多少異なっています。ただし、差異は小さく、平均で 12.5 サンプルです。この群遅延は、振幅のみ設計の群遅延 (25 サンプル) の半分です。

hfvt(2) = fvtool(Hd_mnp,Hd_mo, Analysis = 'grpdelay');
legend(hfvt(2),'Magnitude and Phase Design (Low Group Delay)', ...
    'Magnitude Only (Linear Phase, High Group Delay)', Location = 'NorthEast')
hfvt(2).zoom([.3 .56 0 35]);

Figure Figure 6: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains 2 objects of type line. These objects represent Magnitude and Phase Design (Low Group Delay), Magnitude Only (Linear Phase, High Group Delay).

群遅延の違いは位相応答にも見られます。勾配が小さくなるほど、群遅延も小さくなることを示しています。

hfvt(2).Analysis = 'phase';
hfvt(2).zoom([.3 .56 -30 10]);

Figure Figure 6: Phase Response contains an axes object. The axes object with title Phase Response, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Phase (radians) contains 2 objects of type line. These objects represent Magnitude and Phase Design (Low Group Delay), Magnitude Only (Linear Phase, High Group Delay).

チェビシェフローパスフィルターの通過帯域のイコライズ

任意の振幅と位相設計の一般的な応用としては、他にも IIR フィルターの非線形位相応答のイコライズがあります。正規化された通過帯域周波数が 1/16 で、通過帯域リップルが 0.5 dB である 3 次のチェビシェフタイプ I ローパスフィルターを考えます。

Fp = 1/16;  % Passband frequency 
Ap = .5;    % Passband ripples
f = fdesign.lowpass('N,Fp,Ap',3,Fp,Ap);
Hcheby = design(f,'cheby1', SystemObject = true); 
close(hfvt(1));
close(hfvt(2));
hfvt = fvtool(Hcheby);
legend(hfvt, 'Chebyshev Lowpass');

Figure Figure 7: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 2 objects of type line. This object represents Chebyshev Lowpass.

群遅延をプロットします。10 ～ 20 サンプルの範囲の群遅延をもつ通過帯域に、群遅延の大きな歪みがあります。

hfvt(2) = fvtool(Hcheby, Analysis = 'grpdelay');
legend(hfvt(2), 'Chebyshev Lowpass');

Figure Figure 5: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains an object of type line. This object represents Chebyshev Lowpass.

群遅延の歪みを軽減する場合は、IIR フィルターの後に FIR イコライザー $H_{eq} (ω)$ を使用します。理論的には、複合フィルターが理想的なローパスです。複合フィルターの通過帯域応答は $G (ω) = H_{ch} (ω) H_{eq} (ω) = e^{- j g_{d} ω}$ であり、平坦な振幅応答および $g_{d}$ 個のサンプルの一定の群遅延に対する振幅リップルを除去します。対象群 $g_{d}$ は因果フィルター設計に割り当てられた FIR 長さに関連しています。この例では、 $g_{d} = 35$ とするのが合理的な選択です。

要約すると、イコライザー設計には 2 つの帯域があります。

通過帯域では、イコライザーの目的となる周波数応答は $H_{eq} (ω) = e^{- j g_{d} ω} / H_{ch} (ω)$ となるはずです。
阻止帯域では、目的となる応答は $H_{eq} (ω) = 0$ となり、 $H_{ch}$ と一致します。

この 2 帯域設計仕様によって、イコライザーの FIR 近似の焦点が通過帯域と阻止帯域にのみ当たるようになります。周波数領域の残りの部分は don't care 領域とみなされます。

gd = 35;    % Passband Group Delay of the equalized filter (linear phase)
F1 = 0:5e-4:Fp; % Passband
D1 = exp(-1j*gd*pi*F1)./freqz(Hcheby,F1*pi); 

Fst = 3/16; % Stopband
F2 = linspace(Fst,1,100); 
D2 = zeros(1,length(F2));

このイコライザー FIR 仕様の実装に使用できる FIR 設計法はいくつかあります。最小二乗設計と等リップル設計の 2 通りの設計法を使用してパフォーマンスを比較します。

feq = fdesign.arbmagnphase('N,B,F,H',51,2,F1,D1,F2,D2);
Heq_ls = design(feq,'firls', SystemObject = true);      % Least-Squares design
Heq_er = design(feq,'equiripple', SystemObject = true); % Equiripple design

% Create the cascaded filters
Gls = cascade(Hcheby,Heq_ls);
Geq = cascade(Hcheby,Heq_er);

両フィルターについてカスケード接続されたシステムの振幅応答をプロットします。

hfvt = fvtool(Hcheby,Gls,Geq);
legend(hfvt,'Chebyshev Lowpass (no equalization)','Least-Squares Equalization (cascade)', ...
    'Equiripple Equalization (cascade)', Location = 'NorthEast')

Figure Figure 6: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

通過帯域付近をズームインします。イコライズの後、通過帯域リップルは最小二乗設計イコライザーの場合に元のフィルターの 0.5 dB から 0.27 dB に減衰し、等リップル設計イコライザーの場合は 0.16 dB に減衰します。

hfvt(2) = fvtool(Hcheby,Gls,Geq);
legend(hfvt(2),'Chebyshev Lowpass (no equalization)', ...
    'Least-Squares Equalization (cascade)', ...
    'Equiripple Equalization (cascade)', Location = 'NorthEast')
hfvt(2).zoom([0 .1 -0.8 .5]);

Figure Figure 8: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Magnitude (dB) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

ここで、位相 (および群遅延) イコライズに目を向けます。複合群遅延は、通過帯域の 35 サンプル周辺 (対象群遅延) ではほぼ一定です。通過帯域外で複合群遅延が発散しているように見えますが、その領域ではフィルターのゲインがゼロになるため、これは重要ではありません。

hfvt(2).Analysis = 'grpdelay';
hfvt(2).zoom([0 1 0 40]);

Figure Figure 8: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

通過帯域付近をズームインします。通過帯域の群遅延は、8.8 サンプルのピークツーピーク差から、最小二乗イコライザーで 0.51 サンプルにイコライズされ、等リップルイコライザーで 0.62 サンプルにイコライズされています。どちらのイコライザーも等しく良好に動作します。

hfvt(3) = fvtool(Hcheby,Gls,Geq,Analysis='grpdelay');
legend(hfvt(3),'Chebyshev Lowpass (no equalization)', ...
    'Least-Squares Equalization (cascade)', ...
    'Equiripple Equalization (cascade)', Location = 'NorthEast')
hfvt(3).zoom([0 Fp 34 36]);

Figure Figure 9: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains 3 objects of type line. These objects represent Chebyshev Lowpass (no equalization), Least-Squares Equalization (cascade), Equiripple Equalization (cascade).

参照

[1] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, 1989.