function Hk_z = oddHkz(Theta,M_ch,m)

% M_ch: 
% m:
% Theta:

Num_lattice = floor(M_ch/2);                                           % 無損網格個數
Theta = Theta(1:m*Num_lattice);
Theta = reshape(Theta,Num_lattice,[]);

[iRow,iCol] = size(Theta);
if((iRow ~= Num_lattice) || (iCol ~= m))
    disp('Error Theta');
    return;
end

syms w z;
Theta = sym(Theta);
pi = 4*atan(1);
P_puredelay = [pi/4 pi/2*ones(1,m-1)];    % 純延時單元的參數即為初始化的參數
P_puredelay = sym(P_puredelay);

% 第1步：由Theta得到G_k(z)和G_Mk(z)
G_kpz = [cos(Theta(:,1)) zeros(Num_lattice,m-1)];
G_Mkpz = [sin(Theta(:,1)) zeros(Num_lattice,m-1)];
G_oddk = [cos(P_puredelay(1)) zeros(1,m-1)];
G_oddMk = [sin(P_puredelay(1)) zeros(1,m-1)];

j = sqrt(-1);
z = -exp(2*M_ch*j*w);

for k=1:Num_lattice
    for p=1:(m-1)
        G_kpz(k,p+1) = G_kpz(k,p)*cos(Theta(k,p+1)) + z^(-1)*G_Mkpz(k,p)*sin(Theta(k,p+1));
        G_Mkpz(k,p+1) = G_kpz(k,p)*sin(Theta(k,p+1)) - z^(-1)*G_Mkpz(k,p)*cos(Theta(k,p+1));
    end
end
% 以上共有Num_lattice個網格，當M_ch為奇數時，還需要添加一個單獨的網格，用來表示剩余的兩個純延時單元
for p=1:(m-1)
    G_oddk(p+1) = G_oddk(p)*cos(P_puredelay(p+1)) + z^(-1)*G_oddMk(p)*sin(P_puredelay(p+1));
    G_oddMk(p+1) = G_oddk(p)*sin(P_puredelay(p+1)) - z^(-1)*G_oddMk(p)*cos(P_puredelay(p+1));
end

% 第2步：根據對稱性得出所有的2M個G(z)值
% G_k(z) = z^(-(m-1))G_2Mk(z)'   0 =< k =< M-1 --- algorithm
% G_z = [G_kpz(:,m);zeros(Num_lattice+1,1);G_Mkpz(:,m);zeros(Num_lattice+1,1)];
% for k=1:2*M_ch
%     if((k > Num_lattice+1 && k <= M_ch+1) || (k > M_ch+1 && k <= 2*M_ch))
%         G_z(k) = z^(-(m-1))*G_z(2*M_ch-k+1)';
%     else
%         G_z(k) = G_z(k);
%     end
% end
% modified 20071015 ------------------ %
G_z = [G_kpz(:,m);zeros(Num_lattice+1,1);G_Mkpz(:,m);zeros(Num_lattice+1,1)];
for k=1:2*M_ch
    if((k > Num_lattice+1 && k <= M_ch) || (k > M_ch+Num_lattice+1 && k <= 2*M_ch))
        G_z(k) = z^(-(m-1))*G_z(2*M_ch-k+1)';
    else
        G_z(k) = G_z(k);
    end
end

% 最終的G(z)值(包括兩個純延時單元)
% G_z(Num_lattice+1) = G_oddk(m);
% G_z(3*Num_lattice+2) = G_oddMk(m);
% --- modified 20071015 ------------- %
G_z((M_ch-1)/2+1) = G_oddk(m);
G_z(M_ch+(M_ch-1)/2+1) = G_oddMk(m);

% 第3步：由2M個G(z)值得到H(z), 已經表示成H(e^jw)的形式
w0 = (0:0.01:2*pi);
Hz_MF = zeros(size(w0));

Hk_z = sym(zeros(M_ch,1));
for P_ch=1:M_ch
    for q=1:2*M_ch
        syms w;
        Hk_z(P_ch) = Hk_z(P_ch) + 2*cos((2*(P_ch-1)+1)*pi*((q-1)-(2*m*M_ch-1)/2)/(2*M_ch)+(-1)^(P_ch-1)*pi/4) * exp(-(q-1)*j*w)*G_z(q)*(1/sqrt(4*M_ch));
    end
    
    % plot the magnitude-frequency response
    for L=1:length(w0)
        w = w0(L);
        Hz_MF(L) = abs(eval(Hk_z(P_ch)));
    end
    Hz_MF = Hz_MF/max(Hz_MF);
    Hz_MF = 20*log10(Hz_MF);
    plot(w0(1:floor(length(w0)/2))/2/pi,Hz_MF(1:floor(length(w0)/2)));hold on;
end
% Hk_z = vpa(Hk_z);
hold off;axis([0,0.5,-80,5]);title('Analysis section');