clear; close all; global delta_p delta_s; delta_p = 1; delta_s = 1; global wc; wc = 0.5*pi; wd = 0.6*pi; N = 31; M = (N-1)/2; %STEP1: ωiを初期化 %ωiを初期化: 当分割 epsilon = 0.001; if mod(M,2)==1 %M=奇数 wi_left = linspace(0+epsilon,wc,(M+1)/2+1); wi_right = linspace(wd,pi-epsilon, (M+1)/2); else %M=偶数 wi_left = linspace(0,wc,(M+2)/2); wi_right = linspace(wd,pi, (M+2)/2); end wi = [wi_left, wi_right]; %STEP2: Qを計算 num = 0; den = 0; for i=0:M+1 alpha_i = 1; for k=0:M+1 if k==i continue else alpha_i = alpha_i * 1/(cos(wi(k+1))-cos(wi(i+1))); end end num = num + alpha_i * D(wi(i+1)); den = den + ((-1)^i * alpha_i) / W(wi(i+1)); end Q = num / den; %STEP3: H(ωi)を計算し、さらにH(ω)を計算する。 %H(ωi)の計算 Hwi = zeros(size(wi)); for i=0:M+1 Hwi(i+1) = D(wi(i+1)) - (-1)^(i)*Q/W(wi(i+1)); end figure(); plot(wi, Hwi,'o-'); title('H(wi)'); %H(ω)の計算 syms w; H_num = 0; for i=0:M+1 gamma_i = 1; for k=0:M if k==i continue else gamma_i = gamma_i * 1/(cos(wi(k+1)-cos(wi(i+1)))); end end H_num = H_num + (gamma_i / (cos(w)-cos(wi(i+1)))) * Hwi(i+1); end H_den = 0; for i=0:M+1 gamma_i = 1; for k=0:M if k==i continue else gamma_i = gamma_i * 1/(cos(wi(k+1)-cos(wi(i+1)))); end end H_den = H_den + gamma_i / (cos(w)-cos(wi(i+1))); end H = H_num / H_den; %STEP4: E=W•(H-D)の計算 omega = linspace(0,pi,500); E = zeros(size(omega)); for i=1:length(omega) val = W(omega(i)) * (subs(H,w,omega(i))-D(omega(i))); E(i) = vpa(val); end %STEP5: |E|と|Q|の比較 if abs(E)<= Q return else disp('Continue Parks-McCllalen Algorithm') end figure(); plot(omega, E); [E_pks,wi_locs] = findpeaks(E); title('E'); %STEP6: E(w)の極値を与えるwiとその個数を見出す wi2 = omega(wi_locs); L = length(wi2); %STEP7: 極値の個数Lによって場合分けを行う if L==(M+2) disp('L=M+2. Proceed to Step8'); elseif L==(M+3) disp('L=M+3. Proceed to Step8'); else lab = sprintf('L=%d. M=%d. M+2=%d. M+3=%d. Computational error',L,M,M+2,M+3); disp(lab); return; end %STEP8: Qの再計算をする num = 0; den = 0; for i=0:M+1 alpha_i = 1; for k=0:M+1 if k==i continue else alpha_i = alpha_i * 1/(cos(wi2(k+1))-cos(wi2(i+1))); end end num = num + alpha_i * D(wi2(i+1)); den = den + ((-1)^i * alpha_i) / W(wi2(i+1)); end Q2 = num / den; %STEP9: Qに変化があったかを判定する if abs(Q-Q2)>1e-5 disp('Q changed.'); else disp('NO change in Q.'); return end %STEP10: fiの計算 %fiを求める式 Σ_k fk cos(k•wi) = G(wi), i=0 % i=0: f0•cos(0•w0) + f1•cos(1•w0) + ... + fM•cos(M•w0) = G(w0) % i=1: f0•cos(0•w1) + f1•cos(1•w1) + ... + fM•cos(M•w1) = G(w1) % i=M: f0•cos(0•wM) + f1*cos(1•wM) + ... + fM•cos(M•wM) = G(wM) % i=M+1: f0•cos(0•w_M+1) + f1*cos(1•w_M+1) + ... + fM•cos(M•w_M+1) = G(wM+1) % これらは以下の連立方程式に変形できる %(M+2)x(M+1) 行列 (M+1)x1 ベクトル (M+2)x1 ベクトル % [cos(0•w0) cos(1•w0) ... cos(M•w0)] [ f0 ] [ G(w0) ] % [cos(0•w1) cos(1•w1) ... cos(M•w1)] [ f1 ] [ G(w1) ] % [ ] [ ] = [ ] % [cos(0•wM) cos(1•wM) ... cos(M•wM)] [ fM ] [ G(wM) ] % [cos(0•wM+1) cos(1•wM+1) ... cos(M•wM+1)] [ G(wM+1) ] % C • f = Gとすると、 % f = inv(C)•Gで代数的に求まる。 C = zeros(M+2, M+1); for k=0:M+1 for i=0:M C(k+1,i+1) = cos(k*wi2(i+1)); end end G = zeros(M+2, 1); for i=0:M+1 G(i+1) = D(wi2(i+1)) - (-1)^i*Q2/W(wi2(i+1)); end f = C\G; h_head = flip(f(2:end) / 2); h_tail = f(1); h = [h_head; h_tail];%縦ベクトル c = w*zeros(M,1); for k=1:M c(k) = 2*cos(k*w); end c = flip(c); c(end+1) = 1; H_final = transpose(c) * h; H_vals = vpa(subs(H_final, w, omega)); plot(omega,H_vals); big; function res = D(w) global wc; if w<wc res=1; else res=0; end end function res = W(w) global delta_p delta_s; global wc; if w<wc res=1; else res=delta_p/delta_s; end end
計算は最後まで到達したが、ローパスフィルタは得られていない。
Matlabにはfirpm()関数が実装されているので、これを使うのがよいのか。
