以下の３つのコードをマージする。
（理想ディジタルフィルタ（LPF）、カイザー窓関数、ドルフ・チェビシェフ窓関数）
seinzumtode.hatenadiary.jp
seinzumtode.hatenadiary.jp
seinzumtode.hatenadiary.jp
clear; close all; wc = pi/2; T = 1; % N=61; N=121; l = (N-1)/2; hk = []; w0 = pi/8; ks = 1:N; M = (N-1)/2; %インパルス応答hkを計算 for kprime = ks k = kprime-1; hk1 = wc*T/pi*sin((k-l)*wc*T)/(k-l)/wc/T; hk(end+1) = hk1; end hk(l+1) = wc*T/pi; % rho = 40;%dB rho = 120;%dB a = 0.5842*(rho-21)^(0.4)+0.07886*(rho-21); if rho>50 a = 0.1102*(rho-8.7); end %窓関数を計算（カイザー窓） rectangular = 0; hamming = 5.4414; blackmann = 8.885; kaiser = a; n = -M:M; A =@(alpha) alpha * sqrt(1-(2*n/(N-1)).^2); Den =@(alpha) besseli(0, A(alpha)); Num =@(alpha) besseli(0,alpha); w9_r = Den(rectangular)/Num(rectangular);%方形窓の時間特性 w9_h = Den(hamming)/Num(hamming);%ハミング窓の時間特性 w9_b = Den(blackmann)/Num(blackmann);%ブラックマン窓の時間特性 w9_k = Den(kaiser)/Num(kaiser);%カイザー窓の時間特性 %窓関数を計算(ドルフ・チェビシェフ窓) beta = 5;%方形窓の場合1 gamma = (1+cos(pi/(2*M)))/(1+cos(2*beta*pi/(2*M+1))); ns = -M-1:M+1; v_k_n = zeros(M+1,2*M+3); %v0nとv1nを初期化 for k=0:1 if k==0 %v0nの初期化 for n=-M:M if n==0 v_k_n(k+1,find(ns==n)) = 1; end end elseif k==1 %v1nの初期化 for n=-M:M if n==0 v_k_n(k+1,find(ns==n)) = gamma-1; elseif abs(n)==1 v_k_n(k+1, find(ns==n)) = gamma/2; end end end end %vM(n)の計算: M-1まで for k=2:M-1 for n=-k:k v_k_n(k+1, find(ns==n)) = ... +2*(gamma-1)*v_k_n(k, find(ns==n)) ... -v_k_n(k-1, find(ns==n)) ... +gamma*(v_k_n(k, find(ns==(n-1)))+v_k_n(k, find(ns==(n+1)))); end end %vM(M)の計算: 端点の処理を追加 for k=M for n=-k:k if abs(n)==k v_k_n(k+1, find(ns==n)) = ... +2*(gamma-1)*v_k_n(k, find(ns==n)) ... -v_k_n(k-1, find(ns==n)); else v_k_n(k+1, find(ns==n)) = ... +2*(gamma-1)*v_k_n(k, find(ns==n)) ... -v_k_n(k-1, find(ns==n)) ... +gamma*(v_k_n(k, find(ns==(n-1)))+v_k_n(k, find(ns==(n+1)))); end end end %w10(n)の計算 w10 = zeros(size(v_k_n(M+1,:))); for n=-M:M; vM_n = v_k_n(M+1,find(ns==n)); vM_0 = v_k_n(M+1, find(ns==0)); w10(find(ns==n)) = vM_n/vM_0; end w10 = w10(2:end-1); %hkをZ変換して伝達関数を計算 Hz_r = 0; Hz_h = 0; Hz_b = 0; Hz_k = 0; Hz_d = 0; syms z; for n = 0:(length(hk) - 1) Hz_r = Hz_r + w9_r(n+1)*hk(n+1) * z^(-n); Hz_h = Hz_h + w9_h(n+1)*hk(n+1) * z^(-n); Hz_b = Hz_b + w9_b(n+1)*hk(n+1) * z^(-n); Hz_k = Hz_k + w9_k(n+1)*hk(n+1) * z^(-n); Hz_d = Hz_d + w10(n+1)*hk(n+1) * z^(-n); end Hz_r = vpa(simplify(Hz_r)); Hz_h = vpa(simplify(Hz_h)); Hz_b = vpa(simplify(Hz_b)); Hz_k = vpa(simplify(Hz_k)); Hz_d = vpa(simplify(Hz_d)); syms w; %伝達関数H(z)にz=exp(jωT)を代入して周波数応答関数を計算 Hw_r = subs(Hz_r, z, exp(i*w*T)); Hw_h = subs(Hz_h, z, exp(i*w*T)); Hw_b = subs(Hz_b, z, exp(i*w*T)); Hw_k = subs(Hz_k, z, exp(i*w*T)); Hw_d = subs(Hz_d, z, exp(i*w*T)); omegas = linspace(0,pi,100); lab_r = sprintf('Rectangular(N=%d)',N); lab_h = sprintf('Hamming(N=%d)',N); lab_b = sprintf('Blackmann(N=%d)',N); lab_k = sprintf('Kaiser(N=%d)',N); lab_d = sprintf('Dolph-Chebyshev(N=%d)',N); Hw_abs_r = abs(subs(Hw_r, w, omegas)); Hw_abs_h = abs(subs(Hw_h, w, omegas)); Hw_abs_b = abs(subs(Hw_b, w, omegas)); Hw_abs_k = abs(subs(Hw_k, w, omegas)); Hw_abs_d = abs(subs(Hw_d, w, omegas)); plot(omegas, Hw_abs_r,'DisplayName',lab_r); hold on; plot(omegas, Hw_abs_h,'DisplayName',lab_h); plot(omegas, Hw_abs_b,'DisplayName',lab_b); plot(omegas, Hw_abs_k,'DisplayName',lab_k); plot(omegas, Hw_abs_d,'DisplayName',lab_d); legend(); title('振幅特性'); figure(); plot(omegas, 20*log10(Hw_abs_r/max(Hw_abs_r)),'DisplayName',lab_r); hold on; plot(omegas, 20*log10(Hw_abs_h/max(Hw_abs_h)),'DisplayName',lab_h); plot(omegas, 20*log10(Hw_abs_b/max(Hw_abs_b)),'DisplayName',lab_b); plot(omegas, 20*log10(Hw_abs_k/max(Hw_abs_k)),'DisplayName',lab_k); plot(omegas, 20*log10(Hw_abs_d/max(Hw_abs_d)),'DisplayName',lab_d); legend(); title('振幅特性（dB表示）'); ylabel('20 log10 |H(ω)|'); big;
