function [pitch_track,logprob] = wu_main(filename,dispflag,prct)
% function [pitch_track,logprob] = wu_main(filename,dispflag,prct)
%     Implement pitch tracking on the input wav file using the Wu algorithm
%
% INPUT:
%           filename - the name of the input audio file
%           dispflag - flag for plotting the pitch detection (default: 1)
%           prct - percentile (frame-wise) for the unvoiced audio 
%                  (default: 61%)
% OUTPUT:
%           path_org - pitch tracking result f0 (Hz) of Wu algorithm [1]
%           logprob_org - log probability of period (lags) of Wu algorithm [1]
% REFERENCE:
%           [1] Wu, Wang, and Brown, TSAP v11 n3 2003
% EXAMPLE:
%           [pitch_track,logprob] = wu_main('v3.wav');
% NOTE:
%           Reading the input audio file requires audioread package:
%           http://www.ee.columbia.edu/~dpwe/resources/matlab/audioread/
%           The path of audioread, arpath, should be updated properly.
%
% Created by Keansub Lee, kslee@ee.columbia.edu, keansub.lee@gmail.com
% [2011-01-30] Modified by Byung Suk Lee, bsl@ee.columbia.edu
% [2012-09-13] Last modified by Byung Suk Lee, bsl@ee.columbia.edu

if ~exist('dispflag','var'), dispflag = 1; end;
if ~exist('prct','var'), prct = 61; end;
arpath = '../audioread/';
if isempty(strfind(path,arpath)), addpath(arpath); end;

sr = 16000;
nchs = 48;
maxlags = 0.025 * sr; % maxlags = 400 when SR = 16000
% First, filter into subbands from 80 Hz to 800Hz , 55 channels
fmin = 100;bpo=16; q=8; n=2; ftype =2;

%% Read a wav file
[x]= audioread(filename,sr);

%% Calculate logprob
tic;
[logprob] = cal_ac_logprob(x,sr,fmin,bpo,q,n,ftype,nchs,maxlags);
toc;
 
%% HMM-based pitch-tracking
disp('HMM pitch-tracking for Wu algorithm')
tic;
[pitch_track] = pitch_tracking(normalise(logprob),prct, maxlags);
toc;

%% convert period (lag) to f0 (Hz)
pitch_period = pitch_track;
pitch_track = sr./pitch_track;
pitch_track(isinf(pitch_track)) = 0;

%% display
if dispflag == 1
dspnsec = length(x)/sr;
figure;
nhopsec = 10^floor(log10(dspnsec));

[B,F,T] = spectrogram(x,512,256,1024,sr);
subplot(211);
imagesc(T,F,20*log10(abs(B)));caxis([-80 35]);colorbar;
colormap(1-gray);
hold on
plot(0.01*(1:length(pitch_track)),pitch_track,'g.')
set(gca,'XTick',[0:nhopsec:dspnsec]);
set(gca,'XTickLabel',[0:nhopsec:dspnsec]);
ylim([0 1000])
axis xy;ylabel('freq / Hz');
title('Spectrogram of the input file');
legend('Wu output: f0')

subplot(212);imagesc(normalise(logprob));
axis xy;grid;caxis([-1.2039 5.84]);colorbar;
hold on
plot(pitch_period,'g.')
set(gca,'XTick',100*[0:nhopsec:dspnsec]);
set(gca,'XTickLabel',[0:nhopsec:dspnsec]);
ylim([0 maxlags])
ylabel('lags');
title('Summary autocorrelation (SAC) for lags by Wu algorithm');
legend('Wu output: Period')

end
