% lcfilter.mの後半。前半の後ろに挿入し、UTF-8でテキストファイルとして拡張子.mで保存してください。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% フィルタ設計
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i = sqrt(-1);
if !(type == 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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平面での電圧伝達多項式（√電力）の分母になる
% s平面での電圧伝達多項式の共役を得る
mx = [order:-1:0];
mx = (-1 * ones(1, order+1)) .^ mx;
H_s_conj = H_s .* mx;
%--------------------------------------------------------
% 電力伝達多項式にするために、共役同士を掛け算する（もともとあった電力伝達多項式をそのまま用いないのは、もともとの式ではs平面の右側にも根がある式になっているため、実際の伝達多項式としては成立しないので、このようにする）
%--------------------------------------------------------
H_s_square = real(conv(H_s, H_s_conj));
%--------------------------------------------------------
%★★★★★★★★ ここ重要 ★★★★★★★★
%ここまでで得られた答え（H_s_square）は根から多項式を得たものになっているので、実際の電力伝達率ではない。
% そこで電力伝達率 Power_Transfer_Ratioと得られた答えの定数項（DC）を用いて、正しい電力伝達率に変換する（とくにRS ≠ RLのときに重要な考え方）
%--------------------------------------------------------
inverse_dc_att = (1/Power_Transfer_Ratio)/H_s_square(1,end);
H_s_square = H_s_square * inverse_dc_att;
% 不要な計算誤差残差を取り去る
mask = mod([length(H_s_square):-1:1],2);
H_s_square = H_s_square .* mask;
%--------------------------------------------------------
% これによりH_s_squareが、RS = RLのときの電力伝達率を基準の1Wとしたときの電力伝達多項式として得られる
%--------------------------------------------------------
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Elliptic フィルタ設計
% 以下で必要となる変数: elpnum, elpden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H_s = poly(elpden);
% s平面での電圧伝達多項式（分母）の共役を得る
mx = [order:-1:0];
mx = (-1 * ones(1, order+1)) .^ mx;
H_s_conj = H_s .* mx;
H_s_square = real(conv(H_s, H_s_conj));
H_s_num = poly(elpnum);
% s平面での電圧伝達多項式（分子）の共役を得る
order_num = length(H_s_num)-1;
mx = [order_num:-1:0];
mx = (-1 * ones(1, order_num+1)) .^ mx;
H_s_num_conj = H_s_num .* mx; %奇数次項はゼロなのでやってもやらなくても同じ
H_s_num_square = real(conv(H_s_num, H_s_num_conj));
% H_s_num_squareの長さをH_s_squareにそろえる
if (length(H_s_square) == length(H_s_num_square))
vectone = [];
else
vectone = zeros(1, length(H_s_square)-length(H_s_num_square));
end
H_s_num_square = [vectone H_s_num_square];
inverse_dc_att = (1/Power_Transfer_Ratio) * H_s_num_square(1,end)/H_s_square(1,end);
H_s_square = H_s_square * inverse_dc_att;
% 不要な計算誤差残差を取り去る
mask = mod([length(H_s_square):-1:1],2);
H_s_square = H_s_square .* mask;
%--------------------------------------------------------
% H_s_num_square
% ここまでで電力伝達関数が ---------------- として得られた
% H_s_square
%--------------------------------------------------------
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 正規化された反射電力（反射係数Γの自乗）を計算する（Power_Transfer_Ratioの係数がかかっているので、実際の正規化された反射電力になっている）
% ref_sq_num 1
%Γ^2 = ----------------- = 1 - -----------------
% ref_sq_denom H_s_square
%
% Elliptic フィルタの場合
% ref_sq_num H_s_num_square
%Γ^2 = ----------------- = 1 - -----------------
% ref_sq_denom H_s_square
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if !(type == 4) % Butterworth, Chebyshev, Besselのケース
H_s_num_square = [zeros(1, length(H_s_square)-1) 1];
else % Elliptic case
H_s_num_square = H_s_num_square;
end
%Γ^2の分子を計算
ref_sq_num = H_s_square - H_s_num_square;
% 不要な計算誤差残差を取り去る
ref_sq_num = ref_sq_num .* mask;
%Γ^2の分母
ref_sq_denom = H_s_square;
%--------------------------------------------------------
%以降ではPower_Transfer_Ratioは使用していない（以降の同様の説明部分を参照のこと）
%--------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 電力伝達多項式を用いてプロットする
close all
if PLOT1 % 1 にした場合にプロットし、0だとプロットしない
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = i * logspace(-2, 3, 10001); % s = jωとする
accum_num = polyval(H_s_num_square, s);
accum_denom = polyval(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
%--------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 反射電力多項式を用いてプロットする
if PLOT2 % 1 にした場合にプロットし、0だとプロットしない
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = i * logspace(-2, 3, 10001); % s = jωとする
accum_num = polyval(ref_sq_num, s);
accum_denom = polyval(ref_sq_denom, s);
figure
plotacccum = 10*log10(abs(accum_num)) - 10*log10(abs(accum_denom));
semilogx(abs(s), plotacccum)
xlabel('Normalized Angular Frequency')
ylabel('Reflection Level [dB]')
set(gca, 'FontSize', 18);
botofY = floor(plotacccum(1,1)) - 2;
axis([0.01 100 botofY 2])
grid on
end
%--------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 第2ステップ
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------------------------------------
% ここまでで、反射電力（反射係数Γの自乗）の多項式が得られた。
%
% ref_sq_num 1
%Γ^2 = ----------------- = 1 - -----------------
% ref_sq_denom H_s_square
%
% Elliptic フィルタの場合
% ref_sq_num H_s_num_square
%Γ^2 = ----------------- = 1 - -----------------
% ref_sq_denom H_s_square
%
% このうちの分母、分子それぞれ半数の根をとる（ルートをとることとほぼ同義）ことにより、
% 実際の反射係数Γの多項式を得る
%
% reflectpolynum
% Γ = ----------------------------
% reflectpolydenom
% として計算する。
%--------------------------------------------------------
rootsofnum = roots(ref_sq_num);
rootsofnumforsort = [real(rootsofnum) rootsofnum];
rootsofnumforsort = sortrows(rootsofnumforsort, 1);
rootsofnum = rootsofnumforsort(1:order, 2); % s平面の左半分の根を得る
%=================================================
%
% 半分の根から、反射係数Γの分子多項式を計算する
reflectpolynum = real(poly(rootsofnum));
%=================================================
% 反射係数Γの分母多項式を計算する（H_sと同じになるので、そのまま代入）
reflectpolydenom = H_s;
%=================================================
% これにより反射係数Γの有利関数多項式が以下に得られる
%
% reflectpolynum
% Γ = ----------------------------
% reflectpolydenom
%
%
%★★★★★★★★ ここ重要 ★★★★★★★★
% ここで根から多項式に変換した結果では、Power_Transfer_Ratioは使用していない
%「根から多項式をつくると最上位が1になりPower_Transfer_Ratioが必要なのでは？」
% 「分母はH_sをそのまま使っているけれども、合っている？」という疑問が出てくる
% この理由は反射電力（反射係数Γの自乗）の式、ref_sq_numとref_sq_denomの最高次
% の項の係数が同じ大きさになっているため（分子が1 - H_s_squereであり、分母
% がH_s_squereで、定数項のみが引き算されるから）、Γ^2の分母分子の比と、Γ
% の分母分子の比、それぞれの比は変わらない。
% 反射係数ということで、比なので係数Power_Transfer_Ratioは関係ない。
% また分子であるref_sq_numの根から計算された反射係数分子多項式reflectpolynum
% と、reflectpolydenomであるH_sも、H_s = poly(rootsofpower)で根から計算された
% 多項式、それぞれの最高次の係数は同じになるため、全体として辻褄があう
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT 3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%反射係数多項式を用いてプロットする（正しく動作しているかの確認のため）
if PLOT3 % 1 にした場合にプロットし、0だとプロットしない
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = i * logspace(-2, 3, 10001); % s = jωとする
accum_num = polyval(reflectpolynum, s);
accum_denom = polyval(reflectpolydenom, s);
figure
plotacccum = 20*log10(abs(accum_num)) - 20*log10(abs(accum_denom));
semilogx(abs(s), plotacccum)
xlabel('Normalized Angular Frequency')
ylabel('Reflection Level [dB]')
set(gca, 'FontSize', 18);
botofY = floor(plotacccum(1,1)) - 2;
% axis([0.01 100 botofY 2])
axis([0.01 100 -0.003 0.001])
grid on
end
%--------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ここまでの処理で得られた反射係数有利関数多項式
%
% reflectpolynum
% Γ = ----------------------------
% reflectpolydenom
%
% つづいて入力端からみたインピーダンスの有利関数多項式を
% 1 +/- reflectpolynum/reflectpolydenom
% Zin = Rs --------------------------------------
% 1 -/+ reflectpolynum/reflectpolydenom
% ただしRs = 1
% で計算し、Butterworth, Chebyshev, Besselのケースでは連分数を計算することで、
% またEllipticのケースでは別の計算方法を用いることで、それぞれの素子定数を得る
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ラダー素子定数を計算する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if type == 4
% Ellipticの場合、ゼロ周波数を計算する
elpnum = imag(elpnum);
elpnum = sort(elpnum);
elpnum = i * elpnum(floor(order/2+1):end,1);
end
for a = [1 2]
switch(a)
case 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1回目のループでは、入力インピーダンス有利関数多項式ケース1 (＋Γ)を計算
% 1 + reflectpolynum/reflectpolydenom num
% Zin = Rs -------------------------------------- = -----
% 1 - reflectpolynum/reflectpolydenom den
% ただし Rs = 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num = reflectpolydenom + reflectpolynum;
den = reflectpolydenom - reflectpolynum;
case 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2回目のループでは、入力インピーダンス有利関数多項式ケース2 (ーΓ)を計算
% 1 - reflectpolynum/reflectpolydenom num
% Zin = Rs -------------------------------------- = -----
% 1 + reflectpolynum/reflectpolydenom den
% ただし Rs = 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num = reflectpolydenom - reflectpolynum;
den = reflectpolydenom + reflectpolynum;
end
%----------------------------------------------------
% 分子の高次のべきで係数がゼロの項を削除
num = removehead(num);
% 分母の高次のべきで係数がゼロの項を削除
den = removehead(den);
%----------------------------------------------------
if 0 % Print impedance equation or not
fprintf(char(10))
fprintf('vvvvvvvvvvvvvvvvvvvv Input impedance equation vvvvvvvvvvvvvvvvvvvv')
fprintf(char(10))
num
den
fprintf('^^^^^^^^^^^^^^^^^^^^ END ^^^^^^^^^^^^^^^^^^^^')
fprintf(char(10))
end
if (length(num) > length(den))
fprintf(char(10))
fprintf('===== Calculate T type %d th ladder for RS = 1 ohm and Omega = 1 rad/sec =====', order)
fprintf(char(10))
fprintf('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 and Omega = 1 rad/sec =====', order)
fprintf(char(10))
fprintf('Shunt element is created first')
fprintf(char(10))
IsTtype = 0;
temp = den; % 分母と分子を交換（逆数にする）
den = num;
num = temp;
end
if !(type == 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Butterworth, Chebyshev and Besselフィルタの素子定数
% Zin = num/den
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iter = 1;
while(length(den) > 1)
[quo reminder] = deconv(num, den);
num = den; %◆分母を分子に移動
%----------------------------------------------------
% 余りのうちの高次のべきで係数がゼロの項を削除して、分母の変数に代入
reminder = removehead(reminder);
den = 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 ('is Arm L element of %2.6d H',num(1)/den)
fprintf(char(10))
fprintf ('and load impedance of %2.6d Ohm',num(2)/den)
else
fprintf ('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))
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ellipticフィルタの素子定数を計算する
% Zin = num/den
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iter = 1;
while(length(den) > 1) %分母の次数が0次のみになるまで実行
if (mod(iter,2) == 1)
% -----------------------------------------------
% 単一素子の計算
% T型を例にしてインピーダンスで示しているが、π型の場合はアドミッタンスになるので、ZをYで読み替える
% 連分数展開でできる他のフィルタより、このEllipticのケースのほうが厄介だった
% 共振部分 Zres(Yも同じ) = sL + 1/sC = (s^2 LC + 1)/sC = (s^2 + ω^2)/(s/L)
% があるので、Z1stにs = ω_res = s_zero（ゼロ点による直列共振）を代入して計算したリアクタンス（アドミッタンス）が最初の単一素子の大きさになる
% -----------------------------------------------
% Ellipticのゼロ周波数を代入。ゼロ点であるelpnumは低い周波数から並んでいる。低い周波数から計算
if (floor(iter/2)+1 <= length(elpnum))
omega = elpnum(floor(iter/2)+1, 1);
else
omega = 1;
% 最後のエレメントを計算するときオーバになってしまうので、1を代入しておく
end
% -----------------------------------------------
% L = Z1st(s)にs = ω_res = s_zeroを代入してLを求める（ちゃんと都合よく実数部のみが残る）。
% なおπ型の場合はCが主要素子として求まる
% -----------------------------------------------
% 分子にs = s_zeroを代入して計算
Ztemp_num_calc = obtainZYvalue(num, omega);
% 分母にs = s_zeroを代入して計算
Ztemp_den_calc = obtainZYvalue(den, omega);
SingleElement = Ztemp_num_calc/(Ztemp_den_calc * omega);
% 余計な計算誤差を取り去る
SingleElement = real(SingleElement);
%単一素子SingleElementを除外した残りの共振回路から負荷を見たインピーダンス（もしくはアドミッタンス）
% Zright_remain = Z1st - s * SingleElement
tmp = conv([SingleElement 0], den);
num = num - tmp;
% -----------------------------------------------
% Zright_remain = num/den
% -----------------------------------------------
else
% -----------------------------------------------
% 共振素子の計算（ゼロ）。
% T型を例にしてインピーダンスで示しているが、π型の場合はアドミッタンスになるので、ZをYで読み替える
% 共振部分 Zres(Yも同じ) = sL + 1/sC = (s^2 LC + 1)/sC = (s^2 + ω^2)/(s/L)
% ※ 一応、π型の並列共振も示しておく。Yosc = sC + 1/sL = (s^2 LC + 1)/sL = (s^2 + ω^2)/(s/C)で、π型のアーム並列共振部はCが主要パラメータになる
%
% ここでω^2 = 1/LC。【重要】このωは伝達関数の分子のゼロ点で位置は定数として決まるが、ここのωは虚数ではなく実数
% あくまでも「ω^2 = 1/LC」であり、虚数として考えていない（あとでs = s_zeroと代入するが、そこは虚数）
% 1/Zright_remain = 1/Zres + 1/Znext となるので、1/Znext = 1/Zright_remain - 1/Zres
% Zright_remainは、(s^2 + ω^2)で直列共振するため、(s^2 + ω^2)で分母を因数分解できる。
% 残りのインピーダンスを素子形状がどうであれ、その部分をZtempとすれば、
% Zright_remain = [Ztemp(s^2 + ω^2)]と因数分解できる
%（Ztemp = Zright_remain/(s^2 + ω^2)となり、インピーダンスの分母がゼロになることで直列共振）。
% これにより
% 1/Znext = 1/Zright_remain - 1/Zres = 1/[Ztemp(s^2 + ω^2)] - (s/L)/(s^2 + ω^2)
% 両辺に(s^2 + ω^2)を掛けると
% (s^2 + ω^2)/Znext = 1/Ztemp - (s/L)となるので、
% ① Zright_remain = Ztemp(s^2 + ω^2)から、Ztemp = Zright_remain/(s^2 + ω^2)を計算し、
% sに-jω = s_zero（ω, s_zeroはすでにゼロ点で最初から求まっているので）を代入することで、0 = 1/Ztemp - (s/L)から、
% ② s = s_zero（虚数）でL = s_zero ZtempでLを求めることができる
% ③ Cはω^2 = 1/LCよりC = 1/Lω^2
%これで共振部分 Zres(Yも同じ)の定数が求まる
%
% ④ 1/Znext = 1/Zright_remain - 1/ZresでZnextを計算すればよい
% なお(s^2 + ω^2)/Znext = 1/Ztemp - (s/L)を使える。
% 右辺を変換し、逆数をとり、(s^2 + ω^2)で割れば共振部分を除外した、残りの素子部分のインピーダンスZnextを計算できる
% -----------------------------------------------
% Ellipticのゼロ周波数を代入。ゼロ点であるelpnumは低い周波数から並んでいる。さきに低い周波数から計算
omega = elpnum(round(iter/2), 1);
% -----------------------------------------------
% ステップ① Ztemp = Zright_remain/(s^2 + ω^2)を割り算で計算する（ZでもYでも同じ計算になる）。
% 実際はZright_remainの分子を(s^2 + ω^2)で割り、Ztemp_numを得る
% Zright_remain = num/den
% Ztemp = Zright_remain/(s^2 + ω^2) = Ztemp_num/den
% -----------------------------------------------
[Ztemp_num reminder] = deconv(num, [1 0 abs(omega^2)]);
% reminderに計算誤差が少し残るが使わない
% -----------------------------------------------
% ステップ② L = s x Ztemp(s = omega)でLを求める（複素数になるが、ちゃんと都合よく実数部のみが残る）。なおπ型の場合はCが主要素子として求まる
% -----------------------------------------------
% 分子にs = wを代入して計算
Ztemp_num_calc = obtainZYvalue(Ztemp_num, omega);
% 分母にs = wを代入して計算
Ztemp_den_calc = obtainZYvalue(den, omega);
PrimaryElement = omega * Ztemp_num_calc/Ztemp_den_calc;
PrimaryElement = real(PrimaryElement); %虚数部の計算誤差残差を削除
% -----------------------------------------------
% ステップ③ CをC = 1/Lω^2で計算する。なおπ型の場合はLが二次素子として求まる
% -----------------------------------------------
SecondalyElement = 1/(PrimaryElement*abs(omega^2));
% -----------------------------------------------
% ④ F = 1/Znext = 1/Zright_remain - 1/ZresでZnextを計算する
% T型を例にしてインピーダンスで示しているが、π型の場合はアドミッタンスになるので、ZをYで読み替える
% (s^2 + ω^2)/Znext = 1/Ztemp - (s/L), L = PrimaryElementを使う
% Ztemp_num
% Ztemp = ------------
% den
%
% Fnum s^2 + ω^2 den s
% F = ----- = --------- = --------- - ---
% Fden Znext Ztemp_num L
%として計算していく
% -----------------------------------------------
Fden = PrimaryElement * Ztemp_num; % F = Fnum/Fdenの分母
Fnum_1 = PrimaryElement * den; % F = Fnum/Fdenの分子第1項
Fnum_2 = conv([1 0], Ztemp_num); % F = Fnum/Fdenの分子第2項
% ベクトル計算のため、分子1項、2項の多項式ベクトルの長さをそろえる
len_Fnum_1 = length(Fnum_1);
len_Fnum_2 = length(Fnum_2);
if len_Fnum_1 > len_Fnum_2
Fnum_2 = [zeros(1,(len_Fnum_1 - len_Fnum_2)) Fnum_2];
else
if len_Fnum_1 < len_Fnum_2
Fnum_1 = [zeros(1,(len_Fnum_2 - len_Fnum_1)) Fnum_1];
end
end
Fnum = Fnum_1 - Fnum_2;
% -----------------------------
% s^2 + ω^2 Fnum
% --------- = F = -----
% Znext Fden
% が得られた
% Znext = (s^2 + ω^2)Fden/Fnum
% を計算する
% -----------------------------
num = conv(Fden, [1 0 abs(omega^2)]);
den = Fnum;
% -----------------------------
% den
% Znext = -----
% num
% が得られた
DC_impedance = num(1,end)/den(1,end);
% 因数分解して分母分子で共通になる項を取り去る
% 分子多項式の方が次数が高い
rootsofnum = roots(num);
rootsofden = roots(den);
rootsofnum_remain = rootsofnum;
rootsofden_remain = rootsofden;
% 分子の根と分母の根が同じな場合、その分母の根をゼロにする
% とりあえず誤差（差異）を0.01%とした。多項式が高次になるほど誤差が増える。11次ではだいぶ合わない
for p = 1:length(rootsofnum)
for q = 1:length(rootsofden)
if abs(rootsofnum(p) - rootsofden(q)) < 1e-4
rootsofnum_remain(p) = 0;
rootsofden_remain(q) = 0;
end
end
end
% ゼロになった根を取り去る
tmp = [];
for n = 1:length(rootsofnum_remain)
if !(rootsofnum_remain(n) == 0)
tmp = [tmp; rootsofnum_remain(n)];
end
end
rootsofnum_remain = tmp;
tmp = [];
for n = 1:length(rootsofden_remain)
if !(rootsofden_remain(n) == 0)
tmp = [tmp; rootsofden_remain(n)];
end
end
rootsofden_remain = tmp;
% 再度分母・分子を有理関数にする
num = poly(rootsofnum_remain);
den = poly(rootsofden_remain);
% 根から多項式を得たため、DC項を補正
DC_compensate = DC_impedance* den(1,end)/num(1,end);
num = DC_compensate * num;
end
fprintf('Number %d element ', iter)
if (mod(iter,2) == 1)
if (IsTtype == 1)
fprintf('is Arm L element of %2.6d H',SingleElement)
else
fprintf('is Shunt C element %2.6d F', SingleElement)
end
else
if (IsTtype == 1)
if !(type == 4) %Ellipse以外
fprintf('is Shunt C element %2.6d F', quo(1))
else %Ellipse
fprintf('is Shunt Series L %2.6d H and C %2.6d F',PrimaryElement, SecondalyElement)
end
else
if !(type == 4) %Ellipse以外
fprintf('is Arm L element %2.6d H', quo(1))
else %Ellipse
fprintf('is Arm Parallel C %2.6d F and L %2.6d H',PrimaryElement, SecondalyElement)
end
end
end
fprintf(char(10))
iter = iter + 1;
end
fprintf('Last element ')
if ((IsTtype == 1) && (mod(order,2) == 1))
fprintf ('is Arm L element of %2.6d H',num(1)/den)
fprintf(char(10))
fprintf ('and load impedance of %2.6d Ohm',num(2)/den)
else
fprintf ('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
fprintf(char(10)')
fprintf(char(10)')
end
24. LCフィルタ設計MATLAB/GNU Octave mファイルに戻る
LCフィルタ設計のTOPに戻る
