Next  |  Prev  |  Up  |  Top  |  REALSIMPLE Top

#### Notes regarding the matlab script in Fig.18:

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

Next  |  Prev  |  Up  |  Top  |  REALSIMPLE Top