freqz

デジタル フィルターの周波数応答

構文

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

[h,w] = freqz(B,A,"ctf",n)

[h,w] = freqz({B,A,g},"ctf",n)

[h,w] = freqz(d,n)

[h,w] = freqz(sos,n)

[h,w] = freqz(___,n,"whole")

[h,f] = freqz(___,n,fs)

[h,f] = freqz(___,n,"whole",fs)

h = freqz(___,w)

h = freqz(___,f,fs)

freqz(___)

説明

[h,w] = freqz(b,a,n) は、指定されたデジタル フィルターの周波数応答を返します。分子係数 b と分母係数 a をもつデジタル フィルターを指定します。この関数は、n 点周波数応答ベクトルを h で、対応する角周波数ベクトルを w で返します。

[h,w] = freqz(B,A,"ctf",n) は、分子係数 B と分母係数 A をもつCascaded Transfer Functions (CTF) として表されるデジタル フィルターの n 点周波数応答を返します。 (R2024b 以降)

[h,w] = freqz({B,A,g},"ctf",n) は、デジタル フィルターの n 点周波数応答を CTF 形式で返します。分子係数 B、分母係数 A、およびフィルター セクション全体のスケーリング値 g を使用してフィルターを指定します。 (R2024b 以降)

[h,w] = freqz(d,n) では、デジタル フィルター d に対して n 点の複素周波数応答が返されます。

[h,w] = freqz(sos,n) では、2 次セクション行列 sos に対応する n 点の複素周波数応答が返されます。

[h,w] = freqz(___,n,"whole") は、単位円全体を n 個に分割するサンプル点の周波数応答を返します。

[h,f] = freqz(___,n,fs) では、レート fs でサンプリングした信号をフィルター処理するために設計されたデジタル フィルターに対して周波数応答ベクトル h と対応する物理周波数ベクトル f が返されます。

[h,f] = freqz(___,n,"whole",fs) は、0 から fs の範囲の n 個の点における周波数ベクトルを返します。

h = freqz(___,w) では、w で提供される正規化周波数で評価された周波数応答ベクトル h が返されます。

h = freqz(___,f,fs) では、f で提供される物理周波数で評価された周波数応答ベクトル h が返されます。

出力引数なしで freqz(___) を使用すると、フィルターの周波数応答がプロットされます。

例

伝達関数の周波数応答

次の伝達関数で記述される 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)')

カスケード伝達関数の周波数応答

R2024b 以降

阻止帯域エッジ周波数が 0.4、阻止帯域の減衰量が 50 dB の 40 次ローパス チェビシェフ II 型デジタル フィルターを設計します。CTF 形式の係数を使用して、フィルターの周波数応答をプロットします。

[B,A] = cheby2(40,50,0.4,"ctf");

freqz(B,A,"ctf")

通過帯域エッジ周波数が 0.3 と 0.7、通過帯域リップルが 0.1 dB、阻止帯域の減衰量が 50 dB の 30 次バンドパス楕円デジタル フィルターを設計します。CTF 形式の係数とゲインを使用して、フィルターの周波数応答をプロットします。

[B,A,g] = ellip(30,0.1,50,[0.3 0.7],"ctf");
freqz({B,A,g},"ctf")

FIR フィルターの周波数応答

$$\beta =8$$ のカイザー ウィンドウを使用して次数 80 の FIR ローパス フィルターを設計します。正規化されたカットオフ周波数、$$0.5\pi$$ ラジアン/サンプルを指定します。フィルターの振幅応答と位相応答を表示します。

b = fir1(80,0.5,kaiser(81,8));
freqz(b,1)

designfilt を使用して同じフィルターを設計します。その振幅応答と位相応答を表示します。

d = designfilt("lowpassfir",FilterOrder=80, ...
    CutoffFrequency=0.5,Window={"kaiser",8});
freqz(d)

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)")

2 次セクションの周波数応答

次の伝達関数で記述される 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 次セクションで伝達関数を表します。単位円全体を 2,001 個に分割する点の周波数応答を求めます。

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

sos1 = [b0*[b1 0] [a1 0]];
sos2 = [b2 a2];

[h,w] = freqz([sos1;sos2],'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)')

入力引数

b, a — 伝達関数の係数

ベクトル

伝達関数の係数。ベクトルで指定します。b と a によるこの伝達関数式は次になります。

$$H(z)=\frac{B(z)}{A(z)}=\frac{b_1+b_2{z}^{-1}\cdots +b_n{z}^{-\left(n-1\right)}+b_{n+1}{z}^{-n}}{a_1+a_2{z}^{-1}\cdots +a_m{z}^{-\left(m-1\right)}+a_{m+1}{z}^{-m}}$$

例: b = [1 3 3 1]/6 と a = [3 0 1 0]/3 は、正規化された 3 dB の周波数 0.5π ラジアン/サンプルをもつ 3 次のバタワース フィルターを指定します。

データ型: double | single
複素数のサポート: あり

n — 応答を求める周波数点の数

512 (既定値) | 正の整数スカラー

応答を求める周波数点の数。2 以上の正の整数スカラーとして指定します。n がない場合の既定値は 512 です。最適な結果を得るため、n をフィルター次数よりも大きい値に設定します。

B,A — カスケード伝達関数 (CTF) 係数

スカラー | ベクトル | 行列

R2024b 以降

カスケード伝達関数 (CTF) 係数。スカラー、ベクトル、または行列として指定します。B と A は、それぞれカスケード伝達関数の分子係数と分母係数をリストします。

B のサイズは L 行 (m + 1) 列、A のサイズは L 行 (n + 1) でなければなりません。ここで、

• L はフィルター セクションの数を表します。

• m はフィルターの分子の次数を表します。

• n はフィルターの分母の次数を表します。

カスケード伝達関数の形式と係数行列の詳細については、CTF 形式によるデジタル フィルターの指定を参照してください。

メモ

A(:,1) のいずれかの要素が 1 と等しくない場合、freqz はフィルター係数を A(:,1) で正規化します。この場合、A(:,1) はゼロ以外でなければなりません。

データ型: double | single
複素数のサポート: あり

g — スケール値

スカラー | ベクトル

R2024b 以降

スケール値。実数値のスカラー、または L + 1 個の要素をもつ実数値のベクトルとして指定します。ここで、L は CTF セクションの数です。スケール値は、カスケード フィルター表現のセクション全体にわたるフィルター ゲインの分布を表します。

freqz 関数は、g の指定方法に応じて、scaleFilterSections 関数を使用してフィルター セクションにゲインを適用します。

• スカラー — この関数は、すべてのフィルター セクションにわたってゲインを均一に配分します。

• ベクトル — この関数は、最初の L 個のゲイン値を対応するフィルター セクションに適用し、最後のゲイン値をすべてのフィルター セクションに均一に配分します。

データ型: double | single

d — デジタル フィルター

digitalFilter オブジェクト

デジタル フィルター。digitalFilter オブジェクトで指定します。デジタル フィルターを周波数応答仕様に基づいて生成するには、関数 designfilt を使用します。

例: d = designfilt('lowpassiir','FilterOrder',3,'HalfPowerFrequency',0.5) は、正規化された 3 dB の周波数 0.5π ラジアン/サンプルをもつ 3 次のバタワース フィルターを指定します。

sos — 2 次セクションの係数

行列

2 次セクションの係数。行列として指定します。sos は K 行 6 列の行列で、セクション数 K が 2 以上でなければなりません。セクション数が 2 未満の場合、関数は入力を分子ベクトルとして扱います。sos の各行は 2 次 (双二次) フィルターの係数に対応しています。sos の i 行目は [bi(1) bi(2) bi(3) ai(1) ai(2) ai(3)] に対応しています。

例: s = [2 4 2 6 0 2;3 3 0 6 0 0] は、正規化された 3 dB の周波数 0.5π ラジアン/サンプルをもつ 3 次のバタワース フィルターを指定します。

データ型: double | single
複素数のサポート: あり

fs — サンプル レート

正のスカラー

サンプル レート。正のスカラーで指定します。時間の単位が秒の場合、fs は Hz で表されます。

データ型: double

w — 角周波数

ベクトル

角周波数。ラジアン/サンプル単位のベクトルとして指定します。w は少なくとも 2 つの要素をもたなければなりません。そうでない場合は関数が n として解釈するためです。w = π はナイキスト周波数に相当します。

f — 周波数

ベクトル

周波数。ベクトルとして指定します。f は少なくとも 2 つの要素をもたなければなりません。そうでない場合は、関数が n として解釈するためです。時間の単位が秒の場合、f は Hz で表されます。

データ型: double

出力引数

h — 周波数応答

ベクトル

ベクトルとして返される周波数応答。n を指定している場合は、h の長さは n になります。n を指定していないか、n を空ベクトルとして指定している場合には、h の長さは 512 になります。

freqz への入力が単精度の場合、関数によって周波数応答が単精度演算で計算されます。出力 h は単精度です。

w — 角周波数

ベクトル

ベクトルとして返される角周波数。w は 0 ～π の値を取ります。入力に 'whole' を指定すると、w の値の範囲は 0 ～ 2π になります。n を指定している場合は、w の長さは n になります。n を指定していないか、n を空ベクトルとして指定している場合には、w の長さは 512 になります。

f — 周波数

ベクトル

周波数。ヘルツ表記でベクトルとして返されます。f は 0 ～ fs/2 Hz の範囲の値を取ります。入力に 'whole' を指定している場合、f の値の範囲は 0 ～ fs Hz になります。n を指定している場合は、f の長さは n になります。n を指定していないか、n を空ベクトルとして指定している場合には、f の長さは 512 になります。

詳細

カスケード伝達関数

IIR デジタル フィルターをカスケード セクションに分割すると、数値安定性が向上し、係数の量子化誤差の影響を受けにくくなります。L 個の伝達関数 H₁(z)、H₂(z)、…、H_L(z) による伝達関数 H(z) のカスケード形式は次のとおりです。

$$H(z)={\displaystyle \prod _{l=1}^L{H}_l(z)}=H_1(z)\times H_2(z)\times \cdots \times H_L(z).$$

CTF 形式によるデジタル フィルターの指定

解析、可視化、信号フィルター処理で使用するため、CTF 形式でデジタル フィルターを指定できます。係数 B と A をリストしてフィルターを指定します。スカラーまたはベクトルの g を指定することで、セクション全体にフィルターのスケーリング ゲインを含めることもできます。

フィルター係数

係数を L 行の行列

$$B=\left[\begin{array}{cccc}{b}_{11}& b_{12}& \cdots & b_{1,m+1}\\ b_{21}& b_{22}& \cdots & b_{2,m+1}\\ ⋮& ⋮& \ddots & ⋮\\ b_{L1}& b_{L2}& \cdots & b_{L,m+1}\end{array}\right],A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1,n+1}\\ a_{21}& a_{22}& \cdots & a_{2,n+1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{L1}& a_{L2}& \cdots & a_{L,n+1}\end{array}\right],$$

として指定した場合、L 個のカスケード伝達関数のシーケンスとしてフィルターを指定したものとみなされます。このとき、フィルターの完全な伝達関数は次のようになります。

$$H\left(z\right)=\frac{b_{11}+b_{12}{z}^{-1}+\cdots +b_{1,m+1}{z}^{-m}}{a_{11}+a_{12}{z}^{-1}+\cdots +a_{1,n+1}{z}^{-n}}\times \frac{b_{21}+b_{22}{z}^{-1}+\cdots +b_{2,m+1}{z}^{-m}}{a_{21}+a_{22}{z}^{-1}+\cdots +a_{2,n+1}{z}^{-n}}\times \cdots \times \frac{b_{L1}+b_{L2}{z}^{-1}+\cdots +b_{L,m+1}{z}^{-m}}{a_{L1}+a_{L2}{z}^{-1}+\cdots +a_{L,n+1}{z}^{-n}},$$

ここで、m ≥ 0 はフィルターの "分子次数"、n ≥ 0 は "分母次数" です。

• B と A の両方をベクトルとして指定した場合、基になるシステムは 1 つのセクションから成る IIR フィルター (L = 1) であるとみなされます。このとき、B は伝達関数の分子を表し、A はその分母を表します。

• B がスカラーの場合、フィルターは全極 IIR フィルターのカスケードであるとみなされます。このとき、各セクションの全体的なシステム ゲインは B に等しくなります。

• A がスカラーの場合、フィルターは FIR フィルターのカスケードであるとみなされます。このとき、各セクションの全体的なシステム ゲインは 1/A に等しくなります。

メモ

• 2 次セクションの行列をカスケード伝達関数に変換するには、sos2ctf 関数を使用します。

• 零点-極-ゲイン フィルター表現をカスケード伝達関数に変換するには、zp2ctf 関数を使用します。

係数とゲイン

全体的なスケーリング ゲイン、または係数値から共通因数としてくくりだされた複数のスケーリング ゲインがある場合、係数とゲインを {B,A,g} 形式の cell 配列として指定できます。フィルター セクションのスケーリングは、固定小数点演算を使用して各フィルター セクションの出力の振幅レベルを確実に同等にする場合に特に重要であり、数値精度の制限によるフィルター応答の不正確さを回避するのに役立ちます。

ゲインは、全体的なゲインを表すスカラー、またはセクションのゲインを表すベクトルとして指定できます。

• ゲインがスカラーの場合、その値はすべてのカスケード フィルター セクションに均一に適用されます。

• ゲインがベクトルの場合、カスケード内のフィルター セクションの数 L よりも 1 つ多い要素をもたなければなりません。最初の L 個のスケール値はそれぞれ対応するフィルター セクションに適用され、最後の値はすべてのカスケード フィルター セクションに均一に適用されます。

係数行列とゲイン ベクトルを

$$B=\left[\begin{array}{cccc}{b}_{11}& b_{12}& \cdots & b_{1,m+1}\\ b_{21}& b_{22}& \cdots & b_{2,m+1}\\ ⋮& ⋮& \ddots & ⋮\\ b_{L1}& b_{L2}& \cdots & b_{L,m+1}\end{array}\right],A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1,n+1}\\ a_{21}& a_{22}& \cdots & a_{2,n+1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{L1}& a_{L2}& \cdots & a_{L,n+1}\end{array}\right],g=\left[\begin{array}{ccccc}{g}_1& g_2& \cdots & g_L& g_{\text{S}}\end{array}\right],$$

として指定した場合、フィルター システムの伝達関数は次であるとみなされます。

$$H\left(z\right)=g_{\text{S}}\left(g_1\frac{b_{11}+b_{12}{z}^{-1}+\cdots +b_{1,m+1}{z}^{-m}}{a_{11}+a_{12}{z}^{-1}+\cdots +a_{1,n+1}{z}^{-n}}\times g_2\frac{b_{21}+b_{22}{z}^{-1}+\cdots +b_{2,m+1}{z}^{-m}}{a_{21}+a_{22}{z}^{-1}+\cdots +a_{2,n+1}{z}^{-n}}\times \cdots \times g_L\frac{b_{L1}+b_{L2}{z}^{-1}+\cdots +b_{L,m+1}{z}^{-m}}{a_{L1}+a_{L2}{z}^{-1}+\cdots +a_{L,n+1}{z}^{-n}}\right).$$

ヒント

• スケーリング ゲインを含むフィルターを CTF 形式で取得できます。butter、cheby1、cheby2、ellip などのデジタル IIR フィルター設計関数の出力を使用します。これらの関数で、"ctf" フィルタータイプ引数を指定し、スケール値を取得するために B、A、g を返すように指定します。 (R2024b 以降)

• 既約マルチレート フィルターがある場合は、関数 freqzmr (DSP System Toolbox) を使用して周波数領域でフィルターを解析します。既約マルチレート フィルターの詳細については、Overview of Multirate Filters (DSP System Toolbox)を参照してください。 (R2024a 以降)

• 関数 freqzmr (DSP System Toolbox) には DSP System Toolbox™ が必要です。 (R2024a 以降)

アルゴリズム

デジタル フィルターの周波数応答は、z = e^jω で計算した伝達関数として解釈されます[1]。

freqz では、ユーザーが指定した実数または複素数の分子多項式や分母多項式から、伝達関数が決定されます。次に、デジタル フィルターの複素周波数応答 H (e^jω) が返されます。周波数応答は、使用した構文によって決定されたサンプル点で計算されます。

freqz では、入力引数として周波数ベクトルを指定しなかった場合、通常は FFT アルゴリズムを使用して周波数応答が計算されます。周波数応答は伝達された分子係数と分母係数の比として計算され、希望する長さになるように 0 が付加されます。

周波数ベクトルを入力として指定すると、freqz では、各周波数点で多項式が計算され、分子の応答が分母の応答で除算されます。多項式を計算するために、関数はホーナー法を使用します。

参照

[1] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

[2] Lyons, Richard G. Understanding Digital Signal Processing. Upper Saddle River, NJ: Prentice Hall, 2004.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• freqz 関数に CTF 係数を入力する場合、コードを生成するには、freqz 構文に "ctf" 入力引数を含めなければなりません。

• freqz への最初の入力は、コンパイル時に可変サイズの行列である場合、実行時にベクトルまたは空の配列にしてはなりません。

• 入力 n は、コンパイル時に可変サイズである場合、実行時にスカラーまたは空の配列にしてはなりません。

バージョン履歴

R2006a より前に導入

R2025a: [プロットの作成] ライブ エディター タスクを使用したカスケード伝達関数の可視化

[プロットの作成] ライブ エディター タスクでカスケード伝達関数 (CTF) を使用して、freqz の出力を対話的に可視化できます。

R2024b: カスケード伝達関数を使用したフィルターの解析

freqz 関数は、カスケード伝達関数 (CTF) 形式の入力をサポートします。

R2023a: [プロットの作成] ライブ エディター タスクを使用した関数出力の可視化

[プロットの作成] ライブ エディター タスクを使用すると、freqz の出力を対話的に可視化できます。さまざまなグラフ タイプを選択し、オプションのパラメーターを設定できます。タスクは、ライブ スクリプトの一部となるコードも自動的に生成します。

参考

アプリ

• フィルター アナライザー | フィルター デザイナー

関数

• abs | angle | ctffilt | designfilt | digitalFilter | fft | filter | freqs | impz | invfreqs | logspace | freqzmr (DSP System Toolbox)