Конструктор фильтра задержки произвольной группы

Можно использовать 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 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.BiquadFilter with properties:

                   Structure: 'Direct form II'
             SOSMatrixSource: 'Property'
                   SOSMatrix: [4x6 double]
                 ScaleValues: [5x1 double]
           InitialConditions: 0
    OptimizeUnityScaleValues: true

  Show all properties

Измерьте суммарную групповую задержку в наборе частотных точек от 0 до 1. Измерьте номинальную групповую задержку фильтра с помощью метода measure.

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

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

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

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

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

Выравнивание нижних частот

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

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

Спроектируйте эллиптический фильтр с частотой полосы пропускания $$\mathrm{}\mathrm{}\mathrm{0,2\delta }$$ рад/образец. Измерьте общую задержку группы по полосе пропускания.

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;

Создайте фильтр allpass произвольной групповой задержки 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: [1x201 double]
           B1GroupDelay: [1x201 double]

allpassFilter = design(allpassSpec,'iirlpnorm', 'SystemObject', true)

allpassFilter = 
  dsp.BiquadFilter with properties:

                   Structure: 'Direct form II'
             SOSMatrixSource: 'Property'
                   SOSMatrix: [4x6 double]
                 ScaleValues: [5x1 double]
           InitialConditions: 0
    OptimizeUnityScaleValues: true

  Show all properties

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

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] * пи рад/выборка. Как и в предыдущем примере, существует только одна интересующая полоса, которая соответствует полосе пропускания фильтра. Поскольку вы хотите компенсировать групповую задержку в этой полосе и не заботитесь о результирующих значениях групповой задержки за ее пределами, вы используете многополосную конструкцию и указываете одну полосу.

Спроектируйте полосовой фильтр Чебышева типа-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);

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

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

Результирующий фильтр имеет одну пару полюсов с ограничениями. Изменение групповой задержки в полосе пропускания ([0,3 0,4] * pi rad/выборка) составляет менее 0,2 выборки.

Выравнивание полосы пропускания

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

Спроектируйте полосовой фильтр Чебышева типа-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]);

Создайте фильтр allpass произвольной групповой задержки 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);

Объедините исходный фильтр с компенсационным фильтром для обработки белого шума и оценки групповой задержки на двух выходных каскадах.

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

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