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

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

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

• 通過帯域リップルと阻止帯域リップルは同じでなければなりません。

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

ハーフバンド係数の取得

関数 firhalfband は、FIR ハーフバンド等リップル フィルターの係数を返します。簡単な例として、96 kHz でサンプリングされたデータで処理を行い、通過帯域周波数 22 kHz をもつハーフバンド フィルターを考えます。

Fs  = 96e3;
Fp  = 22e3;
N   = 100;
num = firhalfband(N,Fp/(Fs/2));
zerophase(num,1,linspace(0,Fs/2,512),Fs);

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

fvt = fvtool(num,Fs=Fs);
fvt.Analysis = "impulse";

dsp.FIRHalfbandInterpolator および dsp.FIRHalfbandDecimator

関数 firhalfband には、この他にいくつかの設計オプションがあります。ただし、ストリーミング データを処理する場合はdsp.FIRHalfbandInterpolator System object およびdsp.FIRHalfbandDecimator System object を使用することを推奨します。これら 2 つの System object では、係数を設計するだけでなく、ポリフェーズを効率的に実装できます。これらは倍精度や単精度の浮動小数点データだけでなく、固定小数点データのフィルター処理もサポートします。また、C コードおよび HDL コード生成に加えて、ARM® Cortex® M および ARM® Cortex® A に最適化されたコード生成もサポートします。

halfbandInterpolator = dsp.FIRHalfbandInterpolator(SampleRate=Fs,...
    Specification="Filter order and transition width",...
    FilterOrder=N,TransitionWidth=4000);
fvt = fvtool(halfbandInterpolator,Fs=2*Fs); %#ok<NASGU> 

内挿を実行するには、dsp.FIRHalfbandInterpolator System object™ を使用します。これはマルチレート フィルターであるため、サンプルレートが意味するものを定義することが重要です。この System object やその他すべての System object では、サンプルレートは入力信号のサンプルレートを指します。ただし、FVTool では、サンプルレートはフィルターが処理を行うレートとして定義されます。内挿の場合には、アップサンプリング後にフィルター処理を行います (概念上)。したがって、2 でアップサンプリングされるため、FVTool のサンプル レートを $2\mathrm{Fs}$ として指定する必要があります。

FrameSize = 256;
scope = spectrumAnalyzer(SampleRate=2*Fs);
sine1 = dsp.SineWave(Frequency=10e3,SampleRate=Fs,...
    SamplesPerFrame=FrameSize);
sine2 = dsp.SineWave(Frequency=20e3,SampleRate=Fs,...
    SamplesPerFrame=FrameSize);
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 減衰することに注意してください。これはハーフバンド フィルターによって生じるおおよその減衰量です。フィルターの群遅延を補正することで、入力サンプルに内挿されたサンプルを重ねたプロットを取得できます。フィルターの出力時に入力サンプルが変わっていないことに注意してください。これは、ハーフバンド フィルターのポリフェーズ分岐のいずれかが、入力サンプルを変化させない純粋な遅延分岐であるためです。

grpDel = 50;
n = 0:2:511;
stem(n(1:end-grpDel/2),x(1:end-grpDel/2),"k","filled")
hold on
nu = 0:511;
stem(nu(1:end-grpDel),y(grpDel+1:end))
legend("Input samples","Interpolated samples")

間引きの場合には、フィルター処理後にダウンサンプリング (概念上) するため、dsp.FIRHalfbandDecimator で指定したサンプル レートはフィルターのサンプル レートに対応します。そのため、間引きでは、FVTool で指定した $\mathrm{Fs}$ を係数で乗算する必要はありません。

FrameSize = 256;
FsIn = 2*Fs;
halfbandDecimator = dsp.FIRHalfbandDecimator(SampleRate=FsIn,...
    Specification="Filter order and transition width",...
    FilterOrder=N,TransitionWidth=4000);
fvt = fvtool(halfbandDecimator,Fs=FsIn);%#ok<NASGU> 

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

release(scope);

フィルター係数の取得

関数 tf を使用すると、内挿や間引きからフィルター係数を抽出できます。

num = tf(halfbandInterpolator); % Or num = tf(halfbandDecimator);

さまざまな設計仕様の使用

フィルター次数と遷移幅を指定する代わりに、所定の遷移幅と阻止帯域の減衰量を与えて最小次数のフィルターを設計できます。

Ast = 80; % 80 dB
halfbandInterpolator = dsp.FIRHalfbandInterpolator(SampleRate=Fs,...
    Specification="Transition width and stopband attenuation",...
    StopbandAttenuation=Ast,TransitionWidth=4000);
fvtool(halfbandInterpolator,Fs=2*Fs);

すべての内挿と同様に、絶対単位での通過帯域ゲインは内挿係数と等価です (ハーフバンドの場合は 2)。これは、6.02 dB の通過帯域ゲインに相当します。

フィルター次数と阻止帯域の減衰量を指定することも可能です。

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

内挿と異なり、間引きでは通過帯域に 1 (0 dB) のゲインがあります。

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

ハーフバンドの内挿と間引きを使用すると、合成/解析フィルター バンクを効率的に実装できます。これまでに示したハーフバンド フィルターは、すべてローパス フィルターでした。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"}); %#ok<CLARRSTR> 

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

詳細設計オプション: さまざまな設計アルゴリズムの指定

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

Fs = 44.1e3;
N  = 90;
TW = 1000;
equirippleHBFilter = dsp.FIRHalfbandInterpolator(DesignMethod="Equiripple",...
    Specification="Filter order and transition width",...
    SampleRate=Fs,...
    TransitionWidth=TW,...
    FilterOrder=N); 
kaiserHBFilter = dsp.FIRHalfbandInterpolator(DesignMethod="Kaiser",...
    Specification="Filter order and transition width",...
    SampleRate=Fs,...
    TransitionWidth=TW,...
    FilterOrder=N); 

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

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

fdesign.interpolatorオブジェクトおよびfdesign.decimatorオブジェクトを使用すると、最小二乗線形フィルター、FIR フィルター設計などのその他の設計アルゴリズムを利用できます。任意のフィルター仕様オブジェクトに対して使用できる設計法の一覧を確認するには、関数 designmethods を使用します。

filtSpecs = fdesign.interpolator(2,"halfband","N,TW",N,TW/Fs);
designmethods(filtSpecs,"FIR");

FIR Design Methods for class fdesign.interpolator (N,TW):


equiripple
firls
kaiserwin

阻止帯域の減衰の制御

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

Ast  = 60; % Minimum stopband attenuation
equirippleHBFilter = dsp.FIRHalfbandInterpolator(DesignMethod="Equiripple",...
    Specification="Filter order and stopband attenuation",...
    SampleRate=Fs,...
    StopbandAttenuation=Ast,...
    FilterOrder=N); 
kaiserHBFilter = dsp.FIRHalfbandInterpolator(DesignMethod="Kaiser",...
    Specification="Filter order and stopband attenuation",...
    SampleRate=Fs,...
    StopbandAttenuation=Ast,...
    FilterOrder=N); 
fvt = fvtool(equirippleHBFilter,kaiserHBFilter,Fs=2*Fs);
legend(fvt,"Equiripple design","Kaiser-window design")

最小次数設計

設計仕様を満たすために必要な最小次数のフィルターを設計する場合には、等リップル設計に加えてカイザー ウィンドウ設計も使用できます。カイザー ウィンドウ設計のための実際の次数は、等リップル設計に必要な次数より大きくなりますが、代わりに全体的な阻止帯域の減衰量は改善されます。

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

equirippleHBFilter = dsp.FIRHalfbandDecimator(DesignMethod="Equiripple",...
    Specification="Transition width and stopband attenuation",...
    SampleRate=Fs,...
    TransitionWidth=TW,...
    StopbandAttenuation=Ast); 
kaiserHBFilter = dsp.FIRHalfbandDecimator(DesignMethod="Kaiser",...
    Specification="Transition width and stopband attenuation",...
    SampleRate=Fs,...
    TransitionWidth=TW,...
    StopbandAttenuation=Ast); 
fvt = fvtool(equirippleHBFilter,kaiserHBFilter);
legend(fvt,"Equiripple design","Kaiser-window design")

フィルター設計法の自動選択

"Equiripple" および "Kaiser" に加えて、dsp.FIRHalfbandDecimator System object および dsp.FIRHalfbandInterpolator System object の DesignMethod プロパティを "Auto" に指定することもできます。DesignMethod を "Auto" に設定した場合は、フィルター設計パラメーターに基づき、オブジェクトによって自動的にフィルター設計法が選択されます。

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

上記のフィルター仕様の場合、System object によって等リップル フィルターが設計されることを振幅応答から観察できます。阻止帯域の減衰量が非常に大きい、または遷移幅が非常に狭いなど、設計の制約がかなり厳しい場合は、アルゴリズムによってカイザー ウィンドウ法が自動的に選択されます。カイザー ウィンドウ法は、非常に厳しい仕様のフィルターを設計する場合に最適です。ただし、設計の制約が厳しくない場合は、アルゴリズムによって等リップル設計が実行されます。

以下は、フィルター仕様が厳しすぎるため、等リップル設計を実行できない場合を示します。オブジェクトの DesignMethod プロパティは "Equiripple" に設定されます。これにより、このオブジェクトは等リップル特性を使用してフィルターを設計しようとするため、フィルターの周波数応答からわかるように、設計が収束しません。

Fs = 192e3;
TW = 100; % Transition width
Ast = 180;  % 180 dB minimum attenuation in the stopband
equirippleHBFilter = dsp.FIRHalfbandDecimator(DesignMethod="Equiripple",...
    TransitionWidth=TW,...
    StopbandAttenuation=Ast,...
    SampleRate=Fs);
fvt = fvtool(equirippleHBFilter);

Warning: Final filter order of 10448 is probably too high to optimally meet the constraints.

Warning: Design is not converging.  Number of iterations was 5
1) Check the resulting filter using freqz.
2) Check the specifications.
3) Filter order may be too large or too small.
4) For multiband filters, try making the transition regions more similar in width.
If err is very small, filter order may be too high

legend(fvt,"DesignMethod = Equiripple");

この場合、より適切に収束するフィルターを設計するには、DesignMethod プロパティを "Auto" または "Kaiser" に設定します。このオブジェクトはカイザー ウィンドウ法を使用するハーフバンド フィルターを設計します。この場合、設計されたフィルターは、厳しい設計仕様を正確に満たしていないことに注意してください。しかし、2 つのフィルターの周波数応答を比較するとわかるように、フィルターはこちらの方が適切に収束します。

Fs = 192e3;
TW = 100; % Transition width
Ast = 180;  % 180 dB minimum attenuation in the stopband
autoHBFilter = dsp.FIRHalfbandDecimator(DesignMethod="Auto",...
    TransitionWidth=TW,...
    StopbandAttenuation=Ast,...
    SampleRate=Fs);
fvt = fvtool(autoHBFilter);

Warning: The designed filter does not meet the filter specifications. Increase the Transition Width or reduce the Stopband Attenuation.

legend(fvt,"DesignMethod = Auto");

阻止帯域の減衰が増大する等リップル設計

fdesign.interpolator オブジェクトおよび fdesign.decimator オブジェクトを使用し、関数 design の "StopbandShape" 引数をオプションで指定することで、等リップル設計の阻止帯域の形状を変更することもできます。

Fs = 44.1e3;
TW = 1000/(Fs/2); % Transition width
Ast = 60;  % 60 dB minimum attenuation in the stopband
filtSpecs = fdesign.decimator(2,"halfband","TW,Ast",TW,Ast);
equirippleHBFilter1 = design(filtSpecs,"equiripple",...
    StopbandShape="1/f",StopbandDecay=4,SystemObject=true);
equirippleHBFilter2 = design(filtSpecs,"equiripple",...
    StopbandShape="linear",StopbandDecay=53.333,SystemObject=true);
fvt = fvtool(equirippleHBFilter1,equirippleHBFilter2,...
    Fs=Fs);
legend(fvt,"Stopband decaying as (1/f)^4","Stopband decaying linearly")

ハイパス ハーフバンド フィルター

ハイパス ハーフバンド フィルターは、ローパス ハーフバンド フィルターの 1 つおきの係数の符号を変更することで、ローパス ハーフバンド フィルターから取得できます。あるいは、fdesign.decimator オブジェクトの Type プロパティを "Highpass" に設定することで、ハイパス ハーフバンドを直接設計することもできます。

filtSpecs = fdesign.decimator(2,"halfband",...
    "TW,Ast",TW,Ast,Type="Highpass");
halfbandHPFilter = design(filtSpecs,"equiripple",...
    StopbandShape="linear",StopbandDecay=53.333,SystemObject=true);
fvt = fvtool(halfbandHPFilter,equirippleHBFilter2,Fs=Fs);
legend(fvt,"Highpass halfband filter","Lowpass halfband filter")

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい