function [a1,Dx] = get_tdf_coeff(f0,B,M)
% GET_TDF_COEFF Compute the tunable dispersion filter coefficient.
%
% Parameters:
%     f0: fundamental frequency (Hz)
%     B: Inharmonicity coefficient
%     M: Number of filters in cascade
%     a1: dispersion filter coefficient value
%     Dx: estimated delay at the dc
%
% Example:
%     a1 = get_tdf_coeff(27.5,0.0001,8);
%     b = [a1 1];
%     a = [1 a1];
%
% -------------------------------------------------------------------
% Last modified 6.6.2006 by Jukka Rauhala
% Copyright (c) 2006 by Jukka Rauhala
% Helsinki University of Technology
% Laboratory of Acoustics and Audio signal processing
% www.acoustics.hut.fi
%
% For technical details, see:
%   [1] J. Rauhala and V. Välimäki, "Dispersion Modeling in Waveguide Piano
%       Synthesis Using Tunable Allpass Filters," accepted for publication in 
%       the 9th International Conference on Digital Audio Effects (DAFx-06), 
%       Montreal, Canada, September 2006.
% -------------------------------------------------------------------

% Compute the key index corresponding to f0
ratio = 1.05946309436;
key0 = log(f0*ratio/27.5)/log(ratio);
keyi = round(key0);

% Parameterization coefficients
p1 = [0.012644 0.060623];
p2 = [-0.0082511 1.9737];
k_coeffs = [-0.0017948 -0.023290 -2.9274];
c_coeffs = [polyval(p1,log(M)) polyval(p2,log(M))];

% Calculate parameter D (the phase delay at the dc)
C = polyval((c_coeffs),log(B));
k = polyval((k_coeffs),log(B));
Dexp = -key0*exp(k)+exp(C);
D = (exp(Dexp));
Dx = D*M;

% Determine first-order Thiran allpass filter coefficients
Q = D-1;
a1 = -1*(Q/(Q+2));