ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

指定群遅延が与えられた 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 に対応する周波数値で群遅延を指定します。

h = fdesign.arbgrpdelay('N,F,Gd',N,F,Gd)
 
h =
 
               Response: 'Arbitrary Group Delay'
          Specification: 'N,F,Gd'               
            Description: {3x1 cell}             
    NormalizedFrequency: true                   
            FilterOrder: 8                      
            Frequencies: [0 0.1 1]              
             GroupDelay: [2 3 1]                
                                                
H1 = design(h,'MaxPoleRadius',R, 'SystemObject', true);

% Measure the total group delay at a set of frequency points from 0 to 1.
% Measure the nominal group delay of the filter using the measure method.
Fpoints = 0:0.001:1;
M1 = measure (h,H1,Fpoints);

% Design another filter with an order equal to N2.
h.FilterOrder = N2;
H2 = design(h,'MaxPoleRadius',R, 'SystemObject', true);
M2 = measure(h,H2,Fpoints);

% Plot the measured total group delay minus the nominal group delay.
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')

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

hFVT = fvtool(H1,H2,'Analysis', 'grpdelay');
legend(hFVT, '8th order design','10th order design')

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

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

ローパス イコライズ

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

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

% Design an elliptic filter with a passband frequency of 0.2*pi
% rad/sample. Measure the total group delay over the passband.
Hellip = design(fdesign.lowpass('N,Fp,Ap,Ast',4,0.2,1,40),...
    'ellip', 'SystemObject', true);
wncomp = 0:0.001:0.2;
g = grpdelay(Hellip,wncomp,2); % samples
g1 = max(g)-g;

% Design an 8th order arbitrary group delay allpass filter. Use a
% multiband design and specify a single band.
hgd = fdesign.arbgrpdelay('N,B,F,Gd',8,1,wncomp,g1)
Hgd = design(hgd,'iirlpnorm', 'SystemObject', true);
 
hgd =
 
               Response: 'Arbitrary Group Delay'
          Specification: 'N,B,F,Gd'             
            Description: {4x1 cell}             
    NormalizedFrequency: true                   
            FilterOrder: 8                      
                 NBands: 1                      
          B1Frequencies: [1x201 double]         
           B1GroupDelay: [1x201 double]         
                                                

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

samplesPerFrame = 2048;
wn = (2/samplesPerFrame) * (0:samplesPerFrame-1);
numRealPoints = samplesPerFrame/2 + 1;
htfe = dsp.TransferFunctionEstimator('FrequencyRange','onesided',...
      'SpectralAverages',64);
hplot = 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(Hellip, numRealPoints);
gdComp = grpdelay(Hgd, numRealPoints);
range = wn < wncomp(end);
gdExp = nan(numRealPoints, 1);
gdExp(range) = gdOrig(range) + gdComp(range);

% Stream random samples through filter cascade
Nframes = 300;
for k = 1:Nframes
    x = randn(samplesPerFrame,1);  % Input signal = white Gaussian noise

    y_orig = step(Hellip,x);       % Filter noise with original IIR filter
    y_corr = step(Hgd,y_orig);     % Compensating filter

    Txy = step(htfe,[x, x],[y_orig, y_corr]);
    gdMeas = HelperMeasureGroupDelay(Txy, [], 20);
    step(hplot, [gdMeas, gdExp]);
end

バンドパス イコライズ

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

% Design a bandpass Chebyshev type-1 filter and measure its total group
% delay over the passband.
Hcheby1 = 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(Hcheby1,wncomp,2);
g1 = max(g)-g;

% Design an 8th order arbitrary group delay filter. The pole radius is
% constrained to not exceed 0.95.
hgd = fdesign.arbgrpdelay('N,B,F,Gd',8,1,wncomp,g1);
Hgd = design(hgd,'iirlpnorm','MaxPoleRadius',0.95, 'SystemObject', true);

% Cascade the original filter with the compensation filter to achieve the
% desired group delay equalization.
% Verify by processing white noise and estimating the group delay at the
% two output stages
gdOrig = grpdelay(Hcheby1, numRealPoints);
gdComp = grpdelay(Hgd, numRealPoints);
range = wn > wncomp(1) & wn < wncomp(end);
gdExp = nan(numRealPoints,1); gdExp(range) = gdOrig(range) + gdComp(range);

release(hplot), hplot.YLimits = [0 55];
release(htfe)

% Stream random samples through filter cascade
for k = 1:Nframes
    x = randn(samplesPerFrame,1);  % Input signal = white Gaussian noise

    y_orig = step(Hcheby1,x);       % Filter noise with original IIR filter
    y_corr = step(Hgd,y_orig);     % Compensating filter

    Txy = step(htfe,[x, x],[y_orig, y_corr]);
    gdMeas = HelperMeasureGroupDelay(Txy, [], 20);
    step(hplot, [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 つの帯域 (バンド) を指定します。

% Design a bandstop Chebyshev type-2 filter and measure its total group
% delay over the passbands. Convert the measured group delay to seconds
% because fdesign.arbgrpdelay expects group delay specifications in seconds
% when you specify a sampling frequency.

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

% Design a 14th order arbitrary group delay allpass filter. The pole
% radius is constrained to not exceed 0.95. The group delay specifications
% are given in seconds and the frequency specifications are given in Hertz.

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

% Cascade the original filter with the compensation filter to process
% white noise and estimate the group delay at the two output stages

gdOrig = grpdelay(Hcheby2, numRealPoints);
gdComp = grpdelay(Hgd, 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(hplot),
    hplot.YLimits = [0 40];
    hplot.SampleIncrement = Fs/samplesPerFrame;
    hplot.YLabel = 'Group Delay (samples)';
    hplot.XLabel = 'Frequency (Hz)';
release(htfe)

% Stream random samples through filter cascade
Nframes = 300;
for k = 1:Nframes
    x = randn(samplesPerFrame,1);  % Input signal = white Gaussian noise

    y_orig = step(Hcheby2,x);       % Filter noise with original IIR filter
    y_corr = step(Hgd,y_orig);     % Compensating filter

    Txy = step(htfe,[x, x],[y_orig, y_corr]);
    gdMeas = HelperMeasureGroupDelay(Txy, [], 12);
    step(hplot, [gdMeas, gdExp]);
end

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

この情報は役に立ちましたか?