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

スクリプトを開く

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

ハーフバンドフィルターには 2 つの重要な特徴があります。つまり、通過帯域リップルと阻止帯域リップルは同じでなければならず、通過帯域エッジと阻止帯域エッジの周波数はハーフバンド周波数 Fs/4 (または正規化周波数の pi/2 ラジアン/サンプル) から等距離になければなりません。

ハーフバンド係数の取得

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

Fs  = 96e3;
Fp  = 22e3;
N   = 100;
num = firhalfband(N,Fp/(Fs/2));
fvt = fvtool(num,'Fs',Fs,'Color','white');
fvt.MagnitudeDisplay = 'Zero-phase';

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

fvt.Analysis = 'impulse';

`dsp.FIRHalfbandInterpolator` および `dsp.FIRHalfbandDecimator`

関数 firhalfband には、この他にいくつかの設計オプションがあります。ただし、ほとんどの場合、dsp.FIRHalfbandInterpolator および dsp.FIRHalfbandDecimator で直接操作する方が適切です。これら 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);
fvtool(halfbandInterpolator,'Fs',2*Fs,'Color','white');

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

FrameSize = 256;
scope = dsp.SpectrumAnalyzer('SampleRate',2*Fs,'SpectralAverages',5);
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 で指定した Fs を係数で乗算する必要はありません。

FrameSize = 256;
FsIn = 2*Fs;
halfbandDecimator = dsp.FIRHalfbandDecimator('SampleRate',FsIn, ...
    'Specification','Filter order and transition width', ...
    'FilterOrder',N,'TransitionWidth',4000);
fvtool(halfbandDecimator,'Fs',FsIn,'Color','white');
scope = dsp.SpectrumAnalyzer('SampleRate',Fs,'SpectralAverages',5);
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,'Color','white');

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

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

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

内挿と異なり、間引きでは通過帯域に 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 = dsp.SpectrumAnalyzer('SampleRate',Fs1, ...
    'PlotAsTwoSidedSpectrum',false,'ShowLegend',true,'YLimits', ...
    [-120 30], ...
    'Title', ...
    'Input Signal and Output Signal of Quadrature Mirror Filter');
scope.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);

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

これまでに示した設計は、すべて最適な等リップル設計でした。fdesign.interpolator および fdesign.decimator を使用すると、その他の設計アルゴリズムを利用できます。

Fs = 44.1e3;
N  = 90;
TW = 1000/Fs; % Transition width
filtSpecs = fdesign.interpolator(2,'halfband','N,TW',N,TW);
equirippleHBFilter = design(filtSpecs,'equiripple','SystemObject',true);
leastSquaresHBFilter = design(filtSpecs,'firls','SystemObject',true);
kaiserHBFilter = design(filtSpecs,'kaiserwin','SystemObject',true);

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

fvt = fvtool(equirippleHBFilter,leastSquaresHBFilter,kaiserHBFilter, ...
    'Fs',2*Fs,'Color','white');
legend(fvt,'Equiripple design','Least-squares design', ...
    'Kaiser-window design')

阻止帯域の減衰の制御

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

Ast  = 60; % Minimum stopband attenuation
filtSpecs = fdesign.interpolator(2,'halfband','N,Ast',N,Ast);
equirippleHBFilter = design(filtSpecs,'equiripple','SystemObject',true);
kaiserHBFilter = design(filtSpecs,'kaiserwin','SystemObject',true);
fvt = fvtool(equirippleHBFilter,kaiserHBFilter,'Fs',2*Fs,'Color','white');
legend(fvt,'Equiripple design','Kaiser-window design')

最小次数設計

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

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);
equirippleHBFilter = design(filtSpecs,'equiripple','SystemObject',true);
kaiserHBFilter = design(filtSpecs,'kaiserwin','SystemObject',true);
fvt = fvtool(equirippleHBFilter,kaiserHBFilter,'Fs',Fs,'Color','white');
legend(fvt,'Equiripple Design','Kaiser-window design')

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

カイザーウィンドウフィルターを設計する代わりに、変更した '等リップル' 設計を使用して阻止帯域の減衰量を増やすこともできます。

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,'Color','white');
legend(fvt,'Stopband decaying as (1/f)^4','Stopband decaying linearly')

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

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

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