5623b11Recorder simulator: MATLAB jet-drive model and interactive web appJeremy Magland 1function basis = recorder_basis(P, note)
2% RECORDER_BASIS Modal fit of the passive resonator for one fingering.
3%
4% The series-loop admittance Y(omega) = 1 / (Z_window + Z_bore) from
5% recorder_bore.m relates the acoustic flow U through the mouth to the
6% pressure the jet imposes across it; its peaks are the bore resonances the
7% jet can lock to. Each peak is fitted with one second-order section,
8% Y_n(omega) = i omega B_n / (omega_n^2 - omega^2 + i omega g_n),
9% so that in the time domain U_n'' + g_n U_n' + omega_n^2 U_n = B_n dp/dt.
10% omega_n comes from parabolic interpolation of log|Y| on the grid, g_n from
11% the half-power bandwidth, and B_n = g_n Re(Y(omega_n)) so the fitted peak
12% height matches. The fit discards nothing a recorder player can reach:
13% every local maximum below max_mode_hz is kept, up to max_modes, however
14% weak (a cross-fingered lattice makes the upper resonances weak but they
15% still color the attack).
16%
17% For radiation the basis stores, per mode, the complex flow out of each
18% radiating aperture (the window, each open hole, the foot) per unit loop
19% flow, evaluated at the mode's centre frequency, plus the pressure
20% standing-wave profile along the bore for the mode-shape plot.
22geo = recorder_geometry(P);
23note = min(max(1, round(note)), geo.nNotes);
24fing = geo.fingerings(note, :);
26fgrid = 60:P.grid_hz:P.max_mode_hz;
27out = recorder_bore(P, geo, fing, fgrid);
28mag = abs(out.Y);
30% ---- peak picking with parabolic refinement ----
31fpk = [];
32qpk = [];
33ypk = [];
34nf = length(fgrid);
35for i = 2:nf - 1
36 if ~(mag(i) > mag(i - 1) && mag(i) >= mag(i + 1))
37 continue;
38 end
39 y1 = log(mag(i - 1)); y2 = log(mag(i)); y3 = log(mag(i + 1));
40 den = y1 - 2 * y2 + y3;
41 sh = 0;
42 if abs(den) > 1e-12
43 sh = min(0.5, max(-0.5, 0.5 * (y1 - y3) / den));
44 end
45 fc = fgrid(i) + sh * P.grid_hz;
46 yc = exp(y2 - 0.25 * (y1 - y3) * sh);
47 % local bandwidth from the two neighbouring grid points: for one section,
48 % |Y(f_n + df)|^2 / |Y(f_n)|^2 = g^2 / ((4 pi df)^2 + g^2), so each
49 % neighbour gives g = 4 pi |df| sqrt(r / (1 - r)). Staying local keeps the
50 % estimate off the skirts of adjacent resonances, which a half-power scan
51 % walks onto.
52 gs = [];
53 for s = [-1, 1]
54 df = fgrid(i + s) - fc;
55 r = (mag(i + s) / yc)^2;
56 if r < 0.999
57 gs = [gs, 4 * pi * abs(df) * sqrt(r / (1 - r))];
58 end
59 end
60 if length(gs) == 0
61 g = 2 * pi * P.grid_hz;
62 else
63 g = sum(gs) / length(gs);
64 end
65 fpk = [fpk, fc];
66 qpk = [qpk, 2 * pi * fc / g];
67 ypk = [ypk, yc];
68end
70n = min(length(fpk), round(P.max_modes));
71if n == 0
72 error('No bore resonances below max_mode_hz -- check the geometry.');
73end
74fpk = fpk(1:n);
75qpk = qpk(1:n);
76ypk = ypk(1:n);
78% ---- aperture weights and standing-wave profiles at the mode centres ----
79mo = recorder_bore(P, geo, fing, fpk);
81basis.n = n;
82basis.note = note;
83basis.f = reshape(fpk, n, 1);
84basis.omega = 2 * pi * basis.f;
85basis.Q = reshape(qpk, n, 1);
86basis.g = basis.omega ./ basis.Q;
87basis.B = basis.g .* reshape(ypk, n, 1);
88basis.apX = mo.apX; % nap x 1 axial positions
89basis.apW = mo.W; % nap x n complex flow weights
90basis.profX = mo.profX; % 1 x nsta
91basis.profP = abs(mo.profP); % nsta x n standing-wave magnitudes
92basis.fTarget = geo.fTarget(note);
93basis.Sm = geo.Sm;
94basis.gridF = fgrid;
95basis.gridY = mag;
96end