% Figure 12(c) of "Cancellation model of pitch perception" (JASA 97)

global sr;				% sampling rate
sr = 20800;
level = 90;			% dB
nchan = 2;			% number of BM channels
flo = 1040;			% Hz - lowest BM channel frequency (fhi = sr/2)
lag = 1/flo;

range = 0.0075; 	% s	- range for ACF
windowsize = 0.025;	% s	- window size for ACF


y = synth_ramp;

% scale
level = level - 30;	% 30 dB corresponds to RMS=1 in M&H's model.
y = rms(y, fromdB(level));

% synthesize & apply gammatone filter
[A, B] = MakeERBFilters(sr, nchan, flo);
yy = ERBFilterBank(A, B, y);
yy = flipud(yy);

y = yy(1,:);		% take lowest channel

y = hwrect(y);		% half-wave rectify
y = cfilt(y, lag*sr);
y = hwrect(y);		% half-wave rectify

yy = rdf2(y, range*sr, windowsize*sr);
yy = yy / max(yy);
plot1(yy);
dump('a', yy);

hold on;

y = synth_damp;

% scale
level = level - 30;	% 30 dB corresponds to RMS=1 in M&H's model.
y = rms(y, fromdB(level));

% synthesize & apply gammatone filter
[A, B] = MakeERBFilters(sr, nchan, flo);
yy = ERBFilterBank(A, B, y);
yy = flipud(yy);

y = yy(1,:);		% take lowest channel

y = hwrect(y);		% half-wave rectify
y = cfilt(y, lag*sr);
y = hwrect(y);		% half-wave rectify

yy = rdf2(y, range*sr, windowsize*sr);
yy = yy / max(yy);
plot1_bis(yy, 'r--');
dump('b', yy);
hold off



set(get(gca, 'title'), 'String',...
	'"Cancellation model of pitch perception" - Fig. 12c');