周波数応答

デジタル領域

freqz では、FFT ベースのアルゴリズムを使用してデジタル フィルターの Z 変換周波数応答が計算されます。特に、以下のステートメント

[h,w] = freqz(b,a,p)

では、デジタル フィルターの p 点の複素周波数応答 H(e^jω) が返されます。

$$H(e^{j\omega })=\frac{b(1)+b(2)e^{-j\omega }+\mathrm{...}+b(n+1)e^{-j\omega n}}{a(1)+a(2)e^{-j\omega }+\mathrm{...}+a(m+1)e^{-j\omega m}}$$

最も簡単な形では、freqz は、フィルター係数ベクトル b と a、および周波数応答を計算する点の数を指定する整数 p を受け取ります。freqz では、ベクトル h に複素周波数応答が返され、ベクトル w に実際の周波数点がラジアン/サンプル単位で返されます。

freqz では、サンプリング周波数や、任意の周波数点のベクトルなどの他のパラメーターを受け取ることもできます。次の例は、12 次のチェビシェフ I 型フィルターに対し 256 点の周波数応答を求めるものです。freqz が呼び出されると、サンプリング周波数 fs に 1000 Hz が指定されます。

[b,a] = cheby1(12,0.5,200/500);
[h,f] = freqz(b,a,256,1000);

パラメーター リストにはサンプリング周波数が含まれるため、freqz では周波数応答計算に使用された 0 から fs/2 の間の 256 点の周波数点を含むベクトル f が返されます。

メモ

このツールボックスでは、サンプリング周波数の半分として定義されるナイキスト周波数を単位周波数としています。すべての基本的なフィルター設計関数のカットオフ周波数パラメーターは、ナイキスト周波数によって正規化されます。たとえば、サンプリング周波数が 1000 Hz のシステムでは、300 Hz は 300/500 = 0.6 です。正規化周波数を単位円での角周波数に変換するには、π で乗算します。正規化周波数を Hz に戻すには、サンプル周波数の 1/2 を乗算します。

出力引数を指定せずに freqz を呼び出した場合は、周波数に対する振幅と位相の両方がプロットされます。たとえば、サンプリング周波数が 2000Hz、カットオフ周波数が 400Hz の 9 次のバタワース ローパス フィルターは以下のようになります。

[b,a] = butter(9,400/1000);

このフィルターについて 256 点の複素周波数応答を計算し、freqz を使用して振幅と位相をプロットするには、以下を使用します。

freqz(b,a,256,2000)

freqz では、任意の周波数点のベクトルを受け取って、周波数領域計算に使用することができます。以下に例を示します。

w = linspace(0,pi);
h = freqz(b,a,w);

ベクトル b と a によって定義されたフィルターについて、w の周波数点における複素周波数応答が計算されます。周波数点は、0 から 2π までの範囲の値を取ります。0 から設定したサンプリング周波数までの範囲の周波数ベクトルを指定するには、周波数ベクトルとサンプリング周波数の値の両方をパラメーター リストで指定します。

以下の例では、デジタル周波数応答を計算して表示する方法を説明します。

伝達関数の周波数応答

次の伝達関数で記述される 3 次の IIR ローパス フィルターの振幅応答を計算し、表示します。

$$H(z)=\frac{0.05634(1+z^{-1})(1-1.0166z^{-1}+z^{-2})}{(1-0.683z^{-1})(1-1.4461z^{-1}+0.7957z^{-2})}.$$

分子と分母を多項式の畳み込みとして表します。単位円全体を 2,001 個に分割する点の周波数応答を求めます。

b0 = 0.05634;
b1 = [1  1];
b2 = [1 -1.0166 1];
a1 = [1 -0.683];
a2 = [1 -1.4461 0.7957];

b = b0*conv(b1,b2);
a = conv(a1,a2);

[h,w] = freqz(b,a,'whole',2001);

dB 単位表記の振幅応答をプロットします。

plot(w/pi,20*log10(abs(h)))
ax = gca;
ax.YLim = [-100 20];
ax.XTick = 0:.5:2;
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')

FIR バンドパス フィルターの周波数応答

$$0.35\pi$$ ～ $$0.8\pi$$ ラジアン/サンプルの通過帯域と 3 dB のリップルをもつ FIR バンドパス フィルターを設計します。最初の阻止帯域は $$0$$ ～ $$0.1\pi$$ ラジアン/サンプルで、減衰量 40 dB をもちます。2 番目の阻止帯域は $$0.9\pi$$ ラジアン/サンプルからナイキスト周波数で、減衰量 30 dB をもちます。周波数応答を計算します。線形単位と dB の両方で振幅をプロットします。通過帯域を強調表示します。

sf1 = 0.1;
pf1 = 0.35;
pf2 = 0.8;
sf2 = 0.9;
pb = linspace(pf1,pf2,1e3)*pi;

bp = designfilt("bandpassfir", ...
    StopbandAttenuation1=40,StopbandFrequency1=sf1, ...
    PassbandFrequency1=pf1,PassbandRipple=3, ...
    PassbandFrequency2=pf2,StopbandFrequency2=sf2, ...
    StopbandAttenuation2=30);

[h,w] = freqz(bp,1024);
hpb = freqz(bp,pb);

subplot(2,1,1)
plot(w/pi,abs(h),pb/pi,abs(hpb),'.-')
axis([0 1 -1 2])
legend("Response","Passband",Location="south")
ylabel("Magnitude")

subplot(2,1,2)
plot(w/pi,db(h),pb/pi,db(hpb),".-")
axis([0 1 -60 10])
xlabel("Normalized Frequency (\times\pi rad/sample)")
ylabel("Magnitude (dB)")

ハイパス フィルターの振幅応答

$$0.5\pi$$ ラジアン/サンプルの正規化された 3 dB 周波数をもつ 3 次ハイパス バタワース フィルターを設計します。その周波数応答を計算します。振幅応答をデシベル単位で表してプロットします。

[b,a] = butter(3,0.5,'high');
[h,w] = freqz(b,a);

dB = mag2db(abs(h));

plot(w/pi,dB)
xlabel('\omega / \pi')
ylabel('Magnitude (dB)')
ylim([-82 5])

アナログ領域

関数 freqs では、2 つの入力係数ベクトル b と a で定義されるアナログ フィルターの周波数応答が計算されます。この操作は、freqz の操作と同様です。使用する周波数点の数を指定し、任意の周波数点のベクトルを指定して、フィルターの振幅応答と位相応答をプロットできます。この例では、アナログ周波数応答を計算して表示する方法を説明します。

アナログの IIR ローパス フィルターの比較

カットオフ周波数が 2 GHz である 5 次のアナログ バタワース ローパス フィルターを設計します。$$2\pi$$ 倍にして周波数を秒あたりのラジアン単位に変換します。4096 点でのフィルターの周波数応答を計算します。

n = 5;
wc = 2*pi*2e9;
w = 2*pi*1e9*logspace(-2,1,4096)';

[zb,pb,kb] = butter(n,wc,"s");
[bb,ab] = zp2tf(zb,pb,kb);
[hb,wb] = freqs(bb,ab,w);
gdb = -diff(unwrap(angle(hb)))./diff(wb);

エッジ周波数が同じで、通過帯域リップルが 3 dB である 5 次のチェビシェフ I 型フィルターを設計します。その周波数応答を計算します。

[z1,p1,k1] = cheby1(n,3,wc,"s");
[b1,a1] = zp2tf(z1,p1,k1);
[h1,w1] = freqs(b1,a1,w);
gd1 = -diff(unwrap(angle(h1)))./diff(w1);

エッジ周波数が同じで、阻止帯域の減衰量が 30 dB である 5 次のチェビシェフ II 型フィルターを設計します。その周波数応答を計算します。

[z2,p2,k2] = cheby2(n,30,wc,"s");
[b2,a2] = zp2tf(z2,p2,k2);
[h2,w2] = freqs(b2,a2,w);
gd2 = -diff(unwrap(angle(h2)))./diff(w2);

エッジ周波数が同じで、通過帯域リップルが 3 dB、阻止帯域の減衰量が 30 dB である 5 次の楕円フィルターを設計します。その周波数応答を計算します。

[ze,pe,ke] = ellip(n,3,30,wc,"s");
[be,ae] = zp2tf(ze,pe,ke);
[he,we] = freqs(be,ae,w);
gde = -diff(unwrap(angle(he)))./diff(we);

エッジ周波数が同じである 5 次のベッセル フィルターを設計します。その周波数応答を計算します。

[zf,pf,kf] = besself(n,wc);
[bf,af] = zp2tf(zf,pf,kf);
[hf,wf] = freqs(bf,af,w);
gdf = -diff(unwrap(angle(hf)))./diff(wf);

減衰をデシベルでプロットします。周波数をギガヘルツで表します。フィルターを比較します。

fGHz = [wb w1 w2 we wf]/(2e9*pi);
plot(fGHz,mag2db(abs([hb h1 h2 he hf])))
axis([0 5 -45 5])
grid on
xlabel("Frequency (GHz)")
ylabel("Attenuation (dB)")
legend(["butter" "cheby1" "cheby2" "ellip" "besself"])

サンプル単位で群遅延をプロットします。周波数をギガヘルツで、群遅延をナノ秒で表します。フィルターを比較します。

gdns = [gdb gd1 gd2 gde gdf]*1e9;
gdns(gdns<0) = NaN;
loglog(fGHz(2:end,:),gdns)
grid on
xlabel("Frequency (GHz)")
ylabel("Group delay (ns)")
legend(["butter" "cheby1" "cheby2" "ellip" "besself"])

バタワース フィルターおよびチェビシェフ II 型フィルターには平坦な通過帯域と広い遷移帯域幅があります。チェビシェフ I 型フィルターおよび楕円フィルターは速くロールオフしますが、通過帯域リップルがあります。チェビシェフ II 型設計関数に対する周波数入力は、通過帯域の末尾ではなく阻止帯域の始点を設定します。ベッセル フィルターは、通過帯域に沿って定数近似群遅延をもちます。