遷移帯域が狭い FIR フィルターの効率的な設計

遷移帯域が狭いローパス フィルターの設計

FIR フィルターを使用する欠点の 1 つに、フィルターの遷移帯域に反比例してフィルター次数が増加する傾向にあることが挙げられます。以下の設計仕様 (リップルが線形単位で指定されている場合) を考えます。

Fpass = 2866.5;  % Passband edge
Fstop = 3087;    % Stopband edge
Apass = 0.0174;  % Passband ripple, 0.0174 dB peak to peak
Astop = 66.0206; % Stopband ripple, 66.0206 dB of minimum attenuation
Fs    = 44.1e3;

lowpassSpec = fdesign.lowpass(Fpass,Fstop,Apass,Astop,Fs);

この仕様を満たす通常の線形位相等リップル設計は、以下のように設計できます。

eqripFilter = design(lowpassSpec,'equiripple','SystemObject',true);
cost(eqripFilter)

ans = struct with fields:
                  NumCoefficients: 695
                        NumStates: 694
    MultiplicationsPerInputSample: 695
          AdditionsPerInputSample: 694

必要なフィルター長は 694 タップになります。

内挿 FIR (IFIR) の設計

IFIR 設計アルゴリズムを使用すると、必要な乗数の合計数が減少するため、上述の仕様を効率的に設計することができます。これを行うには、設計問題を 2 つの段に分割します。1 つは多数の乗数を使用せずに厳格な仕様を満たすためにアップサンプルしたフィルター、もう 1 つは前のフィルターのアップサンプル時に作成されたイメージを削除するフィルターです。

interpFilter = design(lowpassSpec,'ifir','SystemObject',true);

明らかに状況は悪くなっています。ここでは、694 乗数のシングル フィルターの代わりに、合計 804 乗数の 2 つのフィルターがあり、乗数が明らかに以前よりも増えています。しかし、2 段目の詳細を見ると、約 5 分の 1 の乗数のみが非ゼロであることがわかります。これにより乗数の実際の合計数が 694 から 208 に削減されます。

cost(interpFilter)

ans = struct with fields:
                  NumCoefficients: 208
                        NumStates: 802
    MultiplicationsPerInputSample: 208
          AdditionsPerInputSample: 206

2 つの設計の応答を比較します。

fvt = fvtool(eqripFilter,interpFilter,'Color','White');
legend(fvt,'Equiripple design', 'IFIR design','Location','Best')

アップサンプリング係数の手動制御

前の例では、乗数の合計数が最小化されるように使用するアップサンプリング係数を自動的に決定しました。この仕様では、最適なアップサンプリング係数は 5 でした。一方、以下の設計オプションを考えてください。

opts = designopts(lowpassSpec,'ifir')

opts = struct with fields:
      FilterStructure: 'dffir'
     UpsamplingFactor: 'auto'
    JointOptimization: 0
         SystemObject: 0

この場合は、アップサンプリング係数を制御できることがわかります。たとえば、5 ではなく 4 でアップサンプリングを行う場合は、以下のようになります。

opts.UpsamplingFactor = 4;
opts.SystemObject = true; 
upfilter = design(lowpassSpec,'ifir',opts);
cost(upfilter)

ans = struct with fields:
                  NumCoefficients: 217
                        NumStates: 767
    MultiplicationsPerInputSample: 217
          AdditionsPerInputSample: 215

これにより合計 217 の非ゼロ乗数の設計が得られます。

同時最適化の使用

IFIR に使用する 2 つのフィルターを同時に設計することができます。この場合は、設計時間がかかるものの乗数の数を大幅に削減することができます (アルゴリズムの性質から、一部のケースでは設計が収束しない可能性があります)。

opts.UpsamplingFactor = 'auto'; % Automatically determine the best factor
opts.JointOptimization = true;
ifirFilter = design(lowpassSpec,'ifir',opts);
cost(ifirFilter)

ans = struct with fields:
                  NumCoefficients: 172
                        NumStates: 726
    MultiplicationsPerInputSample: 172
          AdditionsPerInputSample: 170

この設計の最適なアップサンプリング係数は 6 でした。非ゼロ乗数の数はわずか 152 となります。

マルチレート/多段手法を使用した設計の効率化

これまでに検討した設計では、シングルレート手法を使用しました。これは、MPIS (入力サンプルあたりに必要な乗算) の数が非ゼロ乗数の数と等しくなることを意味します。たとえば、最後に示した設計には 152 MPIS が必要になります。最初の単一の等リップル設計には、694 MPIS が必要でした。

間引きと内挿を組み合わせたマルチレート/多段手法を使用することにより、MPIS 数の少ない効率的な設計を得ることができます。間引きでは、MPIS 数 (平均) は、乗数の数を間引き係数で除算することにより求めることができます。

multistageFilter = design(lowpassSpec,'multistage','SystemObject',true)

multistageFilter = 
  dsp.FilterCascade with properties:

         Stage1: [1x1 dsp.FIRDecimator]
         Stage2: [1x1 dsp.FIRDecimator]
         Stage3: [1x1 dsp.FIRDecimator]
         Stage4: [1x1 dsp.FIRInterpolator]
         Stage5: [1x1 dsp.FIRInterpolator]
         Stage6: [1x1 dsp.FIRInterpolator]
    CloneStages: true

cost(multistageFilter)

ans = struct with fields:
                  NumCoefficients: 396
                        NumStates: 352
    MultiplicationsPerInputSample: 73
          AdditionsPerInputSample: 70.8333

1 段目の乗数は 21 で、間引き係数は 3 です。そのため、MPIS 数は 7 です。2 段目は長さ 45、累積間引き係数は 6 です (つまり、この段の間引き係数を 1 段目の間引き係数で乗算したものです。これはこの段の入力サンプルが 1 段目の入力サンプルの 1/3 のレートで入るためです)。したがって、MPIS の平均数 (マルチレート/多段の全体フィルターの入力の参照) は 45/6=7.5 になります。最後に、3 段目の間引き係数が 1 であることを考えると、この段の MPIS の平均数は 130/6=21.667 になります。3 つの間引きの平均 MPIS の合計は 36.5 です。

内挿の場合のフィルターも、間引きと同じです。さらに、計算コストも同じです。したがって、マルチレート/多段設計全体の合計 MPIS 数は 73 になります。

次に、等リップル設計と以下の設計の応答を比較します。

fvt = fvtool(eqripFilter,multistageFilter,'Color','White', ...
    'NormalizeMagnitudeto1','on');
legend(fvt,'Equiripple design', 'Multirate/multistage design', ...
    'Location','NorthEast')

多段設計の阻止帯域の減衰が、他の設計の約 2 倍であることに注意してください。これは、必要な仕様を満たす上で障害となるエイリアシングを避けるために、間引きで必要とされる 66dB だけ帯域外コンポーネントを減衰させる必要があるためです。同様に、内挿でもこの設計の仕様を満たすために、イメージを 66dB 減衰させる必要があります。

段数の手動制御

ここで得られたマルチレート/多段設計は 6 段から構成されていました。段数は既定で自動的に決定されます。ただし、得られる段数を手動で制御することも可能です。次に例を示します。

lp4stage=design(lowpassSpec,'multistage','NStages',4,'SystemObject',true);
cost(lp4stage)

ans = struct with fields:
                  NumCoefficients: 516
                        NumStates: 402
    MultiplicationsPerInputSample: 86
          AdditionsPerInputSample: 84.5000

この場合の平均 MPIS 数は 86 です。

群遅延

各設計の群遅延を計算することができます。マルチレート/多段設計を使用した場合に、遅延が最も大きくなることに注意してください (これが計算コストの低い設計の欠点です)。IFIR 設計も、単一の等リップル設計より遅延が大きくなりますが、マルチレート/多段設計ほどではありません。

fvt = fvtool(eqripFilter,interpFilter,multistageFilter,...
    'Analysis','grpdelay','Color','White');
legend(fvt, 'Equiripple design','IFIR design',...
    'Multirate/multistage design');

信号のフィルター処理

次に例を使って、IFIR とマルチレート/多段設計では、計算量を大幅に抑えながら単一の等リップル設計と同じ性能が得られることを示します。カスケードに対してフィルター処理を実行するためには、generateFilteringCode(ifirFilter) および generateFilteringCode(multistageFilter) を呼び出す必要があります。ここでは、HelperIFIRFilter と HelperMultiFIRFilter を作成するためにこれは既に行われています。

scope = spectrumAnalyzer('SampleRate',Fs,...
    'PlotAsTwoSidedSpectrum',false,'YLimits', [-90 10],...
    'ShowLegend',true,'ChannelNames',{'Equiripple design',...
    'IFIR design','Multirate/multistage design'});
tic,
while toc < 20
    % Run for 20 seconds
    x = randn(6000,1);
    
    % Filter using single stage FIR filter
    y1 = eqripFilter(x);
    
    % Filter using IFIR filter
    y2 = HelperIFIRFilter(x);
    
    % Filter multi-stage/multi-rate FIR filters
    y3 = HelperMultiFIRFilter(x);
    
    % Compare the output from both approaches
    scope([y1,y2,y3])
end