%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IIRフィルタ (ハイパスフィルタの設計)
%
%                                                      2003.12.10 Katto Lab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

% 変数設定
N = 5;              % Butterworth フィルタの次数
Fs = 1;             % サンプリング周波数 (Hz単位)
W0 = 2;             % アナログハイパスフィルタのカットオフ周波数
cut = W0/(2*pi);    % ディジタルフィルタのカットオフ周波数

% プロトタイプアナログフィルタの作成
%% N 次の Butterworth ローパスフィルタ
[z,p,k] = buttap(N);
[num1,den1] = zp2tf(z,p,k);
%% カットオフ周波数 W0 のハイパスフィルタに変換
[num2,den2] = lp2hp(num1,den1,W0);
%% 周波数特性と正規化
[h0,w0] = freqs(num1,den1);
[h1,w1] = freqs(num2,den2);
w0 = w0 / Fs;
w1 = w1 / Fs;

% 双一次変換(1) (プリワーピングなし)
%% 双一次変換
[bz,az] = bilinear(num2,den2,Fs);
[h2,w2] = freqz(bz,az);

% 双一次変換(2) (プリワーピングあり)
%% プリワーピング
W1 = 2*Fs*tan(cut*(2*pi/Fs)/2);
[num3,den3] = lp2hp(num1,den1,W1);
%% 双一次変換
[bz,az] = bilinear(num3,den3,Fs);
[h3,w3] = freqz(bz,az);

% グラフ表示
figure(1);
plot(w1,abs(h1),'r'); hold on;
plot(w2,abs(h2),'g'); hold on;
plot(w3,abs(h3),'b'); hold on;
axis([0.0, pi, 0.0, 1.2]);
legend('Butterworth HPF', 'IIR (without prewarping)', 'IIR (with prewarping)');
xlabel('normalized frequency');
ylabel('frequency response');

%figure(2);
%subplot(3,1,1);
%plot(w0,abs(h0),'c-'); hold on;
%plot(w1,abs(h1),'r'); hold on;
%axis([0, pi, 0, 1.2]);
%title('Butterworth LPF and HPF'); 
%legend('Butterworth LPF', 'Butterworth HPF');
%xlabel('normalized frequency');
%subplot(3,1,2);
%plot(w2,abs(h2));
%axis([0, pi, 0, 1.2]);
%title('IIR (without prewarping)')
%xlabel('normalized frequency');
%subplot(3,1,3);
%plot(w3,abs(h3));
%axis([0, pi, 0, 1.2]);
%title('IIR (with prewarping)')
%xlabel('normalized frequency');



