Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

指定群遅延が与えられた IIR フィルター設計

この例では、fdesign.arbgrpdelay フィルター デザイナーを使用して、任意の群遅延フィルターを設計する方法を示します。このデザイナーでは、最小 p 次制約付き最適化アルゴリズムを使用して、指定済みの群遅延を満たすオールパス IIR フィルターを設計します。

fdesign.arbgrpdelay は、群遅延イコライズに使用できます。

任意の群遅延フィルター デザイナー

design.arbgrpdelay を使用すると、希望する群遅延応答を指定してオールパス フィルターを設定できます。目的の群遅延は、相対的な意味で指定されます。実際の群遅延は、フィルター次数によって変わります (次数が高いほど遅延が大きい)。ただし、フィルター次数に起因して群遅延のオフセットを減算した場合、設計されたフィルターの群遅延は、目的の群遅延に一致する傾向があります。2 つの異なるフィルター次数を使ったコードの例を、以下に示します。

N = 8;         % Filter order
N2 = 10;       % Alternate filter order
F = [0 0.1 1]; % Frequency vector
Gd = [2 3 1];  % Desired group delay
R = 0.99;      % Pole-radius constraint

オールパス フィルターでは、分子が必ず対応する分母の反転となります。このため、fdesign.arbgrpdelay では、異なる分母と分子の次数を指定できません。

以下のコードは、指定の周波数ポイント F で群遅延値 Gd をもつ任意のシングル バンド群遅延設計を示します。シングル バンドの設計では、ナイキスト間隔全体 [0 1]*pi rad/sample に対応する周波数値で群遅延を指定します。

arbGrpSpec = fdesign.arbgrpdelay('N,F,Gd',N,F,Gd)
arbGrpSpec = 
  arbgrpdelay with properties:

               Response: 'Arbitrary Group Delay'
          Specification: 'N,F,Gd'
            Description: {3x1 cell}
    NormalizedFrequency: 1
            FilterOrder: 8
            Frequencies: [0 0.1000 1]
             GroupDelay: [2 3 1]

arbGrpDelFilter1 = design(arbGrpSpec,'MaxPoleRadius',R, ...
    'SystemObject', true)
arbGrpDelFilter1 = 
  dsp.SOSFilter with properties:

            Structure: 'Direct form II'
    CoefficientSource: 'Property'
            Numerator: [4x3 double]
          Denominator: [4x3 double]
       HasScaleValues: true
          ScaleValues: [0.2506 1 1 1 1.4861]

  Use get to show all properties

0 から 1 までの一連の周波数点において群遅延の合計を測定します。測定法を使用して、フィルターのノミナル群遅延を測定します。

Fpoints = 0:0.001:1;
M1 = measure(arbGrpSpec,arbGrpDelFilter1,Fpoints);

N2 に等しい次数をもつ別のフィルターを設計します。

arbGrpSpec.FilterOrder = N2;
arbGrpDelFilter2 = design(arbGrpSpec,'MaxPoleRadius',R, ...
    'SystemObject', true);
M2 = measure(arbGrpSpec,arbGrpDelFilter2,Fpoints);

測定した群遅延の合計からノミナル群遅延を引いたものをプロットします。

plot(Fpoints, M1.TotalGroupDelay-M1.NomGrpDelay, 'b',...
     Fpoints, M2.TotalGroupDelay-M2.NomGrpDelay, 'g',...
     [0 0.1 1], [2 3 1], 'r');
xlabel('Normalized Frequency (\times\pi rad/sample)');
ylabel('Group delay (samples)'); grid on;
legend('8th order design','10th order design','desired response')

Figure contains an axes object. The axes object with xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (samples) contains 3 objects of type line. These objects represent 8th order design, 10th order design, desired response.

以下のプロットに見られるように、2 つの構成の実際の群遅延は異なります。つまり、この結果では、1 つの設計は目的とする相対群遅延へのより近い近似 (より少ないリップル) とフィルター内での規模の大きい全体的な遅延の間でのバランスを探さなければならないことを意味しています。

fvt = fvtool(arbGrpDelFilter1,arbGrpDelFilter2,'Analysis', 'grpdelay');
legend(fvt, '8th order design','10th order design')

Figure Figure 1: Group delay contains an axes object. The axes object with title Group delay, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Group delay (in samples) contains 2 objects of type line. These objects represent 8th order design, 10th order design.

通過帯域群遅延イコライズ

fdesign.arbgrpdelay は主に、IIR フィルターの非線形位相応答を補正するために使用します。補正フィルターが オールパスであるため、振幅応答に影響を与えることなく、補正する対象のフィルターでカスケードを作成できます。2 つのフィルターのカスケードは IIR フィルターとなるため、線形位相はありません (安定している間)。ただし、フィルター全体の通過帯域内では、ほぼ線形位相応答である応答を得ることができます。

ローパス イコライズ

以下の例では、fdesign.arbgrpdelay を使用して、振幅応答に影響を及ぼすことなく、ローパス楕円フィルターの群遅延応答をイコライズしています。

他のすべての周波数帯域の群遅延を指定せずに (「don't care」領域)、1 つ以上の対象帯域で目的の群遅延値を指定するには、マルチバンド設計を使用します。この例では、対象となるのは 1 つの帯域だけであり、ローパス フィルターの通過帯域に相当します。この帯域で群遅延を補正する必要がありますが、結果として群遅延値が帯域から外れていてもかまいません。

通過帯域周波数が 0.2π ラジアン/サンプルの楕円フィルターを設計します。通過帯域全体で群遅延の合計を測定します。

ellipFilter = design(fdesign.lowpass('N,Fp,Ap,Ast',4,0.2,1,40),...
    'ellip', 'SystemObject', true);
wncomp = 0:0.001:0.2; 
g = grpdelay(ellipFilter,wncomp,2); % samples
g1 = max(g)-g;

任意の 8 次群遅延オールパス フィルターを設計します。マルチバンド設計を使用して、シングル バンドを指定します。

allpassSpec = fdesign.arbgrpdelay('N,B,F,Gd',8,1,wncomp,g1)
allpassSpec = 
  arbgrpdelay with properties:

               Response: 'Arbitrary Group Delay'
          Specification: 'N,B,F,Gd'
            Description: {4x1 cell}
    NormalizedFrequency: 1
            FilterOrder: 8
                 NBands: 1
          B1Frequencies: [0 1.0000e-03 0.0020 0.0030 0.0040 0.0050 0.0060 0.0070 0.0080 0.0090 0.0100 0.0110 0.0120 0.0130 0.0140 0.0150 0.0160 0.0170 0.0180 0.0190 0.0200 0.0210 0.0220 0.0230 0.0240 0.0250 0.0260 0.0270 0.0280 0.0290 ... ] (1x201 double)
           B1GroupDelay: [14.5838 14.5835 14.5824 14.5808 14.5784 14.5753 14.5716 14.5672 14.5622 14.5564 14.5500 14.5429 14.5352 14.5267 14.5176 14.5078 14.4974 14.4862 14.4744 14.4619 14.4488 14.4350 14.4205 14.4054 14.3896 14.3731 ... ] (1x201 double)

allpassFilter = design(allpassSpec,'iirlpnorm', 'SystemObject', true)
allpassFilter = 
  dsp.SOSFilter with properties:

            Structure: 'Direct form II'
    CoefficientSource: 'Property'
            Numerator: [4x3 double]
          Denominator: [4x3 double]
       HasScaleValues: true
          ScaleValues: [0.1653 1 1 1 16.4771]

  Use get to show all properties

補正フィルターを使用して元のフィルターをカスケードし、目標の群遅延イコライズを達成します。ホワイト ノイズを処理し、2 つの出力段階で群遅延を推定して、検証を行います。

samplesPerFrame = 2048;
wn = (2/samplesPerFrame) * (0:samplesPerFrame-1);
numRealPoints = samplesPerFrame/2 + 1;
tfEstimator = dsp.TransferFunctionEstimator('FrequencyRange','onesided',...
      'SpectralAverages',64);
scope = dsp.ArrayPlot('PlotType','Line','YLimits',[0 40],...
      'YLabel','Group Delay (samples)',...
      'XLabel','Normalized Frequency (x pi rad/sample)',...
      'SampleIncrement',2/samplesPerFrame,...
      'Title',['Original (1), Compensated (2), ',...
      'Expected Compensated (3)'], 'ShowLegend', true);

gdOrig = grpdelay(ellipFilter, numRealPoints);
gdComp = grpdelay(allpassFilter, numRealPoints);
range = wn < wncomp(end);
gdExp = nan(numRealPoints, 1); 
gdExp(range) = gdOrig(range) + gdComp(range);

フィルター カスケードにランダムなサンプルをストリーミングします。

Nframes = 300;
for k = 1:Nframes
    x = randn(samplesPerFrame,1);  % Input signal = white Gaussian noise
    
    y_orig = ellipFilter(x);       % Filter noise with original IIR filter
    y_corr = allpassFilter(y_orig);% Compensating filter
    
    Txy = tfEstimator([x, x],[y_orig, y_corr]);
    gdMeas = HelperMeasureGroupDelay(Txy, [], 20);
    scope([gdMeas, gdExp]);
end

バンドパス イコライズ

[0.3 0.4]*pi rad/sample 間隔で、通過帯域領域を含むバンドパス チェビシェフ フィルターの通過帯域群遅延イコライザーを設計します。前述の例と同様に、対象となるのは 1 つの帯域だけであり、フィルターの通過帯域に対応します。この帯域で群遅延を補正する必要がありますが、結果として群遅延値が帯域から外れていてもかまわないため、マルチバンド設計を使用して、シングル バンドを指定します。

バンドパス チェビシェフ タイプ 1 フィルターを設計して、通過帯域全体でフィルターの群遅延の合計を測定します。

bandpassFilter = design(fdesign.bandpass('N,Fp1,Fp2,Ap',4,0.3,0.4,1), ...
    'cheby1', 'SystemObject', true);
wncomp = 0.3:0.001:0.4;
g = grpdelay(bandpassFilter,wncomp,2);
g1 = max(g)-g;

任意の 8 次群遅延フィルターを設計します。極半径は 0.95 を超えないように制限されています。

allpassSpec = fdesign.arbgrpdelay('N,B,F,Gd',8,1,wncomp,g1);
allpassFilter = design(allpassSpec,'iirlpnorm','MaxPoleRadius',0.95, ...
    'SystemObject', true);

補正フィルターを使用して元のフィルターをカスケードし、目標の群遅延イコライズを達成します。ホワイト ノイズを処理し、2 つの出力段階で群遅延を推定して、検証を行います。

gdOrig = grpdelay(bandpassFilter, numRealPoints);
gdComp = grpdelay(allpassFilter, numRealPoints);
range = wn > wncomp(1) & wn < wncomp(end);
gdExp = nan(numRealPoints,1); gdExp(range) = gdOrig(range) + gdComp(range);

release(scope), scope.YLimits = [0 55];
release(tfEstimator);

フィルター カスケードにランダムなサンプルをストリーミングします。

for k = 1:Nframes
    x = randn(samplesPerFrame,1);  % Input signal = white Gaussian noise
    
    y_orig = bandpassFilter(x);    % Filter noise with original IIR filter
    y_corr = allpassFilter(y_orig);% Compensating filter
    
    Txy = tfEstimator([x, x],[y_orig, y_corr]);
    gdMeas = HelperMeasureGroupDelay(Txy, [], 20);
    scope([gdMeas, gdExp]);
end

結果として得られるフィルターには制約された 1 組の極があります。通過帯域 ([0.3 0.4]*pi rad/sample) の群遅延誤差は、0.2 サンプル未満です。

バンドストップ イコライズ

[1 0.4]*pi rad/sample 間隔で、サンプル周波数 1KHz を含むバンドストップ チェビシェフ フィルターの通過帯域群遅延イコライザーを設計します。バンドストップ フィルターには、[0 150] Hz および [200 500] Hz の間隔で 2 つの通過帯域領域があります。これらの帯域で群遅延を補正する必要があるため、マルチバンド設計を使用して、2 つの帯域 (バンド) を指定します。

バンドストップ チェビシェフ タイプ 2 フィルターを設計して、通過帯域全体でフィルターの群遅延の合計を測定します。fdesign.arbgrpdelay はサンプリング周波数を秒単位で指定する群遅延仕様を前提としているので、測定された群遅延を秒に変換します。

Fs = 1e3;
bandstopFilter = design(fdesign.bandstop('N,Fst1,Fst2,Ast',6,150,400,1,...
    Fs), 'cheby2', 'SystemObject', true);
f1 = 0.0:0.5:150; % Hz
g1 = grpdelay(bandstopFilter,f1,Fs).'/Fs; % seconds
f2 = 400:0.5:Fs/2; % Hz
g2 = grpdelay(bandstopFilter,f2,Fs).'/Fs; % seconds
maxg = max([g1 g2]);

任意の 14 次群遅延オールパス フィルターを設計します。極半径は 0.95 を超えないように制限されています。群遅延仕様は秒単位で指定され、周波数仕様はヘルツ単位で指定されています。

allpassSpec = fdesign.arbgrpdelay('N,B,F,Gd',14,2,f1,maxg-g1,f2,...
    maxg-g2,Fs);
allpassFilter = design(allpassSpec,'iirlpnorm','MaxPoleRadius',0.95,...
    'SystemObject', true);

元のフィルターを補正フィルターとカスケード接続して、ホワイト ノイズを処理し、2 つの出力段階で群遅延を推定します。

gdOrig = grpdelay(bandstopFilter, numRealPoints);
gdComp = grpdelay(allpassFilter, numRealPoints);
fcomp = (Fs/samplesPerFrame) * (0:samplesPerFrame-1);
range = (fcomp>f1(1) & fcomp<f1(end)) | (fcomp>f2(1) & fcomp<f2(end));
gdExp = nan(numRealPoints,1); gdExp(range) = gdOrig(range) + gdComp(range);

release(scope), 
    scope.YLimits = [0 40];
    scope.SampleIncrement = Fs/samplesPerFrame;
    scope.YLabel = 'Group Delay (samples)';
    scope.XLabel = 'Frequency (Hz)';
release(tfEstimator)

フィルター カスケードにランダムなサンプルをストリーミングします。

Nframes = 300;
for k = 1:Nframes
    x = randn(samplesPerFrame,1);  % Input signal = white Gaussian noise
    
    y_orig = bandstopFilter(x);    % Filter noise with original IIR filter
    y_corr = allpassFilter(y_orig);% Compensating filter
    
    Txy = tfEstimator([x, x],[y_orig, y_corr]);
    gdMeas = HelperMeasureGroupDelay(Txy, [], 12);
    scope([gdMeas, gdExp]);
end

結果として得られるフィルターには制約された 1 組の極があります。通過帯域の群遅延誤差は、3 サンプル未満となります。