DSP Foundations Notes

January 27, 2015 2 comments

I’ve put together a set of notes introducing Digital Signal Processing. There’s only 25 pages but I hope that its enough to get someone started.

Download DSP Foundations

interpLTspice

November 8, 2014 2 comments
% This function interpolates data returned from the runLTspice function so
% as to make the time vector evenly spaced.
%
% example usage:
%   asc_file = 'filename.asc';
%   volt_curr_list =  runLTspice(asc_file);
%   results = runLTspice(asc_file, volt_curr_list);
%   fs = 10000;
%   interp_results = interpLTspice(results, fs); 
%
% By David Dorran (david.dorran@dit.ie)
function data_struct = interpLTspice(data_struct, fs)
%check if the data is structured as expected
data_ok = 0;
if(isstruct(data_struct))
    data_ok = 1;
    if ~(isfield(data_struct,'data') && isfield(data_struct,'data_name') && isfield(data_struct,'time_vec'))
        data_ok = 0;
    end
end
if(~data_ok)
    error('The variable passed to this function is not a data structure returned from runLTspice and cannot be processed')
end

if~(isnumeric(fs) && length(fs) ==  1 && fs > 0)
    error('The second parameter used in this function should be a positive number representing the sampling frequency')
end



try t = data_struct.time_vec(1):1/fs:data_struct.time_vec(end);
catch
    error(sprintf('A sampling rate of %f resulted in too much data (memory limits exceeded).\n\nTry a different sampling frequency using interpLTspcice(data_struct,fs);', fs));
end
for k = 1:length(data_struct.data)
    data_struct.data{k} =  interp1(data_struct.time_vec, data_struct.data{k}, t);
end
data_struct.time_vec = t;

runLTspice demo

October 29, 2014 2 comments
%% Demonstration of the runLTspice function
%
% Demo video at http://youtu.be/ax9H4eKuZv4
%
% Download runLTspice from:
%         https://dadorran.wordpress.com/2014/10/29/runltspice-m/
% 
% Download this code from https://dadorran.wordpress.com
%
% LTspice schematic available from:
%       https://www.dropbox.com/s/y7j1vmyscvppsn4/RC.asc?dl=0
%

%% read the voltage and currents from the circuit
result = runLTspice('C:\Temp\RC.asc', 'V(n002)');

plot(result.time_vec, result.data{1})
xlabel('Time(seconds)')
ylabel('Voltage')
title(result.data_name{1})

%% apply a step input and plot the capacitor current and voltage wavforms
fs = 200;
step = [zeros(1, fs/2) ones(1, fs/2)];

res = runLTspice('C:\Temp\RC.asc', 'V(n002), I(C1)' , 'V1', step, fs);

plot(res.time_vec, res.data{1})
hold on
plot(res.time_vec, res.data{2}*1000,'r')
xlabel('Time(seconds)')
ylabel('Amplitude')
hold off
legend([ res.data_name{1} ' (v)'],  [res.data_name{2} ' (mA)'])

%% apply a noisy square wave and compare the supply voltage with capacitor voltage
sq_wav = [step step step step];
high_f_noise = filter(1, [1 0.9], randn(1, length(sq_wav))*0.1);
ns_sq_wav = sq_wav + high_f_noise;

r = runLTspice('C:\Temp\RC.asc', 'V(n001), V(n002)' , 'V1', ns_sq_wav, fs);

plot(r.time_vec, r.data{1})
hold on
plot(r.time_vec, r.data{2},'r')
legend('supply voltage','smoothed voltage')
hold off

runLTspice.m

October 29, 2014 Leave a comment
% runLTspice is a function which can be used to read/obtain the voltage and/or
% current waveforms produced when an LTspice schematic (.asc file) is
% run/simulated.
%
% The function also allows you to replace the voltage and current sources
% which are used in the LTspice schematic (.asc file) with signals created in matlab.
%
% obtaining data usage:
%   volt_curr_list =  runLTspice('filename.asc'); %returns a list of all voltages and currents in the schematic.
%   results = runLTspice('filename.asc', volt_curr_list); %returns the data of the voltages and currents specified in parameter 2 (parameter 2 is a comma separated list of voltages and currents that the user wants returned)
%   plot(results.time_vec, results.data{1}); % plot the first set of data received
%   xlabel('Time (seconds)')
%   title(results.data_name{1})
%
% alter source data usage:
%  [volt_curr_list sources_list] = runLTspice('filename.asc');
%  fs = 100; %100 Hz sampling rate
%  step_input = [zeros(1, fs*1) ones(1, fs*10)]; % unit step input of duration 11 seconds 
%  first_source = sources_list(1: findstr(',',sources_list)-1);
%  if(~length(first_source))
%      first_source = sources_list;
%  end
%  results = runLTspice('filename.asc', volt_curr_list, first_source, step_input, fs);
%
% Written by: David Dorran (david.dorran@dit.ie)
%
% Acknowledgements:
% This function uses the LTspice2Matlab function written by Paul Wagner to capture data from a raw LTspice file.
% see - http://www.mathworks.com/matlabcentral/fileexchange/23394-fast-import-of-compressed-binary-raw-files-created-with-ltspice-circuit-simulator

function [varargout] = runLTspice(asc_file, varargin)
if(~exist('LTspice2Matlab'))
    error(sprintf('This function requires the function LTspice2Matlab written by Paul Wagner \n- Download from: http://www.mathworks.com/matlabcentral/fileexchange/23394-fast-import-of-compressed-binary-raw-files-created-with-ltspice-circuit-simulator'))
end
read_data = 0;
display_LTspice_variables = 0;
incorrect_output_parameters = 0;
incorrect_input_parameters = 0;
num_sources = 0;
if(nargin == 1)
    % just need to let the user know what voltages/currents can specified
    % in parameter 2 (parameter 2  is a comma separated list of
    % voltage/currents the the user requires) and what voltage current
    % sources can be modified with user specified data.
    if(nargout == 0)
        display_LTspice_variables = 1;
    elseif( nargout > 2)
        error(sprintf('You must specify either none, 1 or 2 return parameters when obtaining the voltage/current wavforms and sources from an LTspice asc file.\n\nE.g.\n\trunLTspice(''filename.asc'');\nor\n\tvolt_curr_list = runLTspice(''filename.asc'');\nor\n\t[volt_curr_list sources_list] = runLTspice(''filename.asc'');'))
    end
elseif(nargin == 2)
    % No need to modify the asc file with voltage/current pwl sources as
    % the user just wants to be able to capture data from the asc file
    read_data = 1;

    if(~isstr(varargin{1}))
        error(sprintf('The second parameter must be string containing a comma separated list of LTspice voltage/currents you would like returned. \n\nType ''runLTspice(''%s'')'' to get a list of valid values. \n\nType ''help runLTspice'' for usage examples. ', asc_file));
    end
    if(nargout  > 1)
        incorrect_output_parameters = 1;
    end
elseif(nargin > 2)
    read_data = 1;
    if(nargout  > 1)
        incorrect_output_parameters = 1;
    end
    mynargin = nargin - 2;
    if(mod(mynargin, 3))
        incorrect_input_parameters = 1;
    else
        num_sources = mynargin/3;
        for k = 1 : num_sources;
            if ~isstr(varargin{1+(k-1)*3+1}) || isstr(varargin{1+(k-1)*3+2}) || isstr(varargin{1+(k-1)*3+3}) || length(varargin{1+(k-1)*3+3}) > 1
                incorrect_input_parameters = 1;
            end
            sources2update{k} = varargin{1+(k-1)*3+1};
        end
    end
end

if(incorrect_output_parameters)
    error(sprintf('You can only specify one output parameter when retrieving data\n\nType ''help runLTspice'' for example usage.'))
end

if(incorrect_input_parameters)
    error(sprintf('When specifying data for either a supply voltage or current source three parameters are required for each source \ni.e. the source name (a string), the data (a vector/array of numbers), the sampling rate (a number).\n\ne.g. \n\tresults = runLTspice(''%s'', ''V(n002), V(n005) '', ''V1'', [0:0.1:10], 2, ''V2'', ones(1, 20), 4);\n\nNOTE: The values voltages V(n002), V(n005) and sources V1 and V2 must exist in ''%s''.\n\nTo determine the valid voltage/currents and sources in your LTspice asc file run the following command:\n\t runLTspice(''%s'')', asc_file,asc_file, asc_file))
end

if(~exist(asc_file))
    error(['Could not find the LTspice file ' asc_file ])
end

asc_path = fileparts(asc_file);
if(~length(asc_path))
    asc_path = pwd;
end
asc_path(end+1) = '\';

raw_file = regexprep(asc_file, '.asc$', '.raw');

this_file_path = which('runLTspice');
this_dir_path = this_file_path(1:end-length('runLTspice.m'));

config_file = [this_dir_path 'runLTspice_config.txt'];
if(exist(config_file))
    ltspice_path = fileread(config_file);
else
    ltspice_path = 'C:\Program Files\LTC\LTspiceIV\';
end

% See if the LTspice exe can be found in the path ltspice_path. If not get
% the user to enter the path. Once a valid path is entered then store it
% in the config file
exe_file = [ltspice_path 'scad3.exe'];
path_needs_updating = 0;
while ~exist(exe_file)
    disp(['The LTspice execuatble scad3.exe could not be found in ' ltspice_path '.']);
    if(~path_needs_updating)
        disp('Once the correct location is entered it will be recorded for future use of this function.')
    end
    disp(' ')
    ltspice_path = input('Please enter the folder where the executable exists => ','s');
    ltspice_path= strtrim(ltspice_path);
    if(ltspice_path(end) == '/')
        ltspice_path(end) ='\';
    end
    if(ltspice_path(end) ~= '\')
        ltspice_path(end+1) ='\';
    end

    exe_file = [ltspice_path 'scad3.exe'];
    path_needs_updating = 1;
end

if(path_needs_updating)
    disp('Updating LTspice exe location ...');
    try fp = fopen(config_file,'w');
        fprintf(fp, '%s',ltspice_path );
        fclose(fp);
        disp('Update Complete.')
    catch
        disp(['Could not open the configuration file ' config_file ' to record the exe location.']);
        disp('You will need to enter the path to the LTspice manually every time you use this function')
    end
end


% Parse the .asc file to idenitfy the current and voltage sources in the
% schematic. Also sote the line numbers where the values of these sources
% are located - these values will be updated with pwl files in the event
% that the user wishes to update the sources in the schematic
fid = fopen(asc_file);

k=1;
tline{k} = fgetl(fid);
voltage_found = 0;
current_found = 0;
voltage_names = {};
current_names = {};
symbol_num = 0;
voltage_symbol_nums = [];
current_symbol_nums = [];
while ischar(tline{k})
    k = k + 1;
    tline{k} = fgetl(fid);

    if(length(findstr('SYMBOL ',tline{k})))
        voltage_found = 0;
        current_found = 0;
        symbol_num = symbol_num + 1;
    end
    if(length(findstr('SYMATTR Value',tline{k})))
        symbol_value_line(symbol_num) = k;
    end
    if(length(findstr('SYMBOL voltage',tline{k})) || length(findstr('SYMBOL Misc\\signal ',tline{k}))|| length(findstr('SYMBOL Misc\\cell ',tline{k}))|| length(findstr('SYMBOL Misc\\battery ',tline{k})))
        voltage_found = 1;
    end
    if(length(findstr('SYMBOL current',tline{k})))
        current_found = 1;
    end
    if(voltage_found && length(findstr('SYMATTR InstName', tline{k})))
        voltage_names{end+1} = tline{k}(length('SYMATTR InstName')+2:end);
        voltage_symbol_nums(end+1) = symbol_num;
    end
    if(current_found && length(findstr('SYMATTR InstName', tline{k})))
        current_names{end+1} = tline{k}(length('SYMATTR InstName')+2:end);
        current_symbol_nums(end+1) = symbol_num;
    end
end
num_asc_file_lines = k;
source_names = [voltage_names current_names];
source_symbol_nums  = [voltage_symbol_nums  current_symbol_nums] ;
fclose(fid);

%% Just need to let the user know what sources can be used and the voltages/currents that can be read. No need to update the asc file
if(read_data == 0)

    sim_result = test_asc_simulation(asc_file,raw_file);
    if(sim_result == 0)
        error(sprintf('There was a problem running the LTspice file %s. \n\n Please open the file in LTspice and verify that it can be simulated', asc_file));
    end  
    raw_data = LTspice2Matlab(raw_file);
    voltage_current_list = '';
    for k = 1: raw_data.num_variables
        voltage_current_list = [voltage_current_list raw_data.variable_name_list{k} ', '];
    end
    voltage_current_list(end-1:end) = []; %remove trailing comma
    
    source_list = '';
    for k = 1: length(source_names)
        source_list = [source_list source_names{k} ', '];
    end
    source_list(end-1:end) = [];
    
    if(nargout)
        varargout{1} = voltage_current_list;
        if(nargout ==2)
            varargout{2} = source_list;
        end
    else
      disp(sprintf('List of voltage and currents that can be read:\n\t%s\n\nList of voltage and/or current sources that can be written to:\n\t%s', voltage_current_list, source_list))
    end
        
    return;
end

%identify the variables the user wants returned
vars2read = strtrim(strread(varargin{1}, '%s','delimiter',','));
if(~length(vars2read))
    error('The list of voltages and currents to read must be contained in a comma separated string.');
end
    
%% no need to change the pwl file so just run the specified asc file and return the results
if(num_sources == 0)
    sim_result = test_asc_simulation(asc_file,raw_file);
    if(sim_result == 0)
        error(sprintf('There was a problem running the LTspice file %s. \n\n Please open the file in LTspice and verify that it can be simulated', asc_file));
    end  
    raw_data = LTspice2Matlab(raw_file);
    
    varargout = grab_specified_data(vars2read, raw_data);
    
    return;
end

%% Need to update asc file and create pwl files for each source specified if you reach this stage

% calculate the simulation duration so as to match the duration of the longest
%(in time) source signal
sim_duration = 0;
for k = 1: num_sources
  duration = length(varargin{1+(k-1)*3+2})/varargin{1+(k-1)*3+3}; % = length(data)/fs  
  if(duration > sim_duration)
      sim_duration = duration;
  end
  source_found = 0;
    for m = 1: length(source_names)
        if strcmp(sources2update{k}, source_names{m})
            sources2update_value_line(k) = symbol_value_line(source_symbol_nums(m));
            source_found = 1;
        end
    end
    if(~source_found)
        error(sprintf('Unable to find the source %s in the schematic.\n\nType runLTspice(''%s''); at the command line to see a list of valid sources', sources2update{k}, asc_file));
    end
end

%ceate pwl files for each souce
for k = 1: num_sources
    % create pwl files to use to drive the voltage sources
    tmp_pwl_file{k} = [asc_path 'runLTspice_tmp_pwl_' num2str(k) '.txt'];
    create_pwl(varargin{1+(k-1)*3+2},varargin{1+(k-1)*3+3}, tmp_pwl_file{k});
end

% create a temporary file which will basically be a copy of the asc file passed
% to the function with voltage and current sources updated with pwl files.
% Also the simulation duration will be updated.

tmp_asc_file = [asc_path 'runLTspice_tmp_asc.asc'];
tmp_raw_file = [asc_path 'runLTspice_tmp_asc.raw'];


% update the temporary file with pwl files and the duration the simulation
% is to run.
f_op_id = fopen(tmp_asc_file,'w');
if(f_op_id <= 0)
    error(sprintf('This function needs to be able to write to the folder in which the LTspioce asc file exists. \nCopy the file %s to a folder you can write to and try again. ', asc_file));
end
k = 0;
while( k < num_asc_file_lines - 1)
    k = k + 1;
    res = findstr('!.tran ', tline{k});
    if(length(res))
        res2 = regexp(tline{k}, '.tran (?<info>.*?)(?<settings>(steady|startup|nodiscard|uic|$).*)', 'names');
        res_info = res2(1).info;
        tran_settings = res2(1).settings;
        num_spaces = length(findstr(' ', res_info) );
        if(num_spaces == 1 || num_spaces == 0)
            tran_info = [num2str(sim_duration) ' '];
        else
            tran_info = regexprep(res_info, ' .+? ',  [' ' num2str(sim_duration) ' ']  ,'once');
            
        end
        tline{k}(res(1)+length('!.tran '):end) = [];
        tline{k} = [tline{k} tran_info tran_settings];
    end
    
    index = find(sources2update_value_line==k);
    if(length(index))
            fwrite(f_op_id, ['SYMATTR Value PWL file="' tmp_pwl_file{index} '"' ]);
            fprintf(f_op_id, '\n');
    else
        fwrite(f_op_id, tline{k});
        fprintf(f_op_id, '\n');
    end
end

fclose(f_op_id);
sim_result = test_asc_simulation(tmp_asc_file,tmp_raw_file);
if(sim_result == 0)
    error(sprintf('There was a problem running the LTspice file %s. \n\n Please open the file in LTspice and verify that it can be simulated', asc_file));
end  
raw_data = LTspice2Matlab(tmp_raw_file);

varargout = grab_specified_data(vars2read, raw_data);
delete([asc_path 'runLTspice_tmp_*'])
return;

    function result = test_asc_simulation(test_filename,test_rawfile)
        % Will need to check if a raw file was created/updated to make sure
        % that the simulation ran successfully
        orig_timestamp =0;
        if(exist(test_rawfile))
            d = dir(test_rawfile);
            orig_timestamp = datenum(d.date);
        end
        % run/simulate the spice file
        system_command = ['"' exe_file '" -b -run "' test_filename '"'];
        system(system_command);
        result = 1;
        if(exist(test_rawfile))
            d = dir(test_rawfile);
            if( datenum(d.date) < orig_timestamp)
                result = 0; %unsuccessful simulation
            end
        else
            result = 0;; %unsuccessful simulation
        end
    end


    function create_pwl(data, fs, filename)

        file_id = fopen(filename,'w');
        if(file_id < 1)
            error(sprintf('This function needs to be able to write to the folder in which the LTspioce asc file exists. \nCopy the file %s to a folder you can write to and try again. ', asc_file));
        end
        for k = 1: length(data)
           fprintf(file_id, '%6.6f %6.6f ' , (k-1)/fs, data(k)); 
        end   
        fclose(file_id);
    end

    function data_details = grab_specified_data(vars, raw_info)
            for k = 1: length(vars)
                found(k) = 0;
                for m = 1:raw_info.num_variables
                    if(strcmp(vars{k},raw_info.variable_name_list{m}))
                        found(k) = m;
                    end
                end
                if(found(k))
                    data_details{1}.data_name{k} = vars{k};
                    data_details{1}.data{k} = raw_info.variable_mat(found(k),:);
                else

                    data_details{1}.data_name{k} = [vars{k} ' not found!'];
                    data_details{1}.data{k} = [];
                end
            end
            data_details{1}.time_vec = raw_info.time_vect;
    end

end

locate_peaks matlab function

October 6, 2014 Leave a comment
function indices = locate_peaks(ip)
%function to find peaks
%a peak is any sample which is greater than its two nearest neighbours
    index = 1;
    num = 2;
    indices = [];
    for k = 1 + num : length(ip) - num
        seg = ip(k-num:k+num);
        [val, max_index] = max(seg);
        if max_index == num + 1
            indices(index) = k;
            index = index + 1;
        end;
    end;
Categories: matlab code

Linear Phase Filters – why they are used

October 1, 2014 Leave a comment
%% Linear phase filters - preserve shape of a filtered signal
% This is the code used during a youtube video presentation dealing with linear phase filters
%      Search for linear phase at http://youtube.com/ddorran
%
% Code available from https://dadorran.wordpress.com
% 
close all ; clear all; clc
fs = 100;
T = 1/fs; %sampling interval
N = 2000; %length of signal being synthesised
n = 0:N-1; %samples of the signal
t = n*T;
 
plot_range = [N/2-100:N/2+100];
%% synthesise a signal
x = cos(2*pi*10*t) + 0.5*cos(2*pi*20*t + 1.4); 
subplot(2,1,1);
plot(t(plot_range),x(plot_range))
xlabel('Time (seconds)');
ylabel('Amplitude')
title('Synthesised Signals') 
axis tight
 
% Add some noise
ns = randn(1,length(x)+100)*2;
    %filter the noise to synthesise band limited noise
[b a] = butter(5, [0.28 0.33],'bandpass');
ns_filtered = filter(b,a,ns);
    %add noise to clean signal
x_ns = x +ns_filtered(end-length(x)+1:end);
hold on
noisy_x = plot(t(plot_range), x_ns(plot_range),'r');
legend('clean signal', 'noisy signal')
 
%% Plot frequency Content of Noisy Signal
subplot(2,1,2)
X_ns = fft(x_ns);
fax = [0:N-1]/(N/2); % normalised frequency axis
plot(fax(1:N/2), abs(X_ns(1:N/2))/(N/2)) ; %plot first half of spectrum
xlabel('frequency ( x \pi rads/sample)')
ylabel('Magnitude')
title('Magnitude Spectrum of Noisy Signal')
 

%% Filter out the noise using an IIR filter (non-linear phase)
[b_iir a_iir] = cheby1(10, 0.5, [0.27 0.34], 'stop');
y_iir = filter(b_iir,a_iir, x_ns);

[H_iir w] = freqz(b_iir,a_iir); %determine frequency response
subplot(2,1,2);
hold on
plot(w/pi, abs(H_iir),'r')
legend('|X(\omega)|','|H(\omega)|')

pause
Y_iir = fft(y_iir);
plot(fax(1:N/2), abs(Y_iir(1:N/2))/(N/2),'g') ; %plot first half of spectrum
legend('|X(\omega)|','|H(\omega)|','|Y(\omega)|')

pause
subplot(2,1,1)
non_linear_y = plot(t(plot_range),y_iir(plot_range),'g')
legend('clean signal', 'noisy signal','filtered signal')
pause 
set(noisy_x,'visible', 'off')

 
%% Examine the magnitude and phase response of the IIR filter
figure(2)
subplot(2,1,1)
plot(w/pi,abs(H_iir))
xlabel('frequency ( x \pi rads/sample)')
ylabel('Magnitude')
title('Magnitude Response of filter')
subplot(2,1,2)
plot(w/pi,angle(H_iir))
xlabel('frequency ( x \pi rads/sample)')
ylabel('Phase Shift')
title('Phase Response of filter')
 
%% Now filter using an FIR filter (with linear phase)
b_fir = fir1(100,  [0.27 0.34],'stop');
a_fir = 1;
y_fir = filter(b_fir,a_fir, x_ns);

figure(1)
subplot(2,1,1)
plot(t(plot_range),y_fir(plot_range),'k')
legend('clean signal', 'noisy signal','filtered signal (non-linear)','filtered signal (linear)')

[H_fir, w ]= freqz(b_fir,a_fir);
subplot(2,1,2)
plot(w/pi, abs(H_fir),'k')
legend('|X(\omega)|','|H(\omega) Non-linear|','|Y(\omega)|','|H(\omega)| linear')


 
%% Compare the frequency responses of the two filter design approaches
figure(2)
subplot(2,1,1)
hold on
plot(w/pi,abs(H_fir),'g')
legend('non-linear filter','linear filter')
subplot(2,1,2)
hold on
plot(w/pi,angle(H_fir),'g')
legend('non-linear filter','linear filter')
pause

%% Why does linear phase preserve the shape??
close all
clear all; clc;
fs = 1000;
t = 0:1/fs:2;
x1 = cos(2*pi*3*t-pi/2);
x2 = cos(2*pi*5*t-(pi/2)/3*5);
 
pause
subplot(3,1,1)
plot(t,x1)
subplot(3,1,2)
plot(t,x2)
subplot(3,1,3)
plot(t,x1+x2,'g')
hold on

pitch/period tracking using autocorrelation

September 24, 2014 Leave a comment
%% 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)