ドルフ・チェビシェフ窓の時間特性・周波数特性

N=9（M=4）の場合、 v_k(n)は以下のようになる。

k/n -4(=-M) -3 -2 -1 0 1 2 3 4(=M)
0 0 0 0 0 1 0 0 0 0
1 0 0 0 γ/2 γ-1 γ/2 0 0 0
2 0 0 γ^2/2 2γ^2 2(γ^2-2γ+1) 2γ^2 γ^2/2 0 0
3 0 □ □ □ □ □ □ □ 0
4=M □ □ □ □ □ □ □ □ □ 
clear; close all;

N = 129;
% N = 9;
M = (N-1)/2;
T = 1;
% beta = 1;%方形窓の場合1
beta = 3;%方形窓の場合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

%ホントは以下が正しいと思われるが端点の変化が急峻となる。
figure();
w10 = w10(2:end-1);
ns = ns(2:end-1);
plot(ns, w10);
title('ドルフ・チェビシェフ窓（-M≤n≤M）');

%端点を削除したバージョン（ホントはだめだが傾向を見たいだけなので...）
figure();
w10 = w10(3:end-2);
ns = ns(3:end-2);
plot(ns,w10);
title('ドルフ・チェビシェフ窓（-(M-1)≤n≤(M-1)）');
big;


端点の計算がおかしい。たぶん漸化式の最後のM行目の計算で(n-1),(n+1)を使っているのが問題と思われる。

端点を削除すると以下になる。

M行目の計算で、そもそも(n-1),(n+1)の点を使わないようにすると以下になる。急激な変化はおこらなくなるが、なんか違う気がする。

βがパラメータとなっているので、この値を１から上げてみたら端点で滑らかに接続するようになった。以下はβ=3の場合

周波数特性

clear; close all; 

N=127;
% % N = 9;
M = (N-1)/2;
T = 1;
% beta = 1;%方形窓の場合1
beta = 3;%方形窓の場合1
gamma = (1+cos(pi/(2*M)))/(1+cos(2*beta*pi/(2*M+1)));

%チェビシェフ多項式
syms w;
chebyshev = @(c1,c0) 2*w*c1 - c0;
c0 = 1;
c1 = w;
cs = [c0 c1];
for n=2:M
    n
  cn_prev = cs(n-1);
  cn = cs(n);
  cn_next = chebyshev(cn, cn_prev);
  cs(end+1) = simplify(cn_next);
end

% epsilon = 0.0001;
epsilon = 0.0001;
omega = linspace(0, pi-epsilon, 200);
T = 1;
CM = cs(end);
A = gamma*cos(omega*T+gamma-1);
B = 2*gamma-1;
Denom = [];
for idx = 1:length(A)
    Denom(end+1) = subs(CM, w, A(idx));
end
Num = subs(CM,w,round(B,4));
W10 = eval(Denom/Num);
lab = sprintf('ドルフ・チェビシェフ窓の振幅応答(N=%d)',N);
plot(omega, 20*log10(abs(W10)/max(abs(W10))));
title(lab);
ylabel('20 log10 |W|');
xlabel('ωT');

big;