d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1function [fr, db] = wv_spectrum(x, fs, fmax, npts)
2% WV_SPECTRUM Magnitude spectrum of x in decibels, for plotting.
3%
4% Tapered with wv_window and zero-padded, then reduced to npts points up to
5% fmax by keeping the largest value in each bin of the plot, so that narrow
6% peaks survive the thinning. 0 dB is the largest peak.
7x = x(:);
8n = numel(x);
9w = wv_window(n);
10nfft = 2 ^ nextpow2(2 * n);
11X = abs(fft(w .* x, nfft));
12X = X(1:nfft/2+1);
13X = X / max(max(X), 1e-300);
14f = (0:nfft/2)' * fs / nfft;
15fr = linspace(0, fmax, npts);
16db = -200 * ones(1, npts);
17bw = fr(2) - fr(1);
18idx = floor(f / bw + 0.5) + 1;
19% real() because log10 is typed as possibly complex, which would keep this
20% loop off numbl's compiled path
21v = real(20 * log10(max(X, 1e-300)));
22for k = 1:numel(f)
23 j = idx(k);
24 if j >= 1 && j <= npts && v(k) > db(j)
25 db(j) = v(k);
26 end
27end
28db(db < -160) = -160;
29end