Matlab listing: oboeanal.m

The following Matlab script illustrates use of the findpeaks function above to determine the pitch of an oboe tone (given the general location of the correct spectral peak) and configure a spectrum analysis using the rectangular, Hamming, and Blackman windows. This script was used to create Figures 3.16 through 3.18.

diary './dia/oboeanal.dia';
diary on;

zpf = 5; % desired min zero-padding factor in all FFTs
fn = '../wav/oboe.ff.C4B4.wav';
[x,fs,nbits] = wavread(fn);
nx = length(x);
dur = nx/fs;
disp(sprintf(['Read %s, fs = %f, nbits = %d, ',...
              'length = %d samples = %0.1f sec'],...
             fn,fs,nbits,nx,dur));
soundsc(x,fs);
M = fs/10; % 100 ms window for pitch estimation
N = 2^nextpow2(M*zpf) % zero pad
Xzp = fft(x,N);

figure(K-1);
fmax = 500; % maximum freq in Hz
kmax = fmax*N/fs; % maximum freq in bins
prange = 1:kmax+1;
Xdb = dbn(Xzp(prange)); % dB normalized to 0 max
f = fs*[0:N-1]/N/1000; % frequency axis in kHz
plot(f(prange),Xdb);
ylabel('Magnitude (dB)');
xlabel('Normalized Frequency (cycles/sample)');
[amps,bins] = findpeaks(Xdb,1);
k0 = bins(1)-1;
f0 = k0*fs/N;
nP = fs/f0; % samples per period (used below)
disp(sprintf('Estimated pitch = %f Hz',f0));

% Test pitch estimate:
pause(ceil(dur)); % let previous sound finish
T = 1/fs;
t = 0:T:dur;
y = sawtooth(2*pi*f0*t);
sound(y,fs)

winnames={'none','boxcar','none','hamming','none','blackman'};

% Configure window analysis

for K=[2 4 6]
  wintype = winnames{K};
  M = ceil(K*nP); % min Blackman length
  cmd = sprintf('w = %s(M);',wintype);
  eval(cmd);
  Nw = 2^nextpow2(M*zpf) % zero pad
  xw = x(1:M) .* w;
  Xwzp = fft(xw,Nw);
  %fmax = 5*f0; % maximum freq in Hz
  fmax = fs/4;
  kmax = fmax*Nw/fs; % maximum freq in bins
  prange = 1:kmax+1;
  Xwdb = dbn(Xwzp(prange)); % dB normalized to 0 max

  figure(K);
  subplot(2,1,1);
  plot(xw,'-k');
  axis tight;
  title(sprintf('Oboe data - %s window, K=%d',wintype,K));
  ylabel('Amplitude');
  xlabel('Time (samples)');
  subplot(2,1,2);
  f = fs*[0:Nw-1]/Nw/1000; % frequency axis in kHz
  plot(f(prange),Xwdb,'-k'); grid on;
  axis([f(prange(1)) f(prange(end)) -90 0]);
  ylabel('Magnitude (dB)');
  xlabel('Frequency (kHz)');

  plotfile = sprintf('../eps/oboe%s.eps',wintype);
  saveplot(plotfile);

end

diary off;