指定群遅延が与えられた IIR フィルター設計
この例では、designfiltを使用して、任意の群遅延フィルターを設計する方法を示します。このデザイナーでは、最小 p 次制約付き最適化アルゴリズムを使用して、指定済みの群遅延を満たすオールパス IIR フィルターを設計します。
任意の群遅延応答を群遅延イコライゼーションに使用できます。
任意の群遅延フィルター デザイナー
arbgrpdelayiir 応答を designfilt と共に使用すると、希望する群遅延応答を指定してオールパス フィルターを設定できます。目的の群遅延は、相対的な意味で指定されます。実際の群遅延は、フィルター次数によって変わります (次数が高いほど遅延が大きい)。ただし、フィルター次数に起因して群遅延のオフセットを減算した場合、設計されたフィルターの群遅延は、目的の群遅延に一致する傾向があります。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
オールパス フィルターでは、分子が必ず対応する分母の反転となります。そのため、arbgrpdelayiir 応答を designfilt と共に使用した場合、分母と分子で異なる次数を指定できません。
以下のコードは、指定の周波数ポイント F で群遅延値 Gd をもつ任意のシングル バンド群遅延設計を示します。シングル バンドの設計では、ナイキスト間隔全体 [0 1]*pi rad/sample に対応する周波数値で群遅延を指定します。
arbGrpDelFilter1 = designfilt("arbgrpdelayiir",FilterOrder=N,Frequencies=F,GroupDelayResponse=Gd, ... MaxPoleRadius=R,SystemObject=true)
arbGrpDelFilter1 =
dsp.SOSFilter with properties:
Structure: 'Direct form II'
CoefficientSource: 'Property'
Numerator: [4×3 double]
Denominator: [4×3 double]
HasScaleValues: true
ScaleValues: [0.2506 1 1 1 1.4861]
Show all properties
指定された周波数ポイントにおける群遅延の合計を測定します。測定法を使用して、フィルターのノミナル群遅延を測定します。
Fpoints = 0:0.001:1; M1 = measure(arbGrpDelFilter1);
N2 に等しい次数をもつ別のフィルターを設計します。
arbGrpSpec.FilterOrder = N2; arbGrpDelFilter2 = designfilt("arbgrpdelayiir",FilterOrder=N2,Frequencies=F,GroupDelayResponse=Gd, ... MaxPoleRadius=R,SystemObject=true)
arbGrpDelFilter2 =
dsp.SOSFilter with properties:
Structure: 'Direct form II'
CoefficientSource: 'Property'
Numerator: [5×3 double]
Denominator: [5×3 double]
HasScaleValues: true
ScaleValues: [0.2026 1 1 1 1 2.0176]
Show all properties
M2 = measure(arbGrpDelFilter2);
測定した群遅延の合計からノミナル群遅延を引いたものをプロットします。
plot(M1.Frequencies,M1.TotalGroupDelay-M1.NomGrpDelay,... M2.Frequencies,M2.TotalGroupDelay-M2.NomGrpDelay,... [0 0.1 1],[2 3 1]); xlabel('Normalized Frequency (\times\pi rad/sample)') ylabel('Group delay (samples)') grid on legend('8th order design','10th order design','desired response')
以下のプロットに見られるように、2 つの構成の実際の群遅延は異なります。つまり、この結果では、1 つの設計は目的とする相対群遅延へのより近い近似 (より少ないリップル) とフィルター内での規模の大きい全体的な遅延の間でのバランスを探さなければならないことを意味しています。
FA = filterAnalyzer(arbGrpDelFilter1,arbGrpDelFilter2,Analysis="groupdelay"); setLegendStrings(FA, ["8th order design","10th order design"])
通過帯域群遅延イコライズ
arbgrpdelayiir 応答は主に、IIR フィルターの非線形位相応答を補正するために使用します。補正フィルターが オールパスであるため、振幅応答に影響を与えることなく、補正する対象のフィルターでカスケードを作成できます。2 つのフィルターのカスケードは IIR フィルターとなるため、線形位相はありません (安定している間)。ただし、フィルター全体の通過帯域内では、ほぼ線形位相応答である応答を得ることができます。
ローパス イコライズ
以下の例では、arbgrpdelayiir 応答を designfilt と共に使用して、振幅応答に影響を及ぼすことなく、ローパス楕円フィルターの群遅延応答をイコライズしています。
他のすべての周波数帯域の群遅延を指定せずに (「don't care」領域)、1 つ以上の対象帯域で目的の群遅延値を指定するには、マルチバンド設計を使用します。この例では、対象となるのは 1 つの帯域だけであり、ローパス フィルターの通過帯域に相当します。この帯域で群遅延を補正する必要がありますが、結果として群遅延値が帯域から外れていてもかまいません。
通過帯域周波数が ラジアン/サンプルの楕円フィルターを設計します。通過帯域全体で群遅延の合計を測定します。
ellipFilter = designfilt("lowpassiir", FilterOrder=4,PassbandFrequency=0.2, ... PassbandRipple=1,StopbandAttenuation=40, ... DesignMethod="ellip",SystemObject=true); wncomp = 0:0.001:0.2; g = grpdelay(ellipFilter,wncomp,2); % samples g1 = max(g)-g;
任意の 8 次群遅延オールパス フィルターを設計します。マルチバンド設計を使用して、シングル バンドを指定します。
allpassFilter = designfilt("arbgrpdelayiir",FilterOrder=8,NumBands=1, ... BandFrequencies1=wncomp,BandGroupDelayResponse1=g1, ... DesignMethod="lpnorm",SystemObject=true);
補正フィルターを使用して元のフィルターをカスケードし、目標の群遅延イコライズを達成します。ホワイト ノイズを処理し、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="Group Delay", ... ChannelNames=["Original", "Compensated","Expected Compensated"], ... 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 = designfilt("bandpassiir", FilterOrder=4, ... PassbandFrequency1=0.3,PassBandFrequency2=0.4, PassbandRipple=1, ... DesignMethod="cheby1",SystemObject=true); wncomp = 0.3:0.001:0.4; g = grpdelay(bandpassFilter,wncomp,2); g1 = max(g)-g;
任意の 8 次群遅延フィルターを設計します。極半径は 0.95 を超えないように制限されています。
allpassFilter = designfilt("arbgrpdelayiir",FilterOrder=8,NumBands=1, ... BandFrequencies1=wncomp,BandGroupDelayResponse1=g1, ... DesignMethod="lpnorm",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 フィルターを設計して、通過帯域全体でフィルターの群遅延の合計を測定します。arbgrpdelayiir はサンプリング周波数を秒単位で指定する群遅延仕様を前提としているため、測定された群遅延を秒に変換します。
Fs = 1e3; bandstopFilter = designfilt("bandstopiir",FilterOrder=6, ... StopbandFrequency1=150,StopbandFrequency2=400,StopbandAttenuation=1, ... SampleRate=Fs,DesignMethod="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 を超えないように制限されています。群遅延仕様は秒単位で指定され、周波数仕様はヘルツ単位で指定されています。
allpassFilter = designfilt("arbgrpdelayiir",FilterOrder=14,NumBands=2, ... BandFrequencies1=f1,BandGroupDelayResponse1=maxg-g1, ... BandFrequencies2=f2,BandGroupDelayResponse2=maxg-g2, ... SampleRate=Fs,DesignMethod="lpnorm",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 サンプル未満となります。
