- The
`del`parameter specifies the number of initial samples to skip over in the input sound (the round-trip delay of the measurement system). It can be a very critical parameter if the filter-design method is sensitive to frequency-response phase. Since`invfreqsmethod`converts the measured frequency response to minimum phase, it is not sensitive to`del`over a wide range. - The frequency-zoom interval
`f1z`to`f2z`should include the resonance peak. Since the measured response is noisier at very low and very high frequencies, frequency-zooming improves the resulting match. - The estimated Q values printed at the end are

and the estimated pole frequencies are

While these estimates could be improved by various means [8], they appear to be more than sufficiently accurate for the application at hand. - The plot overlays in Figures 14 through 16 are plots of the
returned frequency responses
`Hp`(measured) and`Hd`(the model) versus sampled radian frequency`w`.

function [Q,wp,Hp,Hd,w] = invfreqsmethod(h,f1,f2,fs); % % INVFREQSMETHOD - use invfreqs to estimate Q and resonance % frequency from a resonator impulse response. % USAGE: % [Q,wp,w,Hp,Hd] = invfreqsmethod(h,f1,f2,fs); % where % h = impulse response (power of 2 length preferred) % f1 = lowest frequency of interest (Hz) % f2 = highest frequency of interest (Hz) % fs = sampling rate (Hz) % Q = estimated resonator Quality factor % wp = estimated pole frequency (rad/sec) % % invfreqs and invfreqz are VERY sensitive to where time 0 % is defined. Normalize this by converting the impulse % response to its minimum-phase counterpart: h = minphaseir(h); H = fft(h); L0 = length(h); k1 = round(L0*f1/fs); k2 = round(L0*f2/fs); Hp = H(k1:k2); L0p = length(Hp) w = 2*pi*[k1:k2]*fs/L0; wt = 1 ./ w.^2; % Nominal weighting proportional to 1/freq [Bh,Ah] = invfreqs(Hp,w,2,2,wt); Hph = freqs(Bh,Ah,w); % Denominator to canonical form A(s) = s^2 + (wp/Q) s + wp^2: Ahn = Ah/Ah(1); wp = sqrt(Ahn(3)); Q = wp/Ahn(2); fp = wp/(2*pi); disp(sprintf('Pole frequency = %f Hz',fp)); disp(sprintf('Q = %f',Q)); kp = L0*fp/fs - k1 + 1; % bin number of pole (+1 for matlab) k1z = max(1,floor(kp*0.5)); k2z = min(length(w),ceil(kp*1.5)); ndx = [k1z:k2z]; % BiQuad fit using z = exp(s T) ~ 1 + sT approximation: frn = fp/fs; % Normalized pole frequency (cycles per sample) R = 1 - pi*frn/Q; % pole radius theta = 2*pi*frn; % pole angle A = [1, -2*R*cos(theta), R*R]; B = [1 -1]; % zeros guessed by inspection of amp response Hd = freqz(B,A,w/fs); gain = max(abs(Hp))/max(abs(Hd)) B = B*gain; Hd = Hd*gain; |

Download faust_strings.pdf

REALSIMPLE Project — work supported in part by the Wallenberg Global Learning Network .

Released

Center for Computer Research in Music and Acoustics (CCRMA), Stanford University