%Step 1：計算幅值
function y = filtmin(a,w)
n = length(a);
y = cos(w'*(0:n-1)*2*pi)*a ;

%Step 2：調用優化函數
% 繪制初始系數的濾波器頻響圖
a0 = ones(15,1);
incr = 50;
w = linspace(0,0.5,incr);

y0 = filtmin(a0,w);
clf, plot(w,y0,'-:');
drawnow;

% 構造多目標逼近問題
w1 = linspace(0,0.1,incr) ;
w2 = linspace(0.15,0.5,incr);
w0 = [w1 w2];
goal = [1.0*ones(1,length(w1)) zeros(1,length(w2))];
weight = ones(size(goal)); 

% 調用 fgoalattain函數
options = optimset('GoalsExactAchieve',length(goal));
[a,fval,attainfactor,exitflag] = fgoalattain(@filtmin,...
    a0,goal,weight,[],[],[],[],[],[],[],options,w0);
% 繪制優化系數的濾波器頻響圖
y = filtmin(a,w);
hold on, plot(w,y,'r')
axis([0 0.5 -3 3])
xlabel('Frequency (Hz)')
ylabel('Magnitude Response (dB)')
legend('initial', 'final')
grid on