FIR ハーフバンドフィルターの設計

ライブスクリプトを開く

この例では、FIR ハーフバンドフィルターの設計法を示します。ハーフバンドフィルターは、マルチレート信号処理アプリケーションにおいて 2 の係数で内挿または間引きする際に幅広く使用されます。ハーフバンドフィルター係数のほぼ半分はゼロに等しくなるため、多くの場合、ハーフバンドフィルターはポリフェーズ形式で効率的に実装できます。

ハーフバンドフィルターには 2 つの重要な特徴があります。

通過帯域リップルと阻止帯域リップルは同じでなければなりません。
通過帯域エッジと阻止帯域エッジの周波数はハーフバンド周波数 $\frac{Fs}{4}$ Hz (または正規化周波数の $\frac{π}{2}$ ラジアン/サンプル) から等距離になければなりません。

関数 `designHalfbandFIR` を使用したハーフバンドフィルターの設計

関数designHalfbandFIRを使用して、FIR ハーフバンドフィルターを設計します。この関数は、フィルター係数またはサポートされている FIR フィルターオブジェクトの 1 つを返します。

等リップルハーフバンドフィルターの設計

通過帯域周波数エッジを 22 kHz とし、96 kHz でサンプリングされたデータで処理を行うハーフバンドフィルターを考えます。そのサンプルレートに対応するハーフバンド周波数は 24 kHz で、これはサンプルレートの 4 分の 1 に相当します。この設計の遷移幅は 2 kHz、つまり正規化単位で 1/12 ラジアン/秒です。DesignMethod を "equiripple" に設定します。

Fs  = 96e3; % Sample rate
Fp  = 22e3; % Passband edge
Fh  = Fs/4; % Halfband frequency

TW = (Fh-Fp)/Fh; % Normalized transition width
N   = 100; % Filter order

num = designHalfbandFIR(FilterOrder=N,TransitionWidth=TW,Passband="lowpass",DesignMethod="equiripple");
zerophase(num,1,linspace(0,Fs/2,512),Fs);

Figure contains an axes object. The axes object with title Zero-Phase Response, xlabel Frequency (kHz), ylabel Amplitude contains an object of type line.

応答でズームインすることで、通過帯域と阻止帯域のピークツーピークリップルが同じであることを確認できます。 $\frac{Fs}{4}$ (24 kHz) 点を中心とした対称性があります。指定したとおりに通過帯域が最大 22 kHz まで拡張し、26 kHz で阻止帯域が始まります。インパルス応答を調べることで、それ以外のすべての係数が 0 に等しいことを検証することもできます。これにより、係数 2 で内挿または間引きするためのフィルター実装を大変効率よく行えます。

impz(num);
xlim([30 70])

Figure contains an axes object. The axes object with title Impulse Response, xlabel n (samples), ylabel Amplitude contains an object of type stem.

ハーフバンドフィルターによるストリーミングデータのフィルター処理

ストリーミングデータを処理する場合は、効率的なフィルター実装を提供するdsp.FIRFilterSystem object、dsp.FIRHalfbandInterpolatorSystem object、およびdsp.FIRHalfbandDecimatorSystem object を使用します。これらの System object は、倍精度や単精度の浮動小数点データだけでなく固定小数点データのフィルター処理もサポートしています。また、C コードおよび HDL コード生成に加えて、ARM® Cortex® M および ARM® Cortex® A に最適化されたコード生成もサポートします。dsp.FIRFilter を使用する場合は、Structure を 'single-rate' に設定して designHalfbandFIR を使用します。dsp.FIRHalfbandInterpolator および dsp.FIRHalfbandDecimator を使用する場合は、Structure をそれぞれ 'interp' または 'decim' に設定します。

FIR ハーフバンドシングルレートフィルターを作成します。

num = designHalfbandFIR(FilterOrder=N,TransitionWidth=TW,Passband="lowpass",DesignMethod="equiripple",Structure="single-rate");
halfbandFilter = dsp.FIRFilter(num)

halfbandFilter = 
  dsp.FIRFilter with properties:

            Structure: 'Direct form'
      NumeratorSource: 'Property'
            Numerator: [0 2.1118e-04 0 -2.4012e-04 0 3.7199e-04 0 -5.4746e-04 0 7.7546e-04 0 -0.0011 0 0.0014 0 -0.0019 0 0.0024 0 -0.0031 0 0.0039 0 -0.0049 0 0.0060 0 -0.0074 0 0.0090 0 -0.0110 0 0.0134 0 -0.0163 0 0.0201 0 -0.0252 0 … ] (1×101 double)
    InitialConditions: 0

  Show all properties

SystemObject パラメーターを true に設定すると、designHalfbandFIR で直接 System object を作成できます。

FIR ハーフバンド内挿器を作成します。

halfbandInterpolator = designHalfbandFIR(FilterOrder=N, TransitionWidth=TW,...
                            Passband="lowpass",DesignMethod="equiripple",...
                            Structure="interp",SystemObject=true);
freqz(halfbandInterpolator,1024,Fs);

Figure contains 2 axes objects. Axes object 1 with title Phase, xlabel Frequency (kHz), ylabel Phase (degrees) contains an object of type line. Axes object 2 with title Magnitude, xlabel Frequency (kHz), ylabel Magnitude (dB) contains an object of type line.

`dsp.FIRHalfbandInterpolator` および `dsp.FIRHalfbandDecimator` の自己設計モード

dsp.FIRHalfbandInterpolator System object および dsp.FIRHalfbandDecimator System object は、フィルター係数を内部的に設計できます。このオプションを使用するには、dsp.FIRHalfbandInterpolator System object または dsp.FIRHalfbandDecimator System object を作成し、使用する設計仕様を Specification プロパティで指定します。Specification が 'Coefficients' に設定されていない場合、オブジェクトは自己設計されます。オブジェクトは、関数 designHalfbandFIR と同じパラメーター名 (FilterOrder、TransitionWidth、および StopbandAttenuation) を使用します。自己設計モードを使用するもう 1 つの利点は、周波数単位の正規化に加えて、任意のサンプルレートのサポートです。

サンプルレート 96 kHz、フィルター次数 50、遷移幅 4 kHz でハーフバンド内挿器を作成します。

Fs = 96e3; % Sample rate
TW = 4e3;  % Transition width
N  = 50;   % Filter order

halfbandInterpolator = dsp.FIRHalfbandInterpolator(SampleRate=Fs,...
    Specification="Filter order and transition width",...
    FilterOrder=N,TransitionWidth=TW);

freqz(halfbandInterpolator,1024,Fs);

オブジェクトがロックされていない限り、設計パラメーターを調整することができます。また、設計はそれに応じて変更されます。

halfbandInterpolator.FilterOrder=N*2;
freqz(halfbandInterpolator,1024,Fs);

関数tfを使用して、FIR ハーフバンド内挿器またはデシメーター System object の係数を抽出します。上の例のフィルター次数は、予想どおり 100 になります。

num = tf(halfbandInterpolator);
length(num)

ans = 101

ハーフバンド内挿を実行するには、入力データを dsp.FIRHalfbandInterpolator System object 経由で渡します。

FrameSize = 256;
sine1 = dsp.SineWave(Frequency=10e3,SampleRate=Fs,SamplesPerFrame=FrameSize);
sine2 = dsp.SineWave(Frequency=20e3,SampleRate=Fs,SamplesPerFrame=FrameSize);

x = sine1() + sine2() + 0.01.*randn(FrameSize,1); % Input signal

y = halfbandInterpolator(x); % Step through the object

フィルターの遅延を補正することで、入力サンプルに内挿サンプルを重ねてプロットします。関数outputDelayを使用して、フィルター遅延を取得します。ハーフバンドフィルターのポリフェーズ分岐の 1 つが入力サンプルを変化させない純粋な遅延分岐であるため、入力サンプルはフィルターの出力において変わっていません。

[D, FsOut] = halfbandInterpolator.outputDelay();
nx = 0:length(x)-1;
ny = 0:length(y)-1;
plot(nx/Fs+D,x,".",MarkerSize=10)
hold on
stem(ny/FsOut,y)
hold off
legend("Input samples","Interpolated samples")
xlim([1e-3, 1.4e-3])

Figure contains an axes object. The axes object contains 2 objects of type line, stem. One or more of the lines displays its values using only markers These objects represent Input samples, Interpolated samples.

spectrumAnalyzerを使用して、出力信号のスペクトル成分をプロットします。内挿の場合には、2 倍にアップサンプリングした後にフィルター処理を行います (概念上)。したがって、2 でアップサンプリングされるため、spectrumAnalyzer のサンプルレートは $2 Fs$ として指定する必要があります。

scope = spectrumAnalyzer(SampleRate=2*Fs);

tic
while toc < 10
    x = sine1() + sine2() + 0.01.*randn(FrameSize,1); %  96 kHz
    y = halfbandInterpolator(x);                      % 192 kHz
    scope(y);
end

release(scope);

スペクトルの複製は約 40 dB 減衰します。これは、ハーフバンドフィルターによって生じるおおよその減衰量です。

間引きの場合、オブジェクトはフィルター処理後にダウンサンプリング (概念上) するため、dsp.FIRHalfbandDecimator では、入力のサンプルレートに対応するようにサンプルレートを指定します。ハーフバンドデシメーターを前のセクションで設計したハーフバンド内挿器とカスケードし、出力のスペクトル成分をプロットします。2 つのブロックの出力のサンプルレートは、入力と同様に Fs です。

halfbandDecimator = dsp.FIRHalfbandDecimator(SampleRate=2*Fs,...
    Specification="Filter order and transition width",...
    FilterOrder=N,TransitionWidth=4000);

scope = spectrumAnalyzer(SampleRate=Fs);
tic
while toc < 10
    x = sine1() + sine2() + 0.01.*randn(FrameSize,1); %  96 kHz
    y = halfbandInterpolator(x);                      % 192 kHz
    xd = halfbandDecimator(y);                        %  96 kHz again
    scope(xd);
end

release(scope);

最小次数設計仕様

フィルター次数を指定する代わりに、遷移幅と阻止帯域の減衰量を指定できます。ハーフバンドフィルター設計アルゴリズムは、仕様を満たす最小のフィルター次数を自動的に求めようとします。

Ast = 80; % 80 dB
num = designHalfbandFIR(StopbandAttenuation=Ast,TransitionWidth=4000/Fs,DesignMethod="equiripple");
minimumOrder = length(num)

minimumOrder = 223

freqz(num,1,1024,Fs)

dsp.FIRHalfbandInterpolator オブジェクトと dsp.FIRHalfbandDecimator オブジェクトでも同じことができます。ただし、Specification プロパティを明示的に "Transition width and stopband attenuation" に設定する必要があります。

halfbandInterpolator = dsp.FIRHalfbandInterpolator(SampleRate=Fs,Specification="Transition width and stopband attenuation",...
            StopbandAttenuation=Ast, TransitionWidth=4000);
freqz(halfbandInterpolator)

$Figure contains 2 axes objects. Axes object 1 with title Phase, xlabel Normalized Frequency (\times\pi rad/sample), ylabel Phase (degrees) contains an object of type line. Axes object 2 with title Magnitude, xlabel Normalized Frequency (\times\pi rad/sample), ylabel Magnitude (dB) contains an object of type line.$

フィルター次数と阻止帯域の減衰量の指定

フィルター次数と阻止帯域の減衰量を指定してフィルターを設計することもできます。

num = designHalfbandFIR(StopbandAttenuation=Ast,FilterOrder=N,DesignMethod="equiripple");
freqz(num,1,1024,Fs);

dsp.FIRHalfbandInterpolator オブジェクトと dsp.FIRHalfbandDecimator オブジェクトでも同じことができます。ただし、Specification プロパティを明示的に "Filter order and stopband attenuation" に設定する必要があります。

halfbandDecimator = dsp.FIRHalfbandDecimator(SampleRate=Fs,...
    Specification="Filter order and stopband attenuation",...
    StopbandAttenuation=Ast,FilterOrder=N);
fvtool(halfbandDecimator,Fs=Fs);

Figure Figure 1: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (kHz), ylabel Magnitude (dB) contains 2 objects of type line.

フィルターバンクでのハーフバンドフィルターの使用

ハーフバンド内挿器とデシメーターを使用して、合成フィルターバンクと解析フィルターバンクを効率的に実装できます。この例でこれまでに説明したハーフバンドフィルターはすべてローパスフィルターでした。1 つの加算器を追加することで、ローパス応答に加えてハイパス応答も得ることができ、この 2 つの応答を使用してフィルターバンクを実装できます。

直交ミラーフィルター (QMF) バンクをシミュレートします。まず、ローパスおよびハイパスのハーフバンドデシメーターを使用して、1 kHz と 3 kHz の正弦波で構成される 8 kHz 信号を 2 つのハーフレート信号 (4 kHz) に分割します。ローパス信号は 1 kHz の正弦波を保持し、ハイパス信号は 3 kHz の正弦波 (ダウンサンプリング後に 1 kHz にエイリアシングされる) を保持します。ハーフバンド内挿器を使用した合成フィルターバンクで信号を再統合します。ハイパス分岐により、エイリアシングされた 1 kHz の正弦波がアップコンバートされ、3 kHz に戻ります。内挿された信号のサンプルレートは 8 kHz です。

Fs1 = 8000; % Units = Hz
Spec = "Filter order and transition width";
Order = 52;
TW = 4.1e2; % Units = Hz

% Construct FIR Halfband Interpolator
halfbandInterpolator = dsp.FIRHalfbandInterpolator( ...
    Specification=Spec,...
    FilterOrder=Order,...
    TransitionWidth=TW,...
    SampleRate=Fs1/2,...
    FilterBankInputPort=true);

% Construct FIR Halfband Decimator
halfbandDecimator = dsp.FIRHalfbandDecimator( ...
    Specification=Spec,...
    FilterOrder=Order,...
    TransitionWidth=TW,...
    SampleRate=Fs1);

% Input
f1 = 1000;
f2 = 3000;
InputWave = dsp.SineWave(Frequency=[f1,f2],SampleRate=Fs1,...
    SamplesPerFrame=1024,Amplitude=[1 0.25]);

% Construct Spectrum Analyzer object to view the input and output
scope = spectrumAnalyzer(SampleRate=Fs1,...
    PlotAsTwoSidedSpectrum=false,ShowLegend=true,...
    YLimits=[-120 30],...
    Title="Input Signal and Output Signal of Quadrature Mirror Filter",...
    ChannelNames={"Input","Output"});

tic
while toc < 10
    Input = sum(InputWave(),2);
    NoisyInput = Input+(10^-5)*randn(1024,1);
    [Lowpass,Highpass] = halfbandDecimator(NoisyInput);
    Output = halfbandInterpolator(Lowpass,Highpass);
    scope([NoisyInput,Output]);
end

release(scope);

カイザーウィンドウ設計

この例でこれまでに示したすべての設計は最適な等リップル設計でした。関数 designHalfbandFIR に加え、dsp.FIRHalfbandDecimator System object および dsp.FIRHalfbandInterpolator System object も、カイザーウィンドウ法を使用してハーフバンドフィルターを設計できます。

等リップルウィンドウフィルター設計とカイザーウィンドウフィルター設計を比較します。前の手順で使用したものと同じ仕様のサンプルレート、フィルター次数、および遷移帯域幅を使用します。

Fs = 44.1e3;
N  = 90;
TW = 1000;
equirippleHBFilter = designHalfbandFIR(DesignMethod="equiripple",...
    TransitionWidth=TW/Fs,FilterOrder=N,...
    SystemObject=true); 
kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",...
    TransitionWidth=TW/Fs,FilterOrder=N,...
    SystemObject=true);

FVTool を使用して設計を比較します。この 2 つの設計により、阻止帯域の最小減衰量とより大きい全体的な減衰量のトレードオフが可能になります。

fvt = fvtool(equirippleHBFilter,kaiserHBFilter,Fs=2*Fs);
legend(fvt,"Equiripple design","Kaiser-window design")

Figure Figure 2: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (kHz), ylabel Magnitude (dB) contains 3 objects of type line. These objects represent Equiripple design, Kaiser-window design.

カイザー設計での阻止帯域の減衰量の指定

代わりに次数と阻止帯域の減衰量を指定することもできます。これにより、全体的な阻止帯域の減衰量と遷移幅とのトレードオフが可能になります。

Ast  = 60; % Minimum stopband attenuation

equirippleHBFilter = designHalfbandFIR(DesignMethod="equiripple",...
    StopbandAttenuation=Ast,FilterOrder=N,...
    SystemObject=true); 
kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",...
    StopbandAttenuation=Ast,FilterOrder=N,...
    SystemObject=true); 
fvt = fvtool(equirippleHBFilter,kaiserHBFilter,Fs=2*Fs);
legend(fvt,"Equiripple design","Kaiser-window design")

Figure Figure 3: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (kHz), ylabel Magnitude (dB) contains 3 objects of type line. These objects represent Equiripple design, Kaiser-window design.

最小次数のカイザー設計

カイザーウィンドウ設計法は、最小次数設計もサポートしています。ターゲット遷移幅と阻止帯域の減衰量を指定すると、設計アルゴリズムは、より小さい次数のカイザー FIR ハーフバンドをターゲットの仕様を満たす範囲で求めようとします。最小次数のカイザー設計は通常、等リップル設計よりも次数が大きくなりますが、その代わりに阻止帯域の減衰量が全体的に優れています。

Fs = 44.1e3;
TW = 1000; % Transition width
Ast = 60;  % 60 dB minimum attenuation in the stopband

equirippleHBFilter = designHalfbandFIR(DesignMethod="equiripple",...
    StopbandAttenuation=Ast,TransitionWidth=TW/Fs,...
    SystemObject=true); 
kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",...
    StopbandAttenuation=Ast,TransitionWidth=TW/Fs,...
    SystemObject=true); 

fvt = fvtool(equirippleHBFilter,kaiserHBFilter);
legend(fvt,"Equiripple design","Kaiser-window design")

Figure Figure 4: 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 Equiripple design, Kaiser-window design.

フィルター設計法の自動設定

カイザー法と等リップル法は設計仕様によって強みが異なります。非常に高い阻止帯域の減衰量や非常に狭い遷移幅などの厳しい設計仕様の場合、等リップル設計法は収束しないことがあります。このような場合には、カイザー法が優れています。

次のコードは、フィルターの仕様が厳しすぎるため、等リップル設計を使用できない場合を示しています。DesignMethod プロパティを "Equiripple" に設定すると、フィルターの周波数応答からわかるように、設計が収束しません。

TW = 0.001; 
Ast = 180;  % 180 dB minimum attenuation in the stopband
equirippleHBFilter = designHalfbandFIR(DesignMethod="equiripple",...
    TransitionWidth=TW,...
    StopbandAttenuation=Ast,...
    SystemObject=true);
fvt = fvtool(equirippleHBFilter);
legend(fvt,"DesignMethod = equiripple");

DesignMethod を "kaiser" に設定して再度設計します。カイザーベースのフィルター設計は厳密な設計仕様を正確に満たしていませんが、2 つのフィルターの周波数応答を比較するとわかるように、こちらのフィルターはより適切に収束します。

kaiserHBFilter = designHalfbandFIR(DesignMethod="kaiser",...
    TransitionWidth=TW,...
    StopbandAttenuation=Ast,...
    SystemObject=true);
fvt = fvtool(kaiserHBFilter);
legend(fvt,"DesignMethod = kaiser");

"equiripple" と "kaiser" に加えて、DesignMethod を "Auto" に設定することもできます。DesignMethod を "Auto" に設定すると、関数 designHalfbandFIR はフィルター設計パラメーターに基づいて設計法を自動的に選択します。これは、DesignMethod パラメーターを指定しない場合に使用される既定の設計法です。Verbose=true の引数を渡すことで、どの設計法が選択されているかを特定できます。以下のコードでは、関数は等リップル設計を自動的に選択します。

TW = 1/44.1; 
Ast = 180; 
autoHBFilter = designHalfbandFIR(DesignMethod="auto",...
                        TransitionWidth=TW,...
                        StopbandAttenuation=Ast, ...
                        Verbose=true);

designHalfbandFIR(TransitionWidth=0.022675736961451247, StopbandAttenuation=180, DesignMethod="kaiser", Passband="lowpass", Structure="single-rate", SystemObject=false)

fvt = fvtool(autoHBFilter);
legend(fvt,"designHalfbandFIR, DesignMethod = auto");

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 an object of type line. This object represents designHalfbandFIR, DesignMethod = auto.

dsp.FIRHalfbandDecimator System object および dsp.FIRHalfbandInterpolator System object は、"auto" 設計法もサポートしています。

Fs = 44.1e3;
TW = 1000; % Transition width
Ast = 60;  % 60 dB minimum attenuation in the stopband
autoHBFilter = dsp.FIRHalfbandDecimator(Specification="Transition width and stopband attenuation",...
    SampleRate=Fs,...
    TransitionWidth=TW,...
    StopbandAttenuation=Ast, ...
    DesignMethod="auto");
fvt = fvtool(autoHBFilter);
legend(fvt,"dsp.FIRHalfbandDecimator, DesignMethod = auto");

Figure Figure 6: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (kHz), ylabel Magnitude (dB) contains 2 objects of type line. This object represents dsp.FIRHalfbandDecimator, DesignMethod = auto.