What am I messing up about this Special Allpass Filter?

Question

Blair 2023-9-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2017316-what-am-i-messing-up-about-this-special-allpass-filter

编辑： Balaji 2024-1-2

I found this allpass matlab code based on this whitepaper's equations:

It calculates the coefficients of a set of cascading allpasses in the form of biquad filters.

function sos = adf(f0,B,df,beta)
% adf.m
%
% Allpass dispersion filter design  
% 
% Ref.
% J.S. Abel, V. Valimaki, and J.O. Smith, "Robust, Efficient Design of
% Allpass Filters for Dispersive String Sound Synthesis," IEEE Signal
% Processing Letters, 2010.
%
% Parameters:
% f0 = fundamental frequency, Hz
% B = inharmonicity coefficent, ratio (e.g., B = 0.0001)
% df = design bandwidth, Hz 
% beta = smoothing factor (e.g., beta = 0.85)
%
% Examples:
% sos = adf(32.702,0.0002,2100,0.85);
% sos = adf(65.406,0.0001,2500,0.85);
% sos = adf(130.81,0.00015,4300,0.85);
%
% The output array sos contains the allpass filter coefficients
% as second-order sections.
% 
% Created: Vesa Valimaki, Dec. 11, 2009
% (based on Jonathan Abel's Matlab code)
%% initialization
% system variables
fs = 44100;		%% sampling rate, Hz
nbins = 2048;	%% number of frequency points
%% design dispersion filter
% period, samples; df delay, samples; integrated delay, radians
tau0 = fs/f0;   %% division needed
pd0 = tau0/sqrt(1 + B*(df/f0).^2);  %% division and sqrt needed
mu0 = pd0/(1 + B*(df/f0)^2);
phi0 = 2*pi*(df/fs)*pd0 - mu0*2*pi*df/fs;
% allpass order, biquads; desired phases, radians
nap = floor(phi0/(2*pi));
phik = pi*[1:2:2*nap-1];
% Newton single iteration
etan = [0:nap-1]/(1.2*nap) * df;  %% division needed
pdn = tau0./sqrt(1 + B*(etan/f0).^2); %% division and sqrt needed
taun = pdn./(1 + B*(etan/f0).^2);
phin = 2*pi*(etan/fs).*pdn;
theta = fs/(2*pi) * (phik - phin + (2*pi/fs)*etan.*taun) ./ (taun - mu0);
% division needed
% compute allpass pole locations
delta = diff([-theta(1) theta])/2;
cc = cos(theta * 2*pi/fs);
eta = (1 - beta.*cos(delta * 2*pi/fs))./(1 - beta);  %% division needed
alpha = sqrt(eta.^2 - 1) - eta;  % sqrt needed
% form second-order allpass sections
temp = [ones(nap,1) 2*alpha'.*cc' alpha'.^2];
sos = [fliplr(temp) temp];

I am unfamiliar with Matlab and am trying to translate it into a realtime language called genexpr from MaxMSP. I think that it isn't too hard to read the translation but something about it isn't right and I was wondering if someone with more knowledge in matlab would be able to see what might be wrong?

Thanks! :).

//biquad filter
gm_tpdf2 (x, b0, b1, b2, a1, a2) { 
	History z1, z2;
	
	y = fixdenorm(b0 * x + z1); //fixdenorm is used to Replace denormal values with zero.
	z1 = b1 * x - a1 * y + z2;
	z2 = b2 * x - a2 * y;
	
	return y;
}
f0 = in2;
B = in3;
df = in4;
beta = in5;
// Initialization
//f0 = 440.0;  // Fundamental frequency in Hz
//B = 0.0001;  // Inharmonicity coefficient
//df = 2500.0; // Design bandwidth in Hz
//beta = 0.85; // Smoothing factor
// Get the current sample rate
fs = samplerate;
// Calculate variables used in pole calculations
tau0 = fs / f0;
pd0 = tau0 / sqrt(1 + B * pow(df / f0, 2));
mu0 = pd0 / (1 + B * pow(df / f0, 2));
phi0 = 2 * pi * (df / fs) * pd0 - mu0 * 2 * pi * df / fs;
// Assuming nap (allpass order) is 4 based on your specification
//nap = 4;
nap = floor(phi0/(2*pi));
// Initial phase values for four biquads, assuming nap = 4
phik1 = pi * ((1*nap) - 1);
phik2 = pi * ((3*nap) - 1);//3 * pi; <was this
phik3 = pi * ((5*nap) - 1);//5 * pi;
phik4 = pi * ((7*nap) - 1);//7 * pi;
// Newton's iteration for calculating theta (pole angles), simplified for one step
//etan1 = 0 / (1.2 * nap) * df;
//etan2 = 1 / (1.2 * nap) * df;
//etan3 = 2 / (1.2 * nap) * df;
//etan4 = 3 / (1.2 * nap) * df;
//etan1 was zero - recheck calculations from paper
etan1 = 1 / (1.2 * nap) * df;
etan2 = 2 / (1.2 * nap) * df;
etan3 = 3 / (1.2 * nap) * df;
etan4 = 4 / (1.2 * nap) * df;
// Calculations for pdn, taun, phin, and theta for each biquad
pdn1 = tau0 / sqrt(1 + B * pow(etan1 / f0, 2));
theta1 = (fs / (2 * pi)) * (phik1 - 2 * pi * (etan1 / fs) * pdn1 + (2 * pi / fs) * etan1 * pdn1 / (pdn1 - mu0));
pdn2 = tau0 / sqrt(1 + B * pow(etan2 / f0, 2));
theta2 = (fs / (2 * pi)) * (phik2 - 2 * pi * (etan2 / fs) * pdn2 + (2 * pi / fs) * etan2 * pdn2 / (pdn2 - mu0));
pdn3 = tau0 / sqrt(1 + B * pow(etan3 / f0, 2));
theta3 = (fs / (2 * pi)) * (phik3 - 2 * pi * (etan3 / fs) * pdn3 + (2 * pi / fs) * etan3 * pdn3 / (pdn3 - mu0));
pdn4 = tau0 / sqrt(1 + B * pow(etan4 / f0, 2));
theta4 = (fs / (2 * pi)) * (phik4 - 2 * pi * (etan4 / fs) * pdn4 + (2 * pi / fs) * etan4 * pdn4 / (pdn4 - mu0));
// Compute allpass pole locations (alpha and eta) for each biquad
delta1 = theta1 / 2;
cc1 = cos(theta1 * 2 * pi / fs);
eta1 = (1 - beta * cos(delta1 * 2 * pi / fs)) / (1 - beta);
alpha1 = sqrt(pow(eta1, 2) - 1) - eta1;
delta2 = theta2 / 2;
cc2 = cos(theta2 * 2 * pi / fs);
eta2 = (1 - beta * cos(delta2 * 2 * pi / fs)) / (1 - beta);
alpha2 = sqrt(pow(eta2, 2) - 1) - eta2;
delta3 = theta3 / 2;
cc3 = cos(theta3 * 2 * pi / fs);
eta3 = (1 - beta * cos(delta3 * 2 * pi / fs)) / (1 - beta);
alpha3 = sqrt(pow(eta3, 2) - 1) - eta3;
delta4 = theta4 / 2;
cc4 = cos(theta4 * 2 * pi / fs);
eta4 = (1 - beta * cos(delta4 * 2 * pi / fs)) / (1 - beta);
alpha4 = sqrt(pow(eta4, 2) - 1) - eta4;
// Calculate biquad coefficients based on pole locations
// For Biquad 1
b0_1 = alpha1 * alpha1;
b1_1 = -2 * alpha1 * cc1;
b2_1 = 1;
a1_1 = -2 * alpha1 * cc1;
a2_1 = alpha1 * alpha1;
// For Biquad 2
b0_2 = alpha2 * alpha2;
b1_2 = -2 * alpha2 * cc2;
b2_2 = 1;
a1_2 = -2 * alpha2 * cc2;
a2_2 = alpha2 * alpha2;
// For Biquad 3
b0_3 = alpha3 * alpha3;
b1_3 = -2 * alpha3 * cc3;
b2_3 = 1;
a1_3 = -2 * alpha3 * cc3;
a2_3 = alpha3 * alpha3;
// For Biquad 4
b0_4 = alpha4 * alpha4;
b1_4 = -2 * alpha4 * cc4;
b2_4 = 1;
a1_4 = -2 * alpha4 * cc4;
a2_4 = alpha4 * alpha4;
//Inputs
x = in1;
//Biquads Cascade
apf1 = gm_tpdf2(x, b0_1, b1_1, b2_1, a1_1, a2_1);
apf2 = gm_tpdf2(apf1, b0_2, b1_2, b2_2, a1_2, a2_2);
apf3 = gm_tpdf2(apf2, b0_3, b1_3, b2_3, a1_3, a2_3);
apf4 = gm_tpdf2(apf3, b0_4, b1_4, b2_4, a1_4, a2_4);
//Outputs
out1 = apf1; //audio out
out2 = b0_2;
out3 = b1_2;
out4 = b2_2;
out5 = a1_2; 
out6 = a2_2;