pitch/period tracking using autocorrelation

%% Using Autocorrelation to track the local period of a signal
% This code is used as part of a youtube video demonstration 
% See http://youtube.com/ddorran
%
% Code available at https://dadorran.wordpress.com       
%
% The following wav file can be downloaded from 
%       https://www.dropbox.com/s/3y25abf1xuqpizj/speech_demo.wav
%% speech analysis example

[ip fs] = wavread('speech_demo.wav');
max_expected_period = round(1/50*fs);
min_expected_period = round(1/200*fs);
frame_len = 2*max_expected_period;

for k = 1 : length(ip)/frame_len -1;
    range = (k-1)*frame_len + 1:k*frame_len;
    frame = ip(range);
    
    %show the input in blue and the selected frame in red
    plot(ip);
    set(gca, 'xtick',[],'position',[ 0.05  0.82   0.91  0.13])
    hold on;
    temp_sig = ones(size(ip))*NaN;
    temp_sig(range) = frame;
    plot(temp_sig,'r');
    hold off
    
    %use xcorr to determine the local period of the frame
    [rxx lag] = xcorr(frame, frame);
    subplot(3,1,3)
    plot(lag, rxx,'r')
    rxx(find(rxx < 0)) = 0; %set any negative correlation values to zero
    center_peak_width = find(rxx(frame_len:end) == 0 ,1); %find first zero after center
    %center of rxx is located at length(frame)+1
    rxx(frame_len-center_peak_width : frame_len+center_peak_width  ) = min(rxx);
%     hold on
%     plot(lag, rxx,'g');
%     hold off
    [max_val loc] = max(rxx);
    period = abs(loc - length(frame)+1); 
    
    title(['Period estimate = ' num2str(period) 'samples (' num2str(fs/period) 'Hz)']);
    set(gca, 'position', [ 0.05  0.07    0.91  0.25])
    
    [max_val max_loc] = max(frame);
    num_cycles_in_frame = ceil(frame_len/period);
    test_start_positions = max_loc-(period*[-num_cycles_in_frame:num_cycles_in_frame]);
    index = find(test_start_positions > 0,1, 'last');
    start_position = test_start_positions(index);
    colours = 'rg';
    
    subplot(3,1,2)
    plot(frame);
    
    set(gca, 'position',[ 0.05 0.47 0.91 0.33])
    pause
    for g = 1 : num_cycles_in_frame
        if(start_position+period*(g) <= frame_len && period > min_expected_period)
            cycle_seg = ones(1, frame_len)*NaN;
            cycle_seg(start_position+period*(g-1):start_position+period*(g))  =...
                            frame(start_position+period*(g-1):start_position+period*(g));
            hold on
            
            plot(cycle_seg,colours(mod(g, length(colours))+1)) %plot one of the available colors
            hold off
        end
    end
    pause
end

%% synthesise a periodic signal to use as a basic demo
fs = 500;
T = 1/fs;
N = 250; % desired length of signal
t = [0:N-1]*T; %time vector 
f1 = 8; f2=f1*2; 
x = sin(2*pi*f1*t-pi/2) + sin(2*pi*f2*t);
plot(t, x)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Synthesised Signal');

%% Determine the autocorrelation function
[rxx lags] = xcorr(x,x);
figure
plot(lags, rxx)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')

%% Illustrate the auto correlation process
%function available from https://dadorran.wordpress.com
illustrate_xcorr(x,x) 

%% Identify most prominent peaks
% Most prominent peak will be at the center of the correlation function
first_peak_loc = length(x) + 1;

% Lots of possible ways to identify second prominent peak. Am going to use a crude approach
% relying on some assumed prior knowledge of the signal. Am going to assume
% that the signal has a minimum possible period of .06 seconds = 30 samples;
min_period_in_samples = 30; 
half_min = min_period_in_samples/2 ;

seq = rxx;
seq(first_peak_loc-half_min: first_peak_loc+half_min) = min(seq);
plot(rxx,'rx');
hold on
plot(seq)

[max_val second_peak_loc] = max(seq);
period_in_samples =  abs(second_peak_loc -first_peak_loc)
period = period_in_samples*T
fundamental_frequency = 1/period

%% Autocorrelation of a noisy signal 
x2 = x + randn(1, length(x))*0.2;
plot(x2)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Noisy Synthesised Signal');

[rxx2 lags] = xcorr(x2,x2);
figure
plot(lags, rxx2)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')

%% Autocorrelation technique can be problematic!
% Consider the following signal
f1 = 8; f2=f1*2; 
x3 = sin(2*pi*f1*t) + 5*sin(2*pi*f2*t);
plot(t, x3)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Synthesised Signal');

[rxx3 lags] = xcorr(x3,x3,'unbiased');
figure
plot(lags, rxx3)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')

seq = rxx3;
seq(first_peak_loc-half_min: first_peak_loc+half_min) = min(seq);
plot(seq)

[max_val second_peak_loc] = max(seq);
period_in_samples =  abs(second_peak_loc -first_peak_loc)

Advertisement

illustrate_xcorr – code for cross correlation demos

% This function illustrates the cross correlation process in action
%
% Usage:
%           fs = 1000;
%             T = 1/fs;
%             N = 500; % desired length of signal
%             t = [0:N-1]*T; %time vector 
%             f1 = 8; f2=f1*2; 
%             x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
%
%           % To step though each sample use the following:
%           illustrate_xcorr(x,x)
%           
%           % to step through using 50 steps use:
%           illustrate_xcorr(x,x, 50)
%
function illustrate_xcorr(x, y, varargin)
if(length(x) > length(y))
    y(end+1:length(x)) = 0; %zero pad so the signals are the same length
else
    x(end+1:length(y)) = 0; %zero pad so the signals are the same length
end

    num_steps = 2*length(x)-1;
if(nargin ==3)
    arg = varargin{1};
    if(isnumeric(arg))
        num_steps = ceil(abs(arg));
    end
end
if(nargin > 3)
    error('See help on this function to see how to use it properly')
end

[rxy lags] = xcorr(x,y); %cross correlate signals

disp('The signal being autocorrelated is shown in blue (two instances)')
disp('As you hit the space bar the lower plot will move into different lag positions')
disp('The correlation function shown in red is updated for each lag position');
disp('keep pressing the space bar to step through the illustration ...');

figure
plot_width = 0.3; plot_height = 0.25;

top_ax_h = subplot(3,1,1);
plot(x)
axis tight
set(top_ax_h, 'visible','off', 'units', 'normalized')
set(top_ax_h,'position', [0.5-plot_width/2 5/6-plot_height/2 plot_width plot_height])

mid_ax_h = subplot(3,1,2);
plot(y)
axis tight
set(mid_ax_h, 'visible','off', 'units', 'normalized')
set(mid_ax_h,'position', [0.5-plot_width/2-plot_width 5/6-3*plot_height/2-0.01 plot_width plot_height])

bottom_ax_h = subplot(3,1,3);
corr_h = plot(lags,rxy,'r');
axis tight
set(bottom_ax_h,'units', 'normalized','Ytick',[])
set(bottom_ax_h,'position', [0.5-plot_width*3/2 0.2-plot_height/2 plot_width*3 plot_height])
set(corr_h, 'Ydata', ones(1, length(rxy))*NaN); %clear the correlation funciton plot once its set up

normalised_shift_size = 2*plot_width/(num_steps-1);
corr_seg_len = length(rxy)/num_steps;
for k = 1 : num_steps
    if(k > 1)
        new_pos = get(mid_ax_h,'position') + [normalised_shift_size 0 0 0];
        set(mid_ax_h,'position', new_pos);
    end
    set(corr_h, 'Ydata', [rxy(1:round(corr_seg_len*k)) ones(1,length(rxy)-round(corr_seg_len*k))*NaN])
    pause
end

requantise

% this function takes a signal ip and modifies it so that
% occupies 2^(num_bits) quantisation levels
%
% ns = rand(1, 1000);
% op = requantise(ns, 2);
% plot(ns)
% hold on 
% plot(op,'r') % youshould be able to clearly see the 4 possible levels the
% new signal occupies
function op = requantise(ip, num_bits)
    num_levels = 2^num_bits;
    quantization_diff = (max(ip)-min(ip))/num_levels;
    quantization_levels = min(ip)+quantization_diff/2:quantization_diff:max(ip)-quantization_diff/2;
    op = zeros(1,length(ip));
    for k = 1: length(ip)
        [min_diff closest_level_index] = min(abs(quantization_levels-  ip(k)));
        op(k) = quantization_levels(closest_level_index);
    end