Произвольная групповая задержка 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)

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. Измерьте задержку именной группы фильтра с помощью метода меры.

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;

Спроектируйте 8-й порядок произвольная групповая задержка allpass фильтр. Используйте многополосный проект и задайте одну полосу.

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: [0 1.0000e-03 0.0020 0.0030 0.0040 0.0050 ... ]
           B1GroupDelay: [14.5838 14.5835 14.5824 14.5808 14.5784 ... ]

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 рад/отсчет) меньше 0,2 выборок.

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

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

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

Спроектируйте 14-й порядок произвольная групповая задержка allpass фильтр. Радиус полюса ограничивается не превысить 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 выборок.