双一次Z変換によるバタワースフィルタを利用したIIRフィルタの設計

ポイント

• プリワーピングによって、求めたいディジタルフィルタの遮断周波数fdをリファレンスとなるバタワースフィルタの遮断周波数faにあらかじめ変換しておく
• 伝達関数G(s)を求めたあと、s->zの双一次Z変換を行い、H(z)を求める

clear; close all;

T = 1e-4; %ディジタルフィルタのサンプリング周期
fs = 1/T; %ディジタルフィルタのサンプリング周波数
fd = 2500; %ディジタルフィルタの遮断周波数[Hz]
wd = 2*pi*fd; %ディジタルフィルタの遮断各周波数

wa = tan(wd*T/2)*2/T;%アナログフィルタの遮断角周波数
fa = wa/(2*pi); %hz: アナログフィルタの遮断周波数
disp('***プリワーピング処理***');
disp(sprintf('IIRディジタルフィルタの設計遮断周波数fd: %g[Hz]',fd));
disp(sprintf('等価アナログフィルタの要求遮断周波数fa: %g[Hz]',fa));

Ns = [15 20 40];

syms s;
Has = [];
figure(1);
lab = sprintf('バタワースフィルタの振幅特性: 遮断周波数%g[Hz]',fa);
title(lab);
xlabel('f[Hz]');
hold on;
figure(2);
title('バタワースフィルタの位相特性');
xlabel('f[Hz]');
hold on;
Hs = s*zeros(size(Ns));
for idx=1:length(Ns)
    N = Ns(idx);
    omega = linspace(0,wa*1.9,100);
    Ha = sqrt(1./(1+(omega/wa).^(2*N)));
    figure(1);
    subplot(211);
    hold on;
    lab1 = sprintf('(N=%d)',N);    
    plot(omega/(2*pi), Ha, 'DisplayName',lab1);
    D = 1;

    phi = 0;
    for k=1:N
        sigma_k = -sin((2*k-1)/(2*N)*pi);
        omega_k =  cos((2*k-1)/(2*N)*pi);
        phi = phi - atan((omega_k-(omega/wa))/sigma_k);
        sk = sigma_k + j*omega_k;
        D = D * (s-sk);
    end
    H = 1/D;
    Hs(idx) = H;
    figure(2);
    subplot(211)
    hold on;
    lab2 = sprintf('(N=%d)',N);
    plot(omega/(2*pi), rad2deg(phi), 'DisplayName',lab2);
end
figure(1);
subplot(211);
lab = sprintf('バタワースフィルタの振幅特性: 遮断周波数%g[Hz]',fa);
title(lab);
xlabel('f[Hz]');
legend();

figure(2);
subplot(211);
title('バタワースフィルタの位相特性');
xlabel('f[Hz]');
legend();


for idx=1:length(Ns)
N = Ns(idx);
H1 = Hs(idx);
syms z;
s_by_z = 2/T*(1-z^(-1))/(1+z^(-1));
Hz1 = subs(H1, s, s_by_z);
Hz1 = vpa(simplify(Hz1));
%ディジタルフィルタのwdでωを規格化
z_exp_jwt = exp(j*(omega/wd)*T);
Hz1_w = subs(Hz1, z, z_exp_jwt);
Hz1_w_abs = abs(Hz1_w);
figure(1);
subplot(212);
hold on;
plot(omega/(2*pi), Hz1_w_abs,'DisplayName',sprintf('N=%d',N));
title(sprintf('等価IIRディジタルフィルタの振幅特性: 遮断周波数%g[Hz]',fd));
xlabel('f[Hz');
legend();


figure(2);
subplot(212);
hold on;
%以下の計算はN=50において終了しない
angle_Hz1_z_exp_jwt = subs(angle(Hz1),z,z_exp_jwt);
unwrap_angle = smoother(angle_Hz1_z_exp_jwt);
plot(omega/(2*pi), rad2deg(unwrap_angle),'DisplayName',sprintf('N=%d',N));
title('等価IIRディジタルフィルタの位相特性');
xlabel('f[Hz]');
legend();
end


legend();
big;



function ys=smoother(xs)
ys = xs;
for idx = 2:length(xs)
    if (xs(idx) - xs(idx-1)) > 0
        disp('curernt value increased!');
        ys(idx:end) = ys(idx:end) - abs(xs(idx)-xs(idx-1));
    end
end
end