Dispersion Filter Design

This Wiki page contains supporting information for the following paper:

Reference

"Robust, Efficient Design of Allpass Filtersfor Dispersive String Sound Synthesis", by Jonathan S. Abel, Vesa Välimäki, and Julius O. Smith III, IEEE SIGNAL PROCESSING LETTERS, VOL. 17, NO. 4, APRIL 2010

Abstract

An efficient allpass filter design method is introduced to match the dispersion characteristics of vibrating stiff strings. The proposed method designs an allpass filter in cascaded biquad form directly from the target group delay, placing the poles at frequencies at which the group delay area function achieves odd integer multiples of , and fixing the pole radii according to a smoothness parameter. The pole frequencies are seen to be roots of quartic polynomials, and an efficient approximation to the desired roots is provided. Design examples show the method to outperform a previous closed-form design. Furthermore, the proposed method can achieve an arbitrarily wide bandwidth of good approximation by increasing the filter order, as the method is numerically robust and yields stable allpass filters.

Matlab for dispersion allpass filter design

File link

http://www.acoustics.hut.fi/publications/papers/spl-adf/adf.m

Text-Format Listing

function sos = adf(f0,B,df,beta)
% adf.m
%
% Allpass dispersion filter design  
% 
% REFERENCE:
% "Robust, Efficient Design of Allpass Filters for 
%  Dispersive String Sound Synthesis," by
%  Jonathan S. Abel, Vesa Valimaki, and Julius O. Smith, 
%  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)
% (ported to CCRMA Wiki by jos 2010-02-23)

%% 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];