% zi_lcfilter.m UTF-8でテキストファイルとし、拡張子.mで保存してください。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LC filter design with zero soure impedance or infinite load impedance
% for Butterworth, Chebyshev, and Bessel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
% 電力伝達多項式をプロットするか
PLOT1 = 0; % 1 にした場合にプロットし、0だとプロットしない
fprintf(char(10))
fprintf('LC filter design with zero soure impedance or infinite load impedance')
fprintf(char(10))
fprintf('for Butterworth, Chebyshev, and Bessel')
fprintf(char(10))
fprintf(char(10))
fprintf('Butterworth = 1')
fprintf(char(10))
fprintf('Chebyshev = 2')
fprintf(char(10))
fprintf('Bessel = 3')
fprintf(char(10))
fprintf(char(10))
%--------------------------------------------------------
type = 0;
st = "Enter filter type number: ";
while !(type == 1 || type == 2 || type == 3)
type = input(st);
end
% フィルタの次数を設定
order = 0;
switch(type)
case 2
st = "Enter filter order (Odd number) 3 to 9: ";
otherwise
st = "Enter filter order 2 to 9: ";
end
while !(order > 1 && order < 10)
order = input(st);
if ((type == 2) && mod(order, 2) == 0)
fprintf(char(10))
fprintf('Error!')
fprintf(char(10))
fprintf('Even order Chebyshev filter can not be generated.')
fprintf(char(10))
fprintf('Because even order should have assigned ripple level attenuation at DC.')
fprintf(char(10))
fprintf('Enter odd order number please.')
fprintf(char(10))
fprintf(char(10))
order = 0;
end
end
% ChebyshevフィルタのリプルをdBで設定
ripple = 0;
switch(type)
case 2
st = "Enter Chebyshev ripple in dB (0.01 to 3): ";
while !(ripple >= 0.01 && ripple <= 3)
ripple = input(st);
end
end
fprintf(char(10)')
fprintf(char(10)')
switch(type)
case 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Butterworth filter カーブ
% Butterworth特性を自乗し、それを逆数とした電力伝達多項式を得る
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i = sqrt(-1);
powerpoly = [1 zeros(1, (2 * order - 1)) 1];
ripple = 1; % PLOT 1でプロットするための一時的な設定
case 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Chebyshev filter カーブ
% Chebyshev特性を自乗し、それを逆数とした電力伝達多項式を得る
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i = sqrt(-1);
% dBで入力したChebyshev リプルを実電力リプル率にする
e = 10^(ripple/10);
e = e - 1;
e = sqrt(e);
% Chebyshev多項式を漸化式として得る
% 漸化式の関係は
% T_n+1(x) = 2x T_n(x) - T_n-1(x)
tminus1 = [zeros(1, order) 1];
tminus0 = [tminus1 0];
tminus0 = tminus0(2:order+2);
for n = 2:order
tn = 2 * [tminus0 0];
tn = tn(2:order+2);
tn = tn - tminus1;
tminus1 = tminus0;
tminus0 = tn;
end
% Tn^2 （電力）を得る
tn2 = conv(tn, tn);
powerpoly = e^2 * tn2;
powerpoly = powerpoly + [zeros(1, order*2) 1];
case 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bessel filter カーブ
% Bessel 特性を自乗し、それを逆数とした電力伝達多項式を得る
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i = sqrt(-1);
% ベッセルフィルタはいきなりs平面での多項式で以下のようにあらわされる
powerpolyP = [];
for n = 0:order
powerpolyP = [powerpolyP (factorial(order + n)/(factorial(order - n)*factorial(n)))*(1/2)^n*i^(order - n)];
end
% 0次の次数で正規化
powerpolyP = powerpolyP/powerpolyP(1, order + 1);
% s平面での多項式の共役を得る
powerpolyN = [];
for n = 0:order
powerpolyN = [powerpolyN (factorial(order + n)/(factorial(order - n)*factorial(n)))*(1/2)^n*(-i)^(order - n)];
end
% 0次の次数で正規化
powerpolyN = powerpolyN/powerpolyN(1, order + 1);
% 自乗して電力伝達多項式を得る
powerpoly = conv(powerpolyP, powerpolyN);
powerpoly = real(powerpoly);
ripple = 1; % PLOT 1でプロットするための一時的な設定
% ω = 1 rad/secにおける減衰量を計算する
mpx = [];
for n = 2*order:-1:0
mpx = [(-1)^n mpx];
end
att = powerpoly * mpx';
att = -10 * log10(att);
fprintf('Attenuation at omega = 1 is %d2.3 dB', att)
fprintf(char(10))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ここまでで伝達多項式powerpolyを得た。
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT 1
% 電力伝達多項式を用いてプロットする
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
if PLOT1 % 1 にした場合にプロットし、0だとプロットしない
s = i * logspace(-2, 3, 10001); % s = jωとする
accum_num = freqresponse(H_s_num_square, s);
accum_denom = freqresponse(H_s_square, s);
plotacccum = 10*log10(abs(accum_num)) - 10*log10(abs(accum_denom));
semilogx(abs(s), plotacccum)
xlabel('Normalized Angular Frequency')
ylabel('Transfer Level [dB]')
grid on
topofY = plotacccum(1,1);
topofY = ceil(topofY/ripple)*ripple + 3 * ripple;
botofY = topofY - 10 * ripple;
yticks([botofY:ripple:topofY]);
set(gca, 'FontSize', 18);
axis([0.01 100 botofY topofY])
figure
semilogx(abs(s), plotacccum)
xlabel('Normalized Angular Frequency')
ylabel('Transfer Level [dB]')
grid on
set(gca, 'FontSize', 18);
axis([0.01 100 -80 5])
end
%--------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% フィルタ設計
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i = sqrt(-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Butterworth, Chebyshev and Besselフィルタ設計
% 以下で必要となる変数: powerpoly, order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% powerpolyのs平面での根の左半分を得る
rootsofpower = roots(powerpoly);
rootsofpower = i * rootsofpower; % 90°回転させ、s平面に変換する
rootsofpowerforsort = [real(rootsofpower) rootsofpower];
rootsofpowerforsort = sortrows(rootsofpowerforsort, 1);
rootsofpower = rootsofpowerforsort(1:order, 2);
H_s = poly(rootsofpower);
H_s = real(H_s); % これがs平面での電圧伝達多項式（√電力）の分母になる
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ラダー素子定数を計算する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num = [];
den = [];
for n = 1 : length(H_s)
if mod(n, 2) == 1
num = [num H_s(1, n)];
den = [den 0];
else
num = [num 0];
den = [den H_s(1, n)];
end
end
num_org = removehead(num);
den_org = removehead(den);
num = num_org;
den = den_org;
for a = 1:2
if a == 1
fprintf('===== Calculate T type %d th ladder for RS = 0 ohm, RS = 1 ohm and Omega = 1 rad/sec =====', order)
fprintf(char(10))
fprintf('From output, arm element is created first.')
fprintf(char(10))
IsTtype = 1;
else
fprintf(char(10))
fprintf('===== Calculate PIE type %d th ladder for RS = 1 ohm, RL = inf ohm and Omega = 1 rad/sec =====', order)
fprintf(char(10))
fprintf('From input, Shunt element is created first.')
fprintf(char(10))
IsTtype = 0;
den = den_org;
num = num_org;
IsTtype = 0;
end
iter = 1;
lc = [];
while(length(den) > 1)
[quo reminder] = deconv(num, den);
lc = [lc quo(1)];
num = den; %◆分母を分子に移動
%----------------------------------------------------
% 余りのうちの高次のべきで係数がゼロの項を削除して、分母の変数に代入
den = removehead(reminder); %◆あまり（分子）を分母に移動
%----------------------------------------------------
fprintf('Number %d element ', iter)
if (mod(iter,2) == 1)
if (IsTtype == 1)
fprintf('is Arm L element of %2.6d H',quo(1))
else
fprintf('is Shunt C element %2.6d F', quo(1))
end
else
if (IsTtype == 1)
fprintf('is Shunt C element %2.6d F', quo(1))
else
fprintf('is Arm L element %2.6d H', quo(1))
end
end
fprintf(char(10))
iter = iter + 1;
end
fprintf('Last element ')
if ((IsTtype == 1) && (mod(order,2) == 1))
fprintf ('at input port is Arm L element of %2.6d H',num(1)/den)
fprintf(char(10))
fprintf ('and source impedance of %2.6d Ohm',num(2)/den)
else
fprintf ('at output port is Shunt C element %2.6d F', num(1)/den)
fprintf(char(10))
fprintf ('and load conductance of %2.6d Siemence (%2.6d Ohm)',num(2)/den, den/num(2))
end
fprintf(char(10)')
end
24. LCフィルタ設計MATLAB/GNU Octave mファイルに戻る
LCフィルタ設計のTOPに戻る
