Произвольная групповая задержка 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 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, он может быть каскадирован с фильтром, который вы хотите компенсировать, не влияя на его амплитудную характеристику. Поскольку каскад двух фильтров является самим БИХ, он не может иметь линейную фазу (будучи стабильным). Однако возможно иметь приблизительно линейный фазовый отклик в полосе пропускания всего фильтра.

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

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

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

Проектируйте эллиптический фильтр с частотой полосы пропускания $$0.2\pi$$ рад/образец. Измерьте общую задержку группы по полосе пропускания.

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

Спроектировать полосу пропускания Чебышевского фильтра 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/sample) составляет менее 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. Спецификации групповой задержки приведены в секундах, а спецификации частоты приведены в 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);

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

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 выборки.