Минимаксная оптимизация
В этом примере показано, как решить нелинейную задачу создания фильтра с помощью минимаксного алгоритма оптимизации, 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')
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');
Минимизируйте значения абсолютного максимума
Мы должны минимизировать значения абсолютного максимума, таким образом, мы устанавливаем options.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')
fig = gcf;
fig.NextPlot = 'replace';