Skip to content

Commit a385036

Browse files
author
Fiete Winter
authored
Fix phase offset of Linkwitz-Riley Allpass Filters and add testing function (#185)
* fix inverted phase of LR Allpass for orders divisible by 4 * add testing script and fix allpass phase offset for z-domain * add LR test script to test_all * Update test_linkwitz_riley.m axis label corrected * fix some comments and trailing newline * workaround for Octave bug
1 parent d6ea498 commit a385036

File tree

3 files changed

+138
-2
lines changed

3 files changed

+138
-2
lines changed

‎SFS_general/linkwitz_riley.m

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -92,10 +92,10 @@
9292
p = p(:); % octave creates row vectors
9393
if strcmp(domain,'z')
9494
z = 1./conj(p);
95-
k = prod(p);
95+
k = prod(p).*(-1).^(n/2);
9696
else
9797
z = -p;
98-
k = 1;
98+
k = (-1).^(n/2);
9999
end
100100
end
101101

‎validation/test_all.m

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,6 +119,10 @@
119119
test_modal_weighting(1);
120120
disp('Hit Enter to continue');
121121
pause
122+
disp('Running "test_linkwitz_riley(1)"');
123+
test_linkwitz_riley(1);
124+
disp('Hit Enter to continue');
125+
pause
122126
disp('Running "test_sphbesselh_zeros(1)"');
123127
test_sphbesselh_zeros(1);
124128
else
@@ -139,6 +143,7 @@
139143
test_imp_25d(0); ...
140144
test_spectrum_signal_conversion(0);
141145
test_modal_weighting(0);
146+
test_linkwitz_riley(0);
142147
test_sphbesselh_zeros(0);
143148
%test_wfs_iir_prefilter(0); ... % needs DSP Tooblox
144149
])

‎validation/test_linkwitz_riley.m

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
function status = test_linkwitz_riley(modus)
2+
%TEST_LINKWITZ_RILEY tests the design of the Linkwitz-Riley Filters
3+
%
4+
% Usage: status = test_linkwitz_riley(modus)
5+
%
6+
% Input parameters:
7+
% modus - 0: numerical
8+
% 1: visual
9+
%
10+
% Output parameters:
11+
% status - true or false
12+
13+
%*****************************************************************************
14+
% The MIT License (MIT) *
15+
% *
16+
% Copyright (c) 2010-2019 SFS Toolbox Developers *
17+
% *
18+
% Permission is hereby granted, free of charge, to any person obtaining a *
19+
% copy of this software and associated documentation files (the "Software"), *
20+
% to deal in the Software without restriction, including without limitation *
21+
% the rights to use, copy, modify, merge, publish, distribute, sublicense, *
22+
% and/or sell copies of the Software, and to permit persons to whom the *
23+
% Software is furnished to do so, subject to the following conditions: *
24+
% *
25+
% The above copyright notice and this permission notice shall be included in *
26+
% all copies or substantial portions of the Software. *
27+
% *
28+
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
29+
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
30+
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
31+
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
32+
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING *
33+
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER *
34+
% DEALINGS IN THE SOFTWARE. *
35+
% *
36+
% The SFS Toolbox allows to simulate and investigate sound field synthesis *
37+
% methods like wave field synthesis or higher order ambisonics. *
38+
% *
39+
% http://sfstoolbox.org sfstoolbox@gmail.com *
40+
%*****************************************************************************
41+
42+
status = false;
43+
44+
%% ===== Checking of input parameters ===================================
45+
nargmin = 1;
46+
nargmax = 1;
47+
narginchk(nargmin,nargmax);
48+
49+
%% ===== Main ============================================================
50+
fs = 44100; % sampling frequency in Hz
51+
fc = 1000; % crossover frequency in Hz
52+
f = logspace(-3,log10(fs/2),10000); % frequeny-axis in Hz
53+
omega = 2*pi*f; % angular frequency
54+
55+
ftype = {'low' 'high' 'all'}; % filter types
56+
Hs = zeros(3,length(f)); % s-domain frequency spectra
57+
Hz = Hs; % z-domain frequency spectra
58+
59+
for Nlr = [2 4 12] % filter orders
60+
61+
for tdx = 1:3
62+
% Laplace-Domain
63+
[zs,ps,ks] = linkwitz_riley(Nlr,2*pi*fc,ftype{tdx},'s');
64+
[soss,gs] = zp2sos(zs,ps,ks,'down','none');
65+
66+
s = [-omega.^2; 1i*omega];
67+
s(3,:) = 1;
68+
69+
if isoctave && strcmp(ftype{tdx}, 'low')
70+
% WORKAROUND for Octave bug (https://savannah.gnu.org/bugs/?51936)
71+
Hs(tdx,:) = gs.*prod( (soss(:,3:-1:1)*s)./(soss(:,4:6)*s),1);
72+
else
73+
Hs(tdx,:) = gs.*prod( (soss(:,1:3)*s)./(soss(:,4:6)*s),1);
74+
end
75+
76+
% z-Domain
77+
[zz,pz,kz] = linkwitz_riley(Nlr,fc./fs*2,ftype{tdx},'z');
78+
[sosz,gz] = zp2sos(zz,pz,kz,'down','none');
79+
80+
z = [exp(0.*omega./fs); exp(-1j.*omega./fs); exp(-2j.*omega./fs)];
81+
Hz(tdx,:) = gz.*prod( (sosz(:,1:3)*z)./(sosz(:,4:6)*z),1);
82+
end
83+
84+
% plotting
85+
if modus
86+
figure;
87+
subplot(2,2,1);
88+
plot_amp(f, Hs);
89+
title(sprintf('Amplitude Spectrum, s-Domain, order=%d', Nlr));
90+
subplot(2,2,2);
91+
plot_phase(f, Hs);
92+
title(sprintf('Phase Spectrum, s-Domain, order=%d', Nlr));
93+
subplot(2,2,3);
94+
plot_amp(f, Hz);
95+
title(sprintf('Amplitude Spectrum, z-Domain, order=%d', Nlr));
96+
subplot(2,2,4);
97+
plot_phase(f, Hz);
98+
title(sprintf('Phase Spectrum, z-Domain, order=%d', Nlr));
99+
end
100+
end
101+
102+
status = true;
103+
104+
end
105+
106+
function plot_amp(f, H)
107+
semilogx( ...
108+
f, db(H(1,:)), 'b', ...
109+
f, db(H(2,:)), 'r', ...
110+
f, db(H(1,:)+H(2,:)), 'g.', ...
111+
f, db(H(3,:)), 'k--' ...
112+
);
113+
legend('Lowpass', 'Highpass', 'Lowpass + Highpass', 'Allpass', 'Location', 'southwest');
114+
xlabel('f / Hz');
115+
ylabel('20 lg |H| / dB')
116+
xlim([f(1), f(end)]);
117+
ylim([-90, 10]);
118+
end
119+
120+
function plot_phase(f, H)
121+
semilogx( ...
122+
f, unwrap(angle(H(1,:))), 'b', ...
123+
f, unwrap(angle(H(2,:))), 'r', ...
124+
f, unwrap(angle(H(1,:)+H(2,:))), 'g.', ...
125+
f, unwrap(angle(H(3,:))), 'k--' ...
126+
);
127+
legend('Lowpass', 'Highpass', 'Lowpass + Highpass', 'Allpass', 'Location', 'southwest');
128+
xlabel('f / Hz');
129+
ylabel('angle(H) / rad')
130+
xlim([f(1), f(end)]);
131+
end

0 commit comments

Comments
 (0)