## Linear Phase Filters – why they are used

%% 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
%
%
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

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);
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
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))
ylabel('Magnitude')
title('Magnitude Response of filter')
subplot(2,1,2)
plot(w/pi,angle(H_iir))
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

%% Using Autocorrelation to track the local period of a signal
% This code is used as part of a youtube video demonstration
%
%
%       https://www.dropbox.com/s/3y25abf1xuqpizj/speech_demo.wav
%% speech analysis example

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
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)



## illustrate_xcorr – code for cross correlation demos

September 24, 2014 1 comment
% 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

Categories: matlab code

## Audio Time Scale Modification – Phase Vocoder Implementation in Matlab

function synthSignal = pl_phaseVocoder_variable_analysis_hop(signal, tsm_factor, winSamps)
% A phase locked vocoder time-scale modification algorithm based upon Jean
% Laroche's 1999 work. The synthesis hop is fixed at one quarter the analysis window (hanning) while the analysis hop is scaled by the time scale factor; this results in more FFT computations than the bonada 2000 approach for slowing down - but provides smoother transitions frame to frame of the magnitude spectrums and, in general, a better quality of output. This code started
% out as the code provided by Tae Hong Park, but has changed significantly
% over the years
%
% David Dorran, Audio Research Group, Dublin Institute of Technology
% david.dorran@dit.ie
% http://eleceng.dit.ie/dorran
% http://eleceng.dit.ie/arg
%

% make sure input is mono and transpose if necessary
[r, c] = size(signal);
if r &amp;gt; 1
signal = signal';
end;
[r, c] = size(signal);
if r &amp;gt; 1
disp('Note :only works on mono signals');
synthSignal = [];
return
end;

% add a few zeros to stop the algorithm failing
if nargin &amp;lt; 3
winSamps = 2048;
end

winSampsPow2 = winSamps;
synHopSamps = winSampsPow2/4;
anHopSamps = round(synHopSamps/tsm_factor);

win = hanning(winSampsPow2);

X = specgram(signal, winSampsPow2, 100,win, winSampsPow2 - anHopSamps);

moduli = abs(X);
phases = angle(X);

[numBins, numFrames ] = size(phases);

syn_phases = zeros(numBins, numFrames); % a holder for synthesis phases

twoPi   = 2*pi;
omega   = twoPi * anHopSamps * [0:numBins-1]'/numBins/2; %the expected phase hop per frame

syn_phases(:,1) = phases(:,1) .* ( synHopSamps/ anHopSamps);

for idx =  2: numFrames
ddx = idx - 1;
deltaPhi = princarg(phases(:,idx) - phases(:,ddx) -omega); %calculate priciple determination of the hetrodyned phase increment
phaseInc = (omega+deltaPhi)/anHopSamps; % phase increment per sample
%locate the peaks
pk_indices = [];
pk_indices  =  locate_peaks(moduli(:,idx));
if(~length(pk_indices))
pk_indices = [1 10 12]; % just in case an odd situation is encountered  e.g. a sequence of zeros
end
%update phase of each peak channel using the phase propagation formula
syn_phases(pk_indices,idx)    = syn_phases(pk_indices,ddx)+synHopSamps*phaseInc(pk_indices); %synthesis phase
%update phase of channels in region of influence
% first calculate angle of rotation
rotation_angles = syn_phases(pk_indices,idx) - phases(pk_indices,idx);
start_point = 1; %initialize the starting point of the region of influence

for k = 1: length(pk_indices) -1
peak = pk_indices(k);
angle_rotation  = rotation_angles(k);
next_peak = pk_indices(k+1);
end_point = round((peak + next_peak)/2);
ri_indices = [start_point : peak-1, (peak+1) : end_point]; %indices of the region of influence
syn_phases(ri_indices,idx) = angle_rotation + phases(ri_indices, idx);
start_point = end_point + 1;
end;
end;

%Make sure that the LHS and RHS of the DFT's of the synthesis frames are a
%complex conjuget of each other
Z = moduli.*exp(i*syn_phases);
Z = Z(1:(numBins),:);
conj_Z = conj(flipud(Z(2:size(Z,1) -1,:)));
Z = [Z;conj_Z];

synthSignal = zeros(round(length(signal)*tsm_factor+length(win)), 1);

curStart = 1;
for idx = 1:numFrames-1
curEnd   = curStart + length(win) - 1;
rIFFT    = real(ifft(Z(:,idx)));
synthSignal([curStart:curEnd]) = synthSignal([curStart:curEnd]) + rIFFT.*win;
curStart = curStart + synHopSamps;
end

%--------------------------------------------------------------------------
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;

%--------------------------------------------------------------------------
function Phase = princarg(Phasein)
two_pi = 2*pi;
a = Phasein/two_pi;
k = round(a);
Phase = Phasein-k*two_pi;



## Audio time-scale modification – VSOLA algorithm in matlab

function op = vsola(ip, tsm_factor, P)
% Implementation of the variable parameter  synchronised overlap add VSOLA algorithm
% This implementation makes use of the standard SOLA algorithm (Rocus and Wilgus, ICASSP 1986) and some efficient paramter settings for the SOLA algorithm (Dorran, Lawlor and Coyle, ICASSP 2003) and (Dorran, Lawlor and Coyle, DAFX 2003)
% Given an input, ip, the longest likely pitch period, P, and  and a time scale modification factor, tsm_factor, return an output, op, that is a time scaled version of the input. The synthesis_overlap_size is the length, in samples of the lowest likely fundametal period of the input.
% for speech synthesis_overlap_size is set to (16/1000)*fs samples and for music synthesis_overlap_size is typically set to (20/1000)*fs samples
%
% David Dorran, Audio Research Group, Dublin Institute of Technology
% david.dorran@dit.ie
% http://eleceng.dit.ie/dorran
% http://eleceng.dit.ie/arg
%

% make sure input is mono and transpose if necessary
[r, c] = size(ip);
if r > 1
ip = ip';
end;
[r, c] = size(ip);
if r > 1
disp('Note :only works on mono signals');
op = [];
return
end;

% initialize the values of analysis_frame_length, analysis_window_offset, synthesis_frame_offset and length_of_overlap
desired_tsm_len = round(length(ip)*tsm_factor);
P = round(P); %longest likely pitch period in samples
Lmax = round(P * 1.5);% found this to a reasonable value for the Lmax- Lmax is the duration over which the correlation function is applied
stationary_length = (P * 1.5); % found this to a reasonable value for the stationary length - This is the max duration that could be discarded/repeated

analysis_window_offset = round((stationary_length - P)/abs(tsm_factor - 1)); % this equation was derived.
synthesis_window_offset = round(analysis_window_offset * tsm_factor);
analysis_frame_length = round(Lmax + synthesis_window_offset);
number_of_analysis_frames = floor((length(ip)- analysis_frame_length)/analysis_window_offset);

if number_of_analysis_frames < 2 %not much time-scaling being done just return the input
op = ip;
return;
end;

%the next two lines just ensure that the last frame finishes at the very end of the signal (not essential)
zpad = zeros(1, (number_of_analysis_frames*analysis_window_offset) + analysis_frame_length - length(ip));

%initialize the output
op = zeros(1, desired_tsm_len);
%initialize the first output frame
op(1 : analysis_frame_length) = ip(1 : analysis_frame_length);

min_overlap = round(Lmax - P); %ensure that there is some minimum overlap
count = 0;

% Loop for the 2nd analysis frame to the number_of_analysis_frames
for m = 1 : number_of_analysis_frames

%grab the mth input frame
ip_frame = ip(analysis_window_offset * m : (analysis_window_offset * m) + analysis_frame_length - 1);

%grab the maximum overlapping segments from the inout frame and the current output
seg_1 = op(round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax : round(synthesis_window_offset*(m-1))+analysis_frame_length -1);
seg_2 = ip_frame(1: Lmax);

%compute the correlation of these segments
correlation   = xcorr(seg_2, seg_1,'coeff');

%Find the best point to overlap (opt_overlap_length) making sure not to exceed the maximum or go below the minimum overlap.
correlation(length(correlation) - Lmax -1: length(correlation)) = -100;
correlation(1: min_overlap) = -100;
[max_correlation, opt_overlap_length] = max(correlation);

if(max_correlation == 0)
opt_overlap_length = Lmax;
end;
%     offset = Lmax - opt_overlap_length;
%     if ((offset + analysis_window_offset -  synthesis_window_offset) >= 0 & (offset + analysis_window_offset -  synthesis_window_offset) <= P)
%         count = count +1;
%     end;

% append mth analysis frame to the current synthesised output using a linear cross fade
ov_len =(round(synthesis_window_offset*m)+analysis_frame_length) - (round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax) + 1;
ov_seg = ov_seg(1:ov_len);
op(round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax: round(synthesis_window_offset*m)+analysis_frame_length) = ov_seg;

end; % end of for loop

% linear cross fade the first segment with the second segment given a certain amount of overlap
% |----------seg1----------|
%			   |---------------seg2-----------------------|
%              |--overlap-|

error(nargchk(3,4,nargin));
if nargin < 4
end
end
% overlap the end of seg_1 with the start of seg_2 using a linear cross-fade
if (length(seg_1) < overlap_length),
seg_2 = seg_2(overlap_length - length(seg_1) + 1: length(seg_2));
overlap_length = length(seg_1);
end; % end of if statement

if (length(seg_2) < overlap_length),
seg_1 = seg_1(length(seg_1) - (overlap_length - overlap_length): length(seg_1));
overlap_length = length(seg_2);
end; % end of if statement

overlapping_region = zeros(1, cross_fade_duration); % initialize the overlapping region
seg_1 = seg_1(1: length(seg_1) - (overlap_length - cross_fade_duration));

op = zeros(1, length(seg_1) + length(seg_2) - cross_fade_duration);

if(overlap_length ~= 1)
end;



## pic18f4620 xc8 compiler bundle

Here’s a small (5.2 MB) bundle that I use to compile and program c code on my pic18f4620 using a pickit2 on my PC. I’ve only tested this on XP.

It doesn’t have all the features of a full blown IDE like MPLABX (no stepping through code I’m afraid) but is a nice concise package that encourages the use of command line scripting.

Getting Started

• Extract the files to a folder on your PC or USB key The extracted folder contains the xc8 compiler files (modified to work with the pic18f4620 only), the pk2cmd software tool which is used to communicate with the PICkit2, and Notepad++, which can be used to write and modify code.

To create a new project do the following:

1. First create a new folder inside the “robosumo_files” folder and copy the build.bat file into that folder.
2.  Run the Notepad++  editor and paste the example code (flashing LED) given below. Save the code in the new folder as a file called “main.c” (don’t include quotes in the filename).
3. Make sure the PICkit2 programmer is correctly connected to the PIC18f4620 and finally double click the build.bat file.
4. If there are no errors in the code and the hardware is set up correctly then the PIC18f4620 will be programmed.

## Example Code – Flashing LED

This is a pared back example program for the PIC18F4620. It blinks an LED connected to pin 19 (RD0). See this circuit

//
// PIC18F4620 example program
// Written by Ted Burke (ted.burke@dit.ie)
// Last update 28-2-2013
//

#include &lt;xc.h&gt;

#pragma config OSC=INTIO67,MCLRE=OFF,WDT=OFF,LVP=OFF,BOREN=OFF

int main(void)
{
// Make RD0 a digital output
TRISD = 0b11111110;

while(1)
{
LATDbits.LATD0 = 1; // Set pin RD0 high
_delay(125000);     // 0.5 second delay
LATDbits.LATD0 = 0; // Set pin RD0 low
_delay(125000);     // 0.5 second delay
}
}

Categories: pic18f, Uncategorized

## pic18f4620 assembly bundle

Here’s a small (3.6 MB) bundle that I use to program my pic18f4620 using a pickit2 with assembly code on my PC. I’ve only tested this on XP.

It doesn’t have all the features of a full blown IDE like MPLABX (no stepping through code I’m afraid) but is a nice concise package that encourages the use of command line scripting.

The bundle also includes some example code.

### Getting Started

The extracted folder contains the mpasmx assembler files, the pk2cmd software tool which is used to communicate with the PICkit2, and Notepad++ text editor, which can be used to write and modify code.

The bundle includes 6 example projects.

To create a new project do the following:

1. First create a new folder inside the “pic18f4620_asm” folder and copy the build.bat file into that folder.
2. Run the Notepad++ editor and paste the example code (flashing LED) given below. Save the code in the new folder as a file called “main.asm” (don’t include quotes in the filename).
3. Make sure the PICkit2 programmer is correctly connected to the PIC18f4620 and finally double click the build.bat file.
4. If there are no errors in the code and the hardware is set up correctly then the PIC18f4620 will be programmed.

Example Code to work with this circuit
This code will flash the LED connected to pin 19 (RD0)

;configure the  assembler directive 'list' so as to set processor to 18f4620 and set the radix used for data expressions to decimal (can be HEX|DEC|OCT)
list p=18f4620, r=DEC
#include <p18f4620.inc>
; configure the micro so that the watchdog timer is off, low-voltage programming is off, master clear is off and the clock works off the internal oscillator
config WDT=OFF, LVP=OFF, MCLRE=OFF, OSC=INTIO67
;The org directive tells the compiler where to position the code in memory
org 0x0000 ;The following code will be programmed in reset address location i.e. This is where the micro jumps to on reset

;--------------------------------------------------------------------------
; Main Program
;--------------------------------------------------------------------------

org 0x0100
Main
clrf TRISD ; All PORT D IO pins are configured as outputs

; Set up timer 0 so it can be used to control flash rate
bcf T0CON, T0CS; use internal instruction cycle clock
bcf T0CON, T08BIT; use 16 bit mode
bcf INTCON, TMR0IF; reset timer overflow flag
bsf T0CON, TMR0ON; turn on timer 0

loop
; check if the timer0 overflow flag has been set
btfss INTCON, TMR0IF
goto loop
bcf INTCON, TMR0IF; reset timer overflow flag

;invert LATD0
movlw 0x1
XORWF LATD,F

goto loop
END


Categories: pic18f, Uncategorized

## Heartrate (BPM) Example Matlab Code

This is the code I used in my youtube video at http://youtu.be/3tdumuwHgxc

% program to determine the BPM of an ECG signal

% count the dominant peaks in the signal (these correspond to heart beats)
% - peaks are defined to be sampels greater than their two nearest neighbours and greater than 1

beat_count = 0;
for k = 2 : length(sig)-1
if(sig(k) &gt; sig(k-1) &amp; sig(k) &gt; sig(k+1) &amp; sig(k) &gt; 1)
%k
%disp('Prominant peak found');
beat_count = beat_count + 1;
end
end

% Divide the beats counted by the signal duration (in minutes)
fs = 100;
N = length(sig);
duration_in_seconds = N/fs;
duration_in_minutes = duration_in_seconds/60;
BPM_avg = beat_count/duration_in_minutes;


## pic18f assembly example 6 – discrete filter implementation

;-----------------------------------------------------------------------------------------------------------------------------
; This assembly code sets up a timer interrupt to be triggered every millisecond on the pic18f4620.
; When the interrupt service routine is called an ADC on channel 0 is performed and this value  is used as an
;    input to the discrete system given by the difference equation y[n] = x[n]-0.3x[n-1]+0.3y[n-1]-0.75y[n-1]
;    The output of the discrete system is then sent to a  DAC (MCP4921) via the SPI protocol.
;
; Test cicruit used is same as at http://eleceng.dit.ie/daq_lite/build.html
;
; This code is just an educational example and no attempt has been made to optimise it in any way. Use it whatever way you want but I don't accept any responsibility for any problems this code might cause you.
;
; NOTE 1: Only the most significant 8 bits of the 10-bit ADC are used. Therefore the
;    the ADC values used are in the range of 0-255.
; NOTE 2: When implementing the discrete system the ADC values are interpreted as a signed signal which varies between 0.9922 (=127/128) and -1.
;    An ADC value of 255 is interpreted as a signal value of 0.9922; an ADC value of 0 is interpreted as a signal value of -1 ;
;    an ADC value of 128 is interpreted as a signal value of 0; An ADC value of 127 is interpreted as a signal value of -0.0078 (1/128)
; NOTE 3: The pic18f does not have a floating point unit and signed decimal values are represented in and 8-bit, fixed point, 2's complement format.
;    A signal value of 0.9922 is stored as a binary value of '01111111' in memory
;    A signal value of -1 is stored as a binary value of '10000000'  in memory
;    A signal value of 0.5 is stored as binary value of '01000000' in memory
;    A signal value of -0.5 is stored as binary value of '11000000' in memory
;    A signal value of 0.75 is stored as binary value of '01100000' in memory
;    A signal value of -0.75 is stored as binary value of '10100000' in memory
;    A signal value of 0.0078 is stored as binary value of '00000001' in memory
;    A signal value of -0.0078 is stored as binary value of '11111111' in memory
;
; by David Dorran (https://dadorran.wordpress.com) May 2014
;-----------------------------------------------------------------------------------------------------------------------------

;configure the  assembler directive 'list' so as to set processor to 18f4620 and set the radix used for data expressions to decimal (can be HEX|DEC|OCT)
list p=18f4620, r=DEC
#include <p18f4620.inc>
; configure the micro so that the watchdog timer is off, low-voltage programming is off, master clear is off and the clock works off the internal oscillator
config WDT=OFF, LVP=OFF, MCLRE=OFF, OSC=INTIO67
;The org directive tells the compiler where to position the code in memory
org 0x0000 ;The following code will be programmed in reset address location i.e. This is where the micro jumps to on reset

;NOTE: When an interrupt is triggered the micro jumps to 0x0008 for high priority and 0x0018 for low priority.
org 0x0008
goto hi_isr ;once a high priority interrupt is triggered jump to hi_isr
org 0x0018
goto low_isr ;once a low priority interrupt is triggered jump to low_isr

;--------------------------------------------------------------------------
; Main Program -------------------------------------------------------------
;--------------------------------------------------------------------------
org 0x0100 ; the following code is placed at address 0x0100 in program memory (Flash memory)
Main
; store the filter variables in a contiguous block in RAM starting at address 0x000
cblock 0x000
; y_1 = y[n-1]; y_2 = y[n-2]; x_1 = x[n-1]; x_2 = x[n-2]
x, y,a1,a2,b0,b1,b2, var1,var2, x_1, y_1, x_2, y_2
endc

; set up b and a coefficients of the discrete system (y[n] = x[n]-0.3x[n-1]+0.3y[n-1]-0.75y[n-1])
movlw .127
movwf b0 ; approx = 1 (exact value is 0.9922)
movlw B'11011010'
movwf b1 ; approx = -0.3 (exact value is -0.3047)
movlw B'00000000'
movwf b2 ; = 0
movlw B'00100110'
movwf a1 ; = approx 0.3 (exact value is 0.3047)
movlw B'10100000'
movwf a2 ; = -0.75

; intialise the 'past value' variables to zero
clrf y_1
clrf y_2
clrf x_1
clrf x_2

call micro_config ; configure pins,timers, ADC etc.

main_loop
nop ; Just wait for a timer interrupt to be triggered
goto main_loop

;----------------------------------------------------------------------------------------------------------
;-------multiply accumulate - multiply var1 by var2 and accumulate result in y   --------------------------
;------   this routine is used a lot to filter the signal -------------------------------------------------
;----------------------------------------------------------------------------------------------------------
mac
; pic18f only has 8X8 unsigned multiply - the following code analyses the result of an 8X8 unsigned multiply to produce a signed result
movf var1, W
mulwf var2 ; var1 * var2 -> PRODH:PRODL (var1 and var interpreted as     btfsc var2, 7 ; Test Sign Bit (NOTE: if var2 is negative (2's complement representation) then its unsigned interpretation is 256-var2 and then the result stored in PROD is var1*(256-var2) = var1*256-var1*var2
btfsc var2, 7 ; Test Sign Bit
subwf PRODH, F ; PRODH = PRODH - var1 (equivalent to PROD = PROD - var1*256)
movf var2, W
btfsc var1, 7 ;Test Sign Bit
subwf PRODH, F ; PRODH = PRODH - var2

; The two most significant bits of PROD are sign bits so shift PRODH to the left and grab the most significant bit of PRODL to get the 8 most significant useful bits.
rlncf PRODH,0 ; W contains PRODH shifted to the left by one
movwf var2 ; using var2 as it is available (not for any other reason)
bsf var2, 0
BTFSS PRODL,7 ; set the most significant bit of PRODL as the least signifcant in W
bcf var2, 0
movf var2,W

addwf y,f ; accumulate the result of multiplication in y. Note not handling overflow. Would also be better to accumulate those 7 bits discarded after the multiplication and round afterwards for a more accurate result

return

;----------------------------------------------------------------------------------------------------------
;-----High Priority Interrupt Service Routine -------------------------------------------------------------
;----------------------------------------------------------------------------------------------------------
org 0x0200
hi_isr
btfss INTCON, TMR0IF
goto end_int
bcf INTCON, TMR0IF; reset interrupt
;An interrupt has been set up to be generated whenever TMR0 goes from FFFFh to 0000h i.e. on overflow.
;Would like the timer overflow flag triggered every 1ms so will set timer to 0xffff - 1000 = 0xFC17
movlw 0x17
movwf TMR0L
movlw 0xFC
movwf TMR0H

clrf y

movf y,W

addlw -128 ; account for dc offset. The input values from ADRESH will be in range 0-255. These values represent a signal that varies betwen 1 and -1: the ADC value of 255 is mapped to 1 and ADC value of 0 is mapped to -1; ADC value of 128 is mapped to 0.

movwf x
movwf var1
movff b0,var2
call mac ; multiply var1 by var2 and accumulate the result in y (mac - multiply accumulate)

movff x_1, var1
movff b1, var2
call mac;

movff x_2, var1
movff b2, var2
call mac;

movff y_1, var1
movff a1, var2
call mac;

movff y_2, var1
movff a2, var2
call mac

; update past values of the inputs and outputs to use on next iteration of the discrete system
movff y_1, y_2
movff y, y_1
movff x_1, x_2
movff x, x_1

movf y,W
addlw .128 ; account for dc offset
movwf y
goto spi_write ;write contents of y to SPI device
end_int
retfie ;return and reset interrupts

;----------------------------------------------------------------------------------------------------------
;-------configure pins, ADC,timers and interrupts on the pic18f4620    -----------------------------------
;----------------------------------------------------------------------------------------------------------
micro_config
; configure LATD0-LATD3 as outputs and RD4-RD7 as inputs
movlw 0xf0
movwf TRISD
clrf LATD ; turn off all LATD output pins

; SPI configuration ------------------------------------------------------------
clrf TRISC
;bcf TRISC,3 ; pin 3 on PORTC used as SPI clock and should be configured as an output to operate in master mode
;bcf TRISC,5 ; pin 5 on PORTC used for SPI serial data out (SDO) and should be configured as an output
movlw B'10110000' ; set up SPI control register see pg 163 for details
movwf SSPCON1;
bsf SSPSTAT,CKE; data transmitted on rising edge
; END SPI config ----------------------------------------------------------

; Set clock frequency (section 2 of the PIC18F4620 Data Sheet)
; Set Fosc = 8MHz, which gives Tcy = 0.5us
movlw B'01110000'
movwf OSCCON

; all channels set up for ADC - no particular reason why
movlw B'00000001'

; set tad so that capacitor on ADC can fully charge - see pg129
movlw B'00100010'

clrf ADCON0; select analog channel 0;

bcf INTCON, GIE ;Disable global interrupts
bcf T0CON, TMR0ON; turn off timer 0
bsf INTCON, TMR0IE; Enable TIMER0 interupt

;set up 1:2 prescaler so that TMR0 SFR is incremented every 2 Tcy  i.e. 1us
bcf T0CON,0
bcf T0CON,1
bcf T0CON,2

bcf T0CON, T0CS; use internal instruction cycle clock
bcf T0CON, T08BIT; use 16 bit mode
bcf T0CON, PSA; turn on prescaler

; setup so that timer0 overflow is triggered quickly initially
; after the first trigger then set it up so that it is triggered at the rate you want in the high_isr routine
movlw 0xFF
movwf TMR0L
movlw 0xFF
movwf TMR0H

bsf INTCON2, TMR0IP; ; Set the timer0 interrupt up as high priority
bsf INTCON, GIE ;Enable global interrupts
bsf T0CON, TMR0ON; turn on timer 0
return

;----------------------------------------------------------------------------------------------------------
;-----WRITE contents in y to DAC (MCP4921) using SPI ------------------------------------------------------
;----------------------------------------------------------------------------------------------------------
spi_write
;two bytes of data sent to spi daq (a 12-bit dac). First nibble of first byte is dac config settings; second nibble is
; the most significant 4 bits of 12bit numerical value. The second byte contains the 8 least signifcant bits of the 12-bit value

bcf LATD, LATD1; // Select SPI DAC chip

movf SSPBUF, W ;WREG reg = contents of SSPBUF (this clears BF) which will be set again after the data is transmitted
; load SSPBUFF with spi dac config data to be sent to SPI device
; First nibble of first byte sent to dac is dac config settings; second nibble is
; the most significant 4 bits of 12bit numerical value
swapf y, f
movlw 0x0f
andwf y, W ; W no holds 4 most signicant bits in lower nibble

iorlw 0x70 ; these are configuration bits for spi dac
movwf SSPBUF

spi_wait1
btfss SSPSTAT, BF ;Has data been received (transmit complete)?
goto spi_wait1 ;No

movf SSPBUF, W ;WREG reg = contents of SSPBUF (this clears BF) which will be set again after the data is transmitted
; load SSPBUFF with y data to be sent to SPI device

movlw 0xf0
andwf y, W ;
movwf SSPBUF

spi_wait2
btfss SSPSTAT, BF ;Has data been received (transmit complete)?
goto spi_wait2 ;No

bsf LATD, LATD1; // Select SPI DAC chip
goto end_int

;----------------------------------------------------------------------------------------------------------
;-----Low Priority Interrupt Service Routine -------------------------------------------------------------
;----------------------------------------------------------------------------------------------------------
org 0x0300
low_isr
nop ; put whatever you would like to do on low priority interrupt here
retfie ;return and reset interrupts

end ; End of ASM code

Categories: pic18f, Uncategorized