БИХ-создание фильтра, учитывая предписанную групповую задержку

Открыть скрипт

В этом примере показано, как спроектировать произвольные фильтры групповой задержки с помощью разработчика фильтра fdesign.arbgrpdelay. Этот разработчик использует наименее ограниченный алгоритм оптимизации, чтобы спроектировать allpass БИХ-фильтры, которые соответствуют предписанной групповой задержке.

fdesign.arbgrpdelay может использоваться в эквализации групповой задержки.

Произвольная групповая задержка Filter Designer

Можно использовать fdesign.arbgrpdelay, чтобы спроектировать фильтр allpass с желаемым ответом групповой задержки. Желаемая групповая задержка задана в относительном смысле. Фактическая групповая задержка зависит от порядка фильтра (чем выше порядок, тем выше задержка). Однако, если вы вычитаете смещение в групповой задержке из-за порядка фильтра, групповая задержка спроектированного фильтра имеет тенденцию совпадать с желаемой групповой задержкой. Следующий код предоставляет пример с помощью двух различных порядка фильтра.

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

Обратите внимание на то, что в фильтре allpass, числитель всегда является обратной версией своего знаменателя. Поэтому вы не можете задать различный числитель и порядки знаменателя в fdesign.arbgrpdelay.

Следующий код показывает одной полосе произвольный проект групповой задержки с желаемыми значениями групповой задержки, Gd, в заданных точках частоты, F. В одном проектах полосы вы задаете групповую задержку по значениям частоты, которые покрывают целый интервал Найквиста [0 1] *pi рад/выборка.

arbGrpSpec = fdesign.arbgrpdelay('N,F,Gd',N,F,Gd) %#ok

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);

% 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(arbGrpSpec,arbGrpDelFilter1,Fpoints);

% Design another filter with an order equal to N2.
arbGrpSpec.FilterOrder = N2;
arbGrpDelFilter2 = design(arbGrpSpec,'MaxPoleRadius',R, ...
    'SystemObject', true);
M2 = measure(arbGrpSpec,arbGrpDelFilter2,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')

Следующий график показывает, что фактическая групповая задержка двух проектов отличается. Значение этого результата состоит в том, что нужно найти компромисс между лучшей подгонкой к желаемой относительной групповой задержке (меньше пульсации) и большей полной задержкой фильтра.

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

Эквализация групповой задержки полосы пропускания

Первичное использование fdesign.arbgrpdelay должно компенсировать нелинейные фазовые отклики БИХ-фильтров. Поскольку фильтр компенсации является allpass, он может быть расположен каскадом с фильтром, который вы хотите компенсировать, не влияя на его ответ величины. Поскольку каскад двух фильтров является самим БИХ-фильтром, он не может иметь линейной фазы (будучи устойчивым). Однако возможно иметь приблизительно линейный фазовый отклик в полосе пропускания полного фильтра.

Эквализация lowpass

Следующий пример использует fdesign.arbgrpdelay, чтобы компенсировать ответ групповой задержки lowpass эллиптический фильтр, не влияя на его ответ величины.

Вы используете многополосный проект, чтобы задать желаемые значения групповой задержки по одной или нескольким полосам интереса при оставлении групповой задержки всех других диапазонов частот незаданной (не заботьтесь об областях). В этом примере существует только одна полоса интереса, который равняется полосе пропускания фильтра lowpass. Вы хотите компенсировать групповую задержку этой полосы и не заботитесь о получившихся значениях групповой задержки за пределами него.

% Design an elliptic filter with a passband frequency of 0.2*pi
% rad/sample. Measure the total group delay over the passband.
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;

% Design an 8th order arbitrary group delay allpass filter. Use a
% multiband design and specify a single band.
allpassSpec = fdesign.arbgrpdelay('N,B,F,Gd',8,1,wncomp,g1) %#ok
allpassFilter = design(allpassSpec,'iirlpnorm', 'SystemObject', true);

allpassSpec = 

  arbgrpdelay with properties:

               Response: 'Arbitrary Group Delay'
          Specification: 'N,B,F,Gd'
            Description: {4x1 cell}
    NormalizedFrequency: 1
            FilterOrder: 8
                 NBands: 1
          B1Frequencies: [1x201 double]
           B1GroupDelay: [1x201 double]

Расположите каскадом исходный фильтр с фильтром компенсации, чтобы достигнуть желаемой эквализации групповой задержки. Проверьте путем обработки белого шума и оценки групповой задержки в этих двух выходных каскадах

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);

% Stream random samples through filter cascade
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 интервал рад/выборки. Как в предыдущем примере, существует только одна полоса интереса, который соответствует полосе пропускания фильтра. Поскольку вы хотите компенсировать групповую задержку этой полосы и не заботитесь о получившихся значениях групповой задержки за пределами него, вы используете многополосный проект и задаете одну полосу.

% Design a bandpass Chebyshev type-1 filter and measure its total group
% delay over the passband.
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;

% Design an 8th order arbitrary group delay filter. The pole radius is
% constrained to not exceed 0.95.
allpassSpec = fdesign.arbgrpdelay('N,B,F,Gd',8,1,wncomp,g1);
allpassFilter = design(allpassSpec,'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(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)

% Stream random samples through filter cascade
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

Получившийся фильтр имеет одну пару ограниченных полюсов. Изменение групповой задержки полосы пропускания ([0.3 0.4] *pi рад/выборка) меньше 0,2 выборок.

Эквализация Bandstop

Спроектируйте эквалайзер групповой задержки полосы пропускания для bandstop Чебышевского фильтра, действующего с частотой дискретизации 1 кГц. Заграждающий фильтр имеет две области полосы пропускания в [0 150] Гц и [200 500] интервалы Гц. Вы хотите компенсировать групповую задержку этих полос, таким образом, вы используете многополосный проект и задаете две полосы.

% 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;
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]);

% 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.

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);

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

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)

% Stream random samples through filter cascade
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

Получившийся фильтр имеет одну пару ограниченных полюсов. Полосы пропускания имеют изменение групповой задержки меньше чем 3 выборок.

Документация DSP System Toolbox
Поддержка
Для просмотра документации необходимо авторизоваться на сайте