## Digital Signal Processing Notes

I’ve made some note relating to Digital Signal Processing and the Discrete Fourier Transform available at http://pzdsp.com/docs.

I intend to add more documents to this over time.

Categories: Uncategorized

## Code to extract music notation from Flowkey app

``````% Matlab code to extract music notation from flowkey screen grab video . See https://youtu.be/4oOMWkMRCCg

ind = 5*25;
im = frame(430:660 ,: ,:);
new_im = rgb2gray(im);
prev_x = sum(new_im); % used to find bes overlap of frames

while 1
ind = ind+5*25;
if(ind > v.NumFrames)
break
end
im = frame(430:660 ,: ,:);
img = rgb2gray(im); %conver to greyscale
x = sum(img);
xc =[]; %stores cross correlation result
seg2 = prev_x(end-600+1:end); %end of previous frame
for w = 1: 300
seg1 = x(1+w:600+w);
xc(w) = sum(seg1.*seg2)./(sqrt(sum(seg1.^2)*sum(seg2.^2))); %cross correlation
end
[val loc] = max(xc); % loc is where the max correlation occurs

[ r c d] = size(im);
new_seg =img(:,201:end);
[ r1 c1 d1] = size(new_seg);
new_im(:,end-600+1-loc+200:end-600-loc+200+c1) = new_seg;
prev_x = x;

end
img = new_im;
% find where the bars occur in the image
for k =1: length(img)
x(k) = sum(img(95:145,k)); %the bar lines cross are clear from rows 75 to 125
end
bars = find(x<5000); %obtained value of 5000 by visually inspecting plot(x)

op_image_max_width = 2000;
op_image_max_height = 1200;

overlap = 10;
start_pt = overlap+1;
prev_end_pt = bars(1);
op_img = img(1,1);
k = 1;
img_num = 1;
while(k<length(bars))
end_pt = bars(k);
if(end_pt-start_pt>op_image_max_width)
[rows cols ] = size(img(:,start_pt-overlap:prev_end_pt+overlap));
[rows1 cols1] = size(op_img);
if(rows+rows1 > op_image_max_height)
imwrite(op_img, ['notation_' num2str(img_num) '.jpg']);
op_img = img(1,1);

img_num = img_num + 1;
end
op_img(end+1:end+rows,1:cols,:) = img(:,start_pt-10:prev_end_pt+10);
start_pt = prev_end_pt;
end
prev_end_pt = end_pt;
k = k + 1;
end

imwrite(op_img, ['notation_' num2str(img_num) '.jpg']);
``````
Categories: Uncategorized

## Arduino sinuoidal waveform generator and signal capture

I have put together an example on how to use the Arduino nano as a signal generator and digital capture device. It runs a pretty low sampling rate 2 kHz but I thought it might be of some use. Check out a video explanation at  https://youtu.be/_31sQn1Cg9g.  Files used are available at http://pzdsp.com/shared/arduino_sin_gen.zip

Categories: Uncategorized

## Virtual Oscilloscope

I’ve created a virtual oscilloscope that hopefully will help students learn how to work with real ones! It doesn’t have too many features so students can focus on the fundamentals that appear on all scopes: scaling, positioning and triggering.

Check it out at pzdsp.com/elab

Categories: Uncategorized

## FFT based FIR filter design

```%% This script deals issues related with 'fft' based FIR filter design
% See youtube video at https://youtu.be/R4RpNG_Botk

% GO TO lines 115 to 133 if you want to skip the preamble

%% First show fft of finite number of impulse response samples  lie on the continuous spectrum
%Design a band reject filter using buit-in fir1 function
b = fir1(22, [0.3913 0.6522],'stop'); %cutoff frequency is (0.3, 0.6 x pi) radians/sample = 0.3,0.6 x fs/2;
h = b; % The non-zero values of the impulse response of any fir filter
%   is equal to the b coefficients
H = fft(h); % Frequency response bin values
bin_freqs = [0:length(H)-1]/length(H)*2*pi;

%% Create the continuous spectrum
% Going to do this three ways - just for the purpose of demonstration.
% The mathematical formula for the frequency response of an mth order FIR is:
%   H(w) = b0 + b1*e^-jw + b2*e^-2jw + b3*e^-2*jw + ... ... + bm*e^(-m*jw)
%
%   This equation comes from the fact that H(w) = H(z) when z= e^jw i.e. when H(z) is
%   evaluated around the unit circle (see https://youtu.be/esZ_6n-qHuU )
%  The transfer function H(z) = b0*Z^0+b1*z^-1+b2*z^-2+b3*z-3+ ... +bmz^-m

% We'll create a 'continuous' spectrum by evaluating this equation for a
% "large" number of frequencies, w.

% Two alternatives to creating a 'continuous' spectrum are to zero pad the
% b coefficients by a "large" amount prior to taking the fft; and to use the built in freqz
% function. I'll use all three approaches here for comparison purposes.

% METHOD 1 - zero-pad method
H_cont_1 = abs(fft([h zeros(1,10000)])); % the value of 10000 could be any relatively large number i.e. large enough to capture the spectral detail of the frequency response
w_cont_1 = [0:length(H_cont_1)-1]/length(H_cont_1)*2*pi;

% METHOD 2 - freqz method
[H_cont_2_complex w_cont_2] = freqz(b,1,length(H_cont_1), 'whole'); % the second parameter of freqz is the a coefficients
H_cont_2 = abs(H_cont_2_complex);

% METHOD 3 -evaluate H(z) around the unit circle method (see
%https://youtu.be/esZ_6n-qHuU)
H_cont_3_complex = [];
k = 0;
w_cont_3 = [0:length(H_cont_1)-1]/length(H_cont_1)*2*pi;
for w = w_cont_3
k = k + 1;
H_cont_3_complex(k) = b(1);
for ind = 2:length(b)
H_cont_3_complex(k) = H_cont_3_complex(k) + b(ind)*exp(-1i*w*(ind-1));
end
end

H_cont_3= abs(H_cont_3_complex);

%% plot 'continuous' frequency response and overlay 'sampled' bin values
figure
plot(w_cont_1, H_cont_1,'LineWidth', 6)
hold on
plot(w_cont_2, H_cont_2,'LineWidth', 3)
plot(w_cont_3, H_cont_3)
plot(bin_freqs, abs(H),'g.','MarkerSize', 16)
ylabel('Magnitude')
legend('Cont. Freq Resp Method 1', 'Cont. Freq Resp Method 2', 'Cont. Freq Resp Method 3', 'Bin Values')
title('Frequency response of filter using a built-in filter design technique')
set(gcf,'Position',[22 416 613 420]);

%% Now let's try and design a filter by starting with the desired frequency response
% First we'll look at the band pass filter
H_desired = [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0];
freq_bins_desired = [0:length(H_desired)-1]/length(H_desired)*2*pi;
b_ifft = ifft(H_desired); %b coefficients of newly designed filter
h_ifft = b_ifft;

% Now compute the continuous spectrum using any of the methods above. (Zero
% padding b_ifft prior to taking the fft is the most straightforward
% appraoch - I think!)
H_cont_desired = abs(fft([b_ifft zeros(1,1000)]));
w_cont_desired = [0:length(H_cont_desired)-1]/length(H_cont_desired)*2*pi;

figure
plot(w_cont_desired, H_cont_desired);
hold on
plot(freq_bins_desired, H_desired,'r.','MarkerSize', 12)

ylabel('Magnitude')
legend('Continuous Frequency Response','DFT bin values of desired frequency response')
title('Frequency response of band pass filter using ''inverse fft'' filter design technique')
set(gcf,'Position',[649   190   610   420]);

%%
% We'll now mimic the 'desired' filter shown in the video i.e. a band reject filter
H_desired = [1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1];
freq_bins_desired = [0:length(H_desired)-1]/length(H_desired)*2*pi;
b_ifft = ifft(H_desired);

% Now compute the continuous spectrum using any of the methods above. (Zero
% padding b_ifft prior to taking the fft is the most straighforward
% appraoch I think)
H_cont_desired = abs(fft([b_ifft zeros(1,1000)]));
w_cont_desired = [0:length(H_cont_desired)-1]/length(H_cont_desired)*2*pi;

figure
plot(w_cont_desired, H_cont_desired);
hold on
plot(freq_bins_desired, H_desired,'r.','MarkerSize', 12)

ylabel('Magnitude')
legend('Continuous Frequency Response','DFT bin values of desired frequency response')
title('Frequency response of band reject filter using ''inverse fft'' filter design technique')
set(gcf,'Position',[649   190   610   420]);

%%
% Now try using a 'practical' linear phase response rather than phase values of zero

H_desired = [1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1];
% An alternative way to create H_desired to ensure first half is a mirror
% image of the second half (for odd lengths only!)
% H_desired = [1 1 1 1 1 0 0 0 1 1 1 1];
% H_desired = [H_desired fliplr(H_desired(2:end))]

phase_diff = pi/length(H_desired)-pi;
phase_desired = [0:floor(length(H_desired)/2)]*phase_diff;
phase_desired = [phase_desired fliplr(phase_desired(2:end))*-1];

freq_bins_desired = [0:length(H_desired)-1]/length(H_desired)*2*pi;
h_ifft = ifft(H_desired.*exp(j*phase_desired));

b_ifft = h_ifft; %b coefficients of newly designed filter
%b_ifft = h_ifft.*hanning(length(h_ifft))'; %windowed b coefficients of newly designed filter

% Now compute the continuous spectrum using any of the methods above. (Zero
% padding b_ifft prior to taking the fft is the most straighforward
% approach - I think!)
H_cont_desired = abs(fft([b_ifft zeros(1,1000)]));
w_cont_desired = [0:length(H_cont_desired)-1]/length(H_cont_desired)*2*pi; %frequency axis

plot(w_cont_desired, H_cont_desired, 'Linewidth', 1);
hold on
plot(freq_bins_desired, H_desired,'r.','MarkerSize', 12)

ylabel('Magnitude')
legend('Continuous Frequency Response','DFT bin values of desired frequency response')
title('Frequency response of band reject filter using ''inverse fft'' filter design technique')
set(gcf,'Position',[649   190   610   420]);

%% Demonstration of filters being a sum of prototype bandpass filters
% H_acc is the accumulation of individual prototype bandpass filter

H_prot = zeros(size(H_desired)); %initialise
H_prot(1) = H_desired(1).*exp(j*phase_desired(1));
H_acc = fft(ifft(H_prot), length(H_cont_desired));

% loop through 11 prototype band pass filters
show_plots = 0;
if(show_plots)
figure
end
for k = 1: 11
prot_bins = [k+1 24-k] %pair of bins associated with positive and negative frequencies
H_prot = zeros(size(H_desired)); %initialise
H_prot(prot_bins) = H_desired(prot_bins).*exp(j*phase_desired(prot_bins));

H_cont_prot = fft(ifft(H_prot), length(H_cont_desired));
H_acc = H_acc + H_cont_prot; %accumulate each individual band pass prototype
if(show_plots)
plot(abs(H_acc))
hold on
plot(abs(H_cont_prot))
hold off
pause(1)
end
end

sum(abs(H_acc) - abs(H_cont_desired))

```

## DSP Foundations Notes

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.

Categories: Uncategorized

## interpLTspice

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

```
Categories: Uncategorized

## runLTspice demo

```%% Demonstration of the runLTspice function
%
% Demo video at http://youtu.be/ax9H4eKuZv4
%
%
%
% 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

```
Categories: Uncategorized

## runLTspice.m

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

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)
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))
else
ltspice_path = 'C:\Program Files\LTC\LTspiceXVII\XVIIx64.exe';
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 ];
path_needs_updating = 0;
while ~exist(exe_file)
disp(['The LTspice execuatble file could not be found. ']);
if(~path_needs_updating)
disp('Once the correct location is entered it will be recorded for future use of this function.')
disp('     Example location 1: C:\Program Files\LTC\LTspiceIV\scad3.exe')
disp('     Example location 2: C:\Program Files\LTC\LTspiceXVII\XVIIx64.exe')
end
disp(' ')
ltspice_path = input('Please enter the location where the executable exists => ','s');
ltspice_path= strtrim(ltspice_path);
if(ltspice_path(end) == '/')
ltspice_path(end) ='\';
end

exe_file = [ltspice_path ];
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

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

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
fs = 0;
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});
fs = varargin{1+(num_sources-1)*3+3};
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 LTsppice asc file exists. \nCopy the file %s to a folder you can write to and try again. ', asc_file));
end
k = 0;
trans_found = 0;
while( k < num_asc_file_lines - 1)
k = k + 1;
res = findstr('!.tran ', tline{k});
if(length(res))
trans_found = 1;
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
sim_duration
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
if(~trans_found)
error(['The simulation must be run in transient mode (rather than AC analysis for example). Select the Transient tab by first selecting Simulate->Edit Simulation Cmd in LTSpice'])
end
raw_data = LTspice2Matlab(tmp_raw_file);

if(fs)
varargout{1} = interpLTspice(varargout{1},fs);
end
%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 -ascii "' test_filename '"'];
%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
error( [ vars{k} ' not found! Check that you spelled the variable correctly. '])
data_details{1}.data{k} = [];
end
end
data_details{1}.time_vec = raw_info.time_vect;
end

end
```
Categories: matlab code

## locate_peaks matlab function

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