Проектирование фильтра нижних частот с узкой полосой пропускания

Одним из недостатков использования КИХ-фильтров является то, что порядок фильтров имеет тенденцию расти обратно пропорционально полосе пропускания перехода фильтра. Рассмотрим следующие проектные характеристики (где пульсации приведены в линейных единицах):

Fpass = 2866.5;  % Passband edge
Fstop = 3087;    % Stopband edge
Apass = 0.0174;  % Passband ripple, 0.0174 dB peak to peak
Astop = 66.0206; % Stopband ripple, 66.0206 dB of minimum attenuation
Fs    = 44.1e3;

lowpassSpec = fdesign.lowpass(Fpass,Fstop,Apass,Astop,Fs);

Регулярная линейно-фазовая равнодействующая конструкция, соответствующая спецификациям, может быть спроектирована с помощью:

eqripFilter = design(lowpassSpec,'equiripple','SystemObject',true);
cost(eqripFilter)

ans = struct with fields:
                  NumCoefficients: 695
                        NumStates: 694
    MultiplicationsPerInputSample: 695
          AdditionsPerInputSample: 694

Требуемая длина фильтра оказывается равной 694 отводам.

Конструкция интерполированной FIR (IFIR)

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

interpFilter = design(lowpassSpec,'ifir','SystemObject',true);

Видимо, мы усугубили ситуацию. Вместо одного фильтра с 694 умножителями мы теперь имеем два фильтра с общим количеством 804 умножителей. Однако тщательное изучение второго этапа показывает, что только один множитель в 5 является ненулевым. Фактическое общее число множителей сократилось с 694 до 208.

cost(interpFilter)

ans = struct with fields:
                  NumCoefficients: 208
                        NumStates: 802
    MultiplicationsPerInputSample: 208
          AdditionsPerInputSample: 206

Давайте сравним ответы двух конструкций:

fvt = fvtool(eqripFilter,interpFilter,'Color','White');
legend(fvt,'Equiripple design', 'IFIR design','Location','Best')

Управление коэффициентом повышения дискретизации вручную

В предыдущем примере мы автоматически определили используемый коэффициент повышения дискретизации таким образом, что общее количество множителей было минимизировано. Оказалось, что для заданных спецификаций оптимальный коэффициент повышения дискретизации составлял 5. Однако если мы рассмотрим варианты конструкции:

opts = designopts(lowpassSpec,'ifir')

opts = struct with fields:
      FilterStructure: 'dffir'
     UpsamplingFactor: 'auto'
    JointOptimization: 0
         SystemObject: 0

мы видим, что можем контролировать повышающий коэффициент дискретизации. Например, если мы хотели увеличить число на 4, а не на 5:

opts.UpsamplingFactor = 4;
opts.SystemObject = true; 
upfilter = design(lowpassSpec,'ifir',opts);
cost(upfilter)

ans = struct with fields:
                  NumCoefficients: 217
                        NumStates: 767
    MultiplicationsPerInputSample: 217
          AdditionsPerInputSample: 215

Мы бы получили конструкцию, которая имеет в общей сложности 217 ненулевых множителей.

Использование совместной оптимизации

Можно проектировать два фильтра, используемых в ИФИР совместно. Таким образом, мы можем сэкономить значительное количество множителей за счет более длительного времени проектирования (из-за природы алгоритма, конструкция также может не сходиться полностью в некоторых случаях):

opts.UpsamplingFactor = 'auto'; % Automatically determine the best factor
opts.JointOptimization = true;
ifirFilter = design(lowpassSpec,'ifir',opts);
cost(ifirFilter)

ans = struct with fields:
                  NumCoefficients: 172
                        NumStates: 726
    MultiplicationsPerInputSample: 172
          AdditionsPerInputSample: 170

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

Использование многоступенчатых/многоступенчатых методов для создания эффективных проектов

Для конструкций, обсуждавшихся до сих пор, использовались односкоростные методы. Это означает, что количество необходимых умножений на входную выборку (MPIS) равно количеству ненулевых умножителей. Например, последний дизайн, который мы показали, требует 152 MPIS. Одноступенчатая равнодействующая конструкция, которую мы начали с требуемого 694 MPIS.

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

multistageFilter = design(lowpassSpec,'multistage','SystemObject',true)

multistageFilter = 
  dsp.FilterCascade with properties:

    Stage1: [1x1 dsp.FIRDecimator]
    Stage2: [1x1 dsp.FIRDecimator]
    Stage3: [1x1 dsp.FIRDecimator]
    Stage4: [1x1 dsp.FIRInterpolator]
    Stage5: [1x1 dsp.FIRInterpolator]
    Stage6: [1x1 dsp.FIRInterpolator]

cost(multistageFilter)

ans = struct with fields:
                  NumCoefficients: 396
                        NumStates: 352
    MultiplicationsPerInputSample: 73
          AdditionsPerInputSample: 70.8333

Первая ступень имеет 21 умножитель, а коэффициент прореживания - 3. Поэтому количество MPIS равно 7. Вторая ступень имеет длину 45 и коэффициент кумулятивного прореживания 6 (то есть коэффициент прореживания этой ступени, умноженный на коэффициент прореживания первой ступени; это происходит потому, что входные выборки на этом этапе уже поступают со скоростью 1/3 скорости входных выборок на первый этап). Среднее число умножений на входную выборку (ссылка на вход всего многоскоростного/многоступенчатого фильтра), таким образом, равно 45/6 = 7,5. Наконец, учитывая, что третья ступень имеет коэффициент прореживания 1, среднее число умножений на вход для этой ступени равно 130/6 = 21,667. Общее число средних MPIS для трех прореживателей составляет 36,5.

Для интерполяторов оказывается, что фильтры идентичны прореживающим. Более того, их вычислительные затраты одинаковы. Поэтому общее количество MPIS для всей многоступенчатой/многоступенчатой конструкции составляет 73.

Теперь мы сравним ответы эквириптной конструкции и этой:

fvt = fvtool(eqripFilter,multistageFilter,'Color','White', ...
    'NormalizeMagnitudeto1','on');
legend(fvt,'Equiripple design', 'Multirate/multistage design', ...
    'Location','NorthEast')

Обратите внимание, что затухание полосы останова для многоступенчатой конструкции примерно вдвое больше, чем для других конструкций. Это связано с тем, что прореживатели должны ослаблять внеполосные компоненты на требуемые 66 дБ, чтобы избежать наложения, которое нарушило бы требуемые спецификации. Аналогичным образом, интерполяторы должны ослаблять изображения на 66 дБ, чтобы соответствовать спецификациям этой проблемы.

Управление количеством этапов вручную

Многоступенчатая/многоступенчатая конструкция, которая была получена, состояла из 6 стадий. Количество этапов определяется автоматически по умолчанию. Однако также можно вручную управлять количеством этапов, которые приводят к этому. Например:

lp4stage=design(lowpassSpec,'multistage','NStages',4,'SystemObject',true);
cost(lp4stage)

ans = struct with fields:
                  NumCoefficients: 516
                        NumStates: 402
    MultiplicationsPerInputSample: 86
          AdditionsPerInputSample: 84.5000

Среднее количество MPIS для этого случая составляет 86

Задержка группы

Мы можем вычислить групповую задержку для каждой конструкции. Обратите внимание, что многоступенчатая/многоступенчатая конструкция вносит наибольшую задержку (это цена за менее дорогостоящую в вычислительном отношении конструкцию). Конструкция IFIR обеспечивает большую задержку, чем одноступенчатая, но меньшую, чем многоступенчатая/многоступенчатая конструкция.

fvt = fvtool(eqripFilter,interpFilter,multistageFilter,...
    'Analysis','grpdelay','Color','White');
legend(fvt, 'Equiripple design','IFIR design',...
    'Multirate/multistage design');

Фильтрация сигнала

Теперь мы показываем на примере, что IFIR и многоступенчатая/многоскоростная конструкция работают сравнимо с одноступенчатой, но требуют гораздо меньшего вычисления. Для выполнения фильтрации каскадов необходимо вызвать generateFilterCode (ifirFilter) и generateFilterCode (multistageFilter). Это уже сделано для создания HelperIFIRFilter и HelperMultiFIRFilter.

scope = dsp.SpectrumAnalyzer('SpectralAverages',50,'SampleRate',Fs,...
    'PlotAsTwoSidedSpectrum',false,'YLimits', [-90 10],...
    'ShowLegend',true','ChannelNames',{'Equiripple design',...
    'IFIR design','Multirate/multistage design'});
tic,
while toc < 20
    % Run for 20 seconds
    x = randn(6000,1);
    
    % Filter using single stage FIR filter
    y1 = eqripFilter(x);
    
    % Filter using IFIR filter
    y2 = HelperIFIRFilter(x);
    
    % Filter multi-stage/multi-rate FIR filters
    y3 = HelperMultiFIRFilter(x);
    
    % Compare the output from both approaches
    scope([y1,y2,y3])
end