Минимаксная оптимизация

Этот пример показов, как решить нелинейную задачу создания фильтра с помощью минимаксного алгоритма оптимизации, fminimax, в Optimization Toolbox™. Обратите внимание, что для запуска этого примера необходимо установить Signal Processing Toolbox™.

Установите параметры конечной точности

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

nbits  = 8;         % How many bits have we to realize filter 
maxbin = 2^nbits-1; % Maximum number expressable in nbits bits
n      = 4;         % Number of coefficients (order of filter plus 1)
Wn     = 0.2;       % Cutoff frequency for filter
Rp     = 1.5;       % Decibels of ripple in the passband
w      = 128;       % Number of frequency points to take

Сначала непрерывный проект

Это непрерывное создание фильтра; используем cheby1но мы также могли бы использовать ellip, yulewalk или remez здесь:

[b1,a1] = cheby1(n-1,Rp,Wn); 

[h,w] = freqz(b1,a1,w); % Frequency response
h = abs(h);             % Magnitude response
plot(w, h)
title('Frequency response using non-integer variables')

Figure contains an axes. The axes with title Frequency response using non-integer variables contains an object of type line.

x = [b1,a1];            % The design variables

Установите границы для коэффициентов фильтра

Теперь мы устанавливаем ограничения на максимальное и минимальное значения:

if (any(x < 0))
%   If there are negative coefficients - must save room to use a sign bit
%   and therefore reduce maxbin
    maxbin = floor(maxbin/2);
    vlb = -maxbin * ones(1, 2*n)-1;
    vub = maxbin * ones(1, 2*n); 
else
%   otherwise, all positive
    vlb = zeros(1,2*n); 
    vub = maxbin * ones(1, 2*n); 
end

Масштабные коэффициенты

Установите самое большое значение, равное maxbin, и соответственно масштабируйте другие коэффициенты фильтра.

[m, mix] = max(abs(x)); 
factor =  maxbin/m; 
x =  factor * x;    % Rescale other filter coefficients
xorig = x;

xmask = 1:2*n;
% Remove the biggest value and the element that controls D.C. Gain
% from the list of values that can be changed. 
xmask(mix) = [];
nx = 2*n;

Установите критерии оптимизации

Использование optimoptions, настроить критерии прекращения на достаточно высокие значения для содействия короткому времени работы. Также включите отображение результатов при каждой итерации:

options = optimoptions('fminimax', ...
    'StepTolerance', 0.1, ...
    'OptimalityTolerance', 1e-4,...
    'ConstraintTolerance', 1e-6, ...
    'Display', 'iter');

Минимизируйте абсолютные максимальные значения

Нам нужно минимизировать абсолютные максимальные значения, поэтому мы устанавливаем опции. MinAbsMax по количеству частотных точек:

if length(w) == 1
   options = optimoptions(options,'AbsoluteMaxObjectiveCount',w);
else
   options = optimoptions(options,'AbsoluteMaxObjectiveCount',length(w));
end

Исключите первое значение для оптимизации

Дискретизируйте и устраните первое значение и выполните оптимизацию по вызову FMINIMAX:

[x, xmask] = elimone(x, xmask, h, w, n, maxbin)
x = 1×8

    0.5441    1.6323    1.6323    0.5441   57.1653 -127.0000  108.0000  -33.8267

xmask = 1×6

     1     2     3     4     5     8

niters = length(xmask); 
disp(sprintf('Performing %g stages of optimization.\n\n', niters));
Performing 6 stages of optimization.
for m = 1:niters
    fun = @(xfree)filtobj(xfree,x,xmask,n,h,maxbin); % objective
    confun = @(xfree)filtcon(xfree,x,xmask,n,h,maxbin); % nonlinear constraint
    disp(sprintf('Stage: %g \n', m));
    x(xmask) = fminimax(fun,x(xmask),[],[],[],[],vlb(xmask),vub(xmask),...
        confun,options);
    [x, xmask] = elimone(x, xmask, h, w, n, maxbin);
end
Stage: 1 
                  Objective        Max     Line search     Directional 
 Iter F-count         value    constraint   steplength      derivative   Procedure 
    0      8              0    0.00329174                                            
    1     17      0.0001845      3.34e-07            1          0.0143     

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
Stage: 2 
                  Objective        Max     Line search     Directional 
 Iter F-count         value    constraint   steplength      derivative   Procedure 
    0      7              0     0.0414182                                            
    1     15        0.01649     0.0002558            1           0.261     
    2     23        0.01544     6.126e-07            1         -0.0282    Hessian modified  

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
Stage: 3 
                  Objective        Max     Line search     Directional 
 Iter F-count         value    constraint   steplength      derivative   Procedure 
    0      6              0     0.0716961                                            
    1     13        0.05943    -1.156e-11            1           0.776     

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
Stage: 4 
                  Objective        Max     Line search     Directional 
 Iter F-count         value    constraint   steplength      derivative   Procedure 
    0      5              0      0.129938                                            
    1     11        0.04278     2.937e-10            1           0.183     

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
Stage: 5 
                  Objective        Max     Line search     Directional 
 Iter F-count         value    constraint   steplength      derivative   Procedure 
    0      4              0     0.0901749                                            
    1      9        0.03867    -4.951e-11            1           0.256     

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.
Stage: 6 
                  Objective        Max     Line search     Directional 
 Iter F-count         value    constraint   steplength      derivative   Procedure 
    0      3              0       0.11283                                            
    1      7        0.05033    -1.249e-16            1           0.197     

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

Проверяйте ближайшие целочисленные значения

Проверьте, дают ли близлежащие значения лучший фильтр.

xold = x;
xmask = 1:2*n;
xmask([n+1, mix]) = [];
x = x + 0.5; 
for i = xmask
    [x, xmask] = elimone(x, xmask, h, w, n, maxbin);
end
xmask = 1:2*n;
xmask([n+1, mix]) = [];
x = x - 0.5;
for i = xmask
    [x, xmask] = elimone(x, xmask, h, w, n, maxbin);
end
if any(abs(x) > maxbin)
  x = xold; 
end

Сравнения частотной характеристики

Сначала строим график частотной характеристики фильтра и сравниваем его с фильтром, где коэффициенты просто округлены вверх или вниз:

subplot(211)
bo = x(1:n); 
ao = x(n+1:2*n); 
h2 = abs(freqz(bo,ao,128));
plot(w,h,w,h2,'o')
title('Optimized filter versus original')

xround = round(xorig)
xround = 1×8

     1     2     2     1    57  -127   108   -34

b = xround(1:n); 
a = xround(n+1:2*n); 
h3 = abs(freqz(b,a,128));
subplot(212)
plot(w,h,w,h3,'+')
title('Rounded filter versus original')

Figure contains 2 axes. Axes 1 with title Optimized filter versus original contains 2 objects of type line. Axes 2 with title Rounded filter versus original contains 2 objects of type line.

fig = gcf;
fig.NextPlot = 'replace';

См. также

Похожие темы