close all;                                              % 모든 생성된 창을 닫음 
clear all;                                              % 열려있는 작업영역(workspace)을 모두 비움 
clc;                                                    % 명령(command) 창을 비움

r=0.88; w0=pi/2;                                        % 필터 파라미터 설정
dw=pi/40;                                               % 주파수 간격 설정
wr=dw*[17:29];                                          % 각 공진기의 공진 주파수 설정 
for k=1:7
   b0=sqrt((1-r)^2*(1-2*r*cos(2*wr(k))+r^2));           % b0 계산 
   den(k,1:3)=[1 -2*r*cos(wr(k)) r^2];                  % 공진기 전달 함수 분모 계수
   num(k,1:3)=[1 0 -1]*b0;                                  % 공진기 전달 함수 분자 계수
end
num(1,1:3)=num(1,1:3)/2; num(7,1:3)=num(1,1:3);         % 공진기 전달 함수 분자 계수
M=512; Hf=zeros(1,M);                                   % 각 공진기 주파수 응답 초기화
for k=1:7
   [Hfr(k,1:M),w]=freqz(num(k,1:3), den(k,1:3),M);      % 각 공진기 주파수 응답 계산 
   Hf=Hf+Hfr(k,1:M);                                    % 대역 통과 필터 주파수 응답 계산
end
h=ifft(Hf);                                             % IDTFT로 임펄스 응답 계산 
n=0:30;                                                 % 시간출 설정

subplot(3,1,1);                                         % 3행 1열 분할 그림 창의 1번 창
plot(w/(2*pi),abs(Hf));                                 % 대역 통과 필터의 진폭 응답을 그림
title('\bf{대역 통과  필터의 진폭 응답}');                % 그림 제목
xlabel('\bf{F}'); grid on;                              % x축 라벨 및 그리드 표시 
subplot(3,1,2);                                         % 3행 1열 분할 그림 창의 2번 창
plot(w/(2*pi),angle(Hf));                               % 대역 통과 필터의 위상 응답을 그림
%title('\bf{대역 통과 필터의 위상 응답}');                % 그림 제목
xlabel('\bf{F}'); grid on;                              % x축 라벨 및 그리드 표시 
subplot(3,1,3);                                         % 3행 1열 분할 그림창의 3번 창
stem(n,h(1:31));                                        % 대역 통과 필터의 임펄스 응답을 그림
title('\bf{대역 통과 필터의 임펄스 응답}');               % 그림 제목
xlabel('\bf{n}'); grid on;                              % x축 라벨 및 그리드 표시 