%《數字信號處理教程——MATLAB釋義與實現》第七章例7.5.3a程序hc753a
% 用等波動法設計低通濾波器
% 電子工業(yè)出版社出版  陳懷琛編著 2004年9月
%
clear
wp = 0.2*pi; ws = 0.3*pi; Rp = 0.25; As = 50;   % 給定設計指標
delta1 = (10^(Rp/20)-1)/(10^(Rp/20)+1);         % 求通帶絕對波動指標
delta2 = (1+delta1)*(10^(-As/20));         %  求阻帶絕對波動指標
weights = [delta2/delta1 1];         % remez 函數要求的加權向量
deltaf = (ws-wp)/(2*pi);            % 估算N需要的過渡帶寬度
N = ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*deltaf)+1) % 估算N
N=N+mod(N-1,2);                 % N必須取奇數 
f = [0 wp/pi ws/pi 1];         %  函數要求的頻率向量
m = [1 1 0 0];                  % remez 函數要求的理想幅特性向量
Asd=10;                         % 為while語句預設一個較小的Asd
while Asd<As
    N = N+2                     % 若達不到要求，增加濾波器長度N，保持奇數
    h = remez(N-1,f,m,weights);        % 用最優(yōu)等波動法設計濾波器系數
    [db,mag,pha,grd,w] = myfreqz(h,[1]);        % 計算濾波器頻率特性
    dw = 2*pi/1000; wsi=ws/dw+1; wpi = wp/dw;  % 求邊界頻率下標
    Asd = -max(db(wsi:1:501));        % 求阻帶波動
end
[A,w1]=amplres(h);                  % 求符幅特性
subplot(2,2,1),stem([1:N],h,'.')    % 畫濾波器脈沖響應
subplot(2,2,2),plot(w,db)           % 畫分貝幅特性
subplot(2,2,3),plot(w,A)            % 畫符幅特性曲線
subplot(2,2,4),plot(w(1:wpi),A(1:wpi)-1,w(wsi:501),A(wsi:501))% 畫誤差波動曲線
set(gcf,'color','w')                % 置圖形背景色為白