close all,clear all,clc;

N=(40-7.95)/2.286/(0.15*pi);
N=round(N)+1                          %N取奇數
omega=[0.15*pi,0.3*pi,0.5*pi,0.65*pi];
wc1=(omega(1)+omega(2))/2;            %理想帶通的截止頻率
wc2=(omega(3)+omega(4))/2;
a=(N-1)/2;                            %滿足線性相位的條件

for n=0:100                           %此處我只畫出理想帶通濾波器101個點的單位樣值
    if n==a                           %響應，當然還可以取其它值
        hd(n+1)=(wc2-wc1)/pi;
    else
    hd(n+1)=(sin((n-a)*wc2)-sin((n-a)*wc1))/pi/(n-a);  %由頻率特性求單位樣值相應的公式
    end
end
n=0:100;
stem(n,hd);

ns=40;
b=0.5842*(ns-21)^0.4+0.07886*(ns-21)
for m=0:N-1
    b1=b*sqrt(1-(1-2*m/(N-1))^2);
    wk(m+1)=besl(b1)/besl(b);       %調用besl函數，求得kaiser窗函數
end
m=0:N-1;
figure;
stem(m,wk);                         %畫出窗函數

for n=0:N-1
    h(n+1)=hd(n+1)*wk(n+1);         %加窗后的帶通濾波器的單位樣值響應
end
figure;
stem(m,h);


k=2*pi/(pi/40)+1;
H=zeros(1,k);                       %對數字帶通濾波器頻率特性初始化
q=0;
for w=-pi:pi/40:pi
    q=q+1;
    for n2=0:N-1
        H(q)=H(q)+h(n2+1)*exp(-j*w*n2);%計算數字濾波器的幅頻特性
    end
end
figure,hold on,grid on;
H=20*log10(abs(H));                    %轉換為對數幅頻特性
w=-pi:pi/40:pi;
w=w/pi;
plot(w,H);