В этом примере показано, как решить проблему проектирования нелинейного фильтра с помощью алгоритма оптимизации minimax. fminimax, в Toolbox™ оптимизации. Следует отметить, что для выполнения этого примера необходимо установить 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');
Нам нужно минимизировать абсолютные максимальные значения, поэтому мы устанавливаем опции. 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';