%《數字信號處理教程——MATLAB釋義與實現》第七章例7.4.2b程序hc742b
% 過渡帶設兩個樣本時的設計程序
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
N=input('N= (N必須為奇數）');wc=pi/3;       % 給出原始數據
N=N+mod(N+1,2);                             % 使N成為奇數
N1=fix(wc/(2*pi/N));N2=N-2*N1-1;      % 劃分通帶和阻帶的樣本點數
T1=input('過渡帶第一樣本T1的值');
T2=input('過渡帶第二樣本T2的值');
A=[ones(1,N1),T1,T2,zeros(1,N2-2),T2,T1,ones(1,N1-1)];	% 符幅特性樣本序列
theta = - pi*[0:N-1]*(N-1)/N;			% 相位特性樣本序列
H=A.*exp(j*theta);				% 頻率特性樣本序列
h=real(ifft(H));				% 反變換求出脈沖序列，去掉運算誤差造成的虛部
wp1=2*pi/N*fix(wc/(2*pi/N)-1);ws1=wp1+4*2*pi/N;
[db,mag,pha,grd,w] = myfreqz(h,[1]);    % 檢驗設計出的濾波器的頻率響應
[Ar,ww,type,L0]=amplres(h);             % 檢驗設計出的濾波器的符幅特性
dw = 2*pi/1000;                         % myfreqz函數中的頻率分辨率
Rp = -min(db(1:fix(wp1/dw)+1))      % 實際的通帶波動
As = -round(max(db(fix(ws1/dw)+1:501)))  % 最小阻帶衰減
l = 0:N-1; wl = (2*pi/N)*l;             % 由頻率樣本下標換成頻率樣本值 
wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];   % 繪制理想幅特性的頻率和幅度數據
subplot(2,2,1);plot(wl(1:N)/pi,A(1:N),'.',wdl,Adl); % 繪圖
axis([0,1,-0.1,1.1]); title('頻率樣本')
xlabel(''); ylabel('A(k)')
set(gca,'XTickMode','manual','XTick',chop([0,2/N,4/N,1],2)) % 畫出刻度線
set(gca,'YTickMode','manual','YTick',[0,1]); grid
subplot(2,2,2); stem(l,h,'.'); axis([-1,N,-0.1,0.5])
title('脈沖響應');ylabel('h(n)');text(N+1,-0.1,'n')
subplot(2,2,3); plot(ww/pi,Ar,wl(1:N)/pi,A(1:N),'.');
axis([0,1,-0.2,1.2]); title('符幅響應')
xlabel('頻率（單位：pi）'); ylabel('Ar(w)')
set(gca,'XTickMode','manual','XTick',chop([0,2/N,4/N,1],2)) % 畫出刻度線
set(gca,'YTickMode','manual','YTick',[0,1]); grid
subplot(2,2,4);plot(w/pi,db); axis([0,1,-80,10]); grid
title('幅度響應'); xlabel('頻率（單位：pi）');
ylabel('分貝數');
set(gca,'XTickMode','manual','XTick',chop([0,2/N,4/N,1],2)) % 畫出刻度線
set(gca,'YTickMode','Manual','YTick',[-As;0]);
set(gca,'YTickLabelMode','manual','YTickLabels',['As';' 0'])
set(gcf,'color','w');                                       % 置圖形背景色為白