%《數字信號處理教程——MATLAB釋義與實現》第七章例7.5.4程序hc754
% 用等波動法設計帶通濾波器
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
ws1 = 0.2*pi; wp1 = 0.35*pi; wp2 = 0.65*pi; ws2 = 0.8*pi;  % 給定設計指標
Rp = 1.0; As = 60;
delta1 = (10^(Rp/20)-1)/(10^(Rp/20)+1);       % 求通帶絕對波動指標
delta2 = (1+delta1)*(10^(-As/20));       %  求阻帶絕對波動指標
deltaH = max(delta1,delta2); deltaL = min(delta1,delta2);
weights = [1 delta2/delta1 1];        % remez 函數要求的加權向量
delta_f =min((ws2-wp2)/(2*pi), (wp1-ws1)/(2*pi));      % 估算N需要的過渡帶寬度
N = ceil((-20*log10(sqrt(delta1*delta2))-13)/(14.6*delta_f)+1) 		% 估算N
f = [0 ws1/pi wp1/pi wp2/pi ws2/pi 1];       %  函數要求的頻率向量
m = [0 0 1 1 0 0];         % remez 函數要求的理想幅特性向量
h = remez(N-1,f,m,weights);      % 用最優等波動法設計濾波器系數
[db,mag,pha,grd,w] = myfreqz(h,[1]);        % 計算濾波器頻率特性
dw=2*pi/1000;					 % 求邊界頻率下標
ws1i=floor(ws1/dw)+1; wp1i = floor(wp1/dw)+1;
ws2i=floor(ws2/dw)+1; wp2i = floor(wp2/dw)+1;
Asd = -max(db(1:1:ws1i))      % 求阻帶波動
N=N+1;
h = remez(N-1,f,m,weights);      % 用最優等波動法設計濾波器系數
[db,mag,pha,grd,w] = myfreqz(h,[1]);        % 計算濾波器頻率特性
Asd = -max(db(1:1:ws1i))      % 求阻帶波動
[Hr,omega,P,L] = amplres(h);
% 畫圖
figure(1); subplot(1,1,1)
subplot(2,2,1); stem([0:1:N-1],h,'.'); title('實際脈沖響應')
axis([0,N-1,-0.4,0.5]);ylabel('h(n)');text(N,-0.4,'n')
set(gca,'XTickMode','manual','XTick',[0,N-1])               % 繪制刻度線
set(gca,'YTickMode','manual','YTick',[-0.4:0.2:0.5])
subplot(2,2,2);plot(w/pi,db);title('幅度響應');
axis([0,1,-80,10]); xlabel(''); ylabel('分貝數')
set(gca,'XTickMode','manual','XTick',f)                     % 繪制刻度線
set(gca,'YTickMode','manual','YTick',[-60,0]);
set(gca,'YTickLabelMode','manual','YTickLabels',['60';' 0']);grid
subplot(2,2,3);plot(omega/pi,Hr);title('振幅響應');
axis([0 1 -0.1 1.1]); xlabel('頻率（單位：pi）'); ylabel('Hr(w)')
set(gca,'XTickMode','manual','XTick',f)                    % 繪制刻度線
set(gca,'YTickMode','manual','YTick',[0,1]);grid
dw = 2*pi/1000; sp_edge1 = ws1/dw+1; sp_edge2 = ws2/dw+1;
subplot(2,2,4);
sb1w = omega(1:1:ws1i)/pi; sb1e = Hr(1:1:ws1i); 
pbw = omega(wp1i:wp2i)/pi; pbe = Hr(wp1i:wp2i)-1;
sb2w = omega(ws2i:501)/pi; sb2e = Hr(ws2i:501); 
plot(sb1w,sb1e,pbw,pbe*(delta2/delta1),sb2w,sb2e);          % 繪制幅度誤差曲線
title('加權誤差');
axis([0,1,-deltaL,deltaL]);
xlabel('頻率（單位：pi）'); ylabel('Hr(w)')
set(gca,'XTickMode','manual','XTick',f)
set(gca,'YTickMode','manual','YTick',[-deltaL,0,deltaL]);
set(gca,'XGrid','on','YGrid','on')
set(gcf,'color','w');                                       % 置圖形背景色為白