Matlab for Passive Reflectance Synthesis Method 1

(Python implementation with upgrades: https://github.com/suth1807/StringDispersion/blob/main/bridge_filter.py)

fs = 8192;  % Sampling rate in Hz (small for display clarity)
fc = 300;   % Upper frequency to look at
nfft = 8192;% FFT size (spectral grid density)

nspec = nfft/2+1;
nc = round(nfft*fc/fs);
f = ([0:nc-1]/nfft)*fs;

% Measured guitar body resonances
F = [4.64 96.52 189.33 219.95]; % frequencies
B = [ 10   10     10     10  ]; % bandwidths

nsec = length(F);

R = exp(-pi*B/fs);     % Pole radii
theta = 2*pi*F/fs;     % Pole angles
poles = R .* exp(j*theta);
A1 = -2*R.*cos(theta);  % 2nd-order section coeff
A2 = R.*R;              % 2nd-order section coeff
denoms = [ones(size(A1)); A1; A2]'
A = [1,zeros(1,2*nsec)];
for i=1:nsec,  
    % polynomial multiplication = FIR filtering:
    A = filter(denoms(i,:),1,A);
end;

Now A contains the (stable) denominator of the desired bridge admittance. We want now to construct a numerator which gives a positive-real result. We'll do this by first creating a passive reflectance and then computing the corresponding PR admittance.

g = 0.9;         % Uniform loss factor

B = g*fliplr(A); % Flip(A)/A = desired allpass

Badm = A - B;
Aadm = A + B;
Badm = Badm/Aadm(1); % Renormalize
Aadm = Aadm/Aadm(1);

% Plots

fr = freqz(Badm,Aadm,nfft,'whole');

nc = round(nfft*fc/fs);
spec = fr(1:nc);
f = ([0:nc-1]/nfft)*fs;
dbmag = db(spec);
phase = angle(spec)*180/pi; 

subplot(2,1,1);
plot(f,dbmag); grid;
title('Synthetic Guitar Bridge Admittance');
ylabel('Magnitude (dB)');
subplot(2,1,2);

plot(f,phase); grid;
ylabel('Phase (degrees)');
xlabel('Frequency (Hz)');