Difference between revisions of "Dispersion Filter Design"

From CCRMA Wiki
Jump to: navigation, search
(Created page with 'Matlab for dispersion filter design. (To be completed.)')
 
Line 1: Line 1:
Matlab for dispersion filter design. (To be completed.)
+
[[Category:Projects]]
 +
 
 +
==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
 +
 
 +
==Matlab for dispersion allpass filter design==
 +
 
 +
  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];

Revision as of 10:49, 23 February 2010


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

Matlab for dispersion allpass filter design

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