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)
xlabel('Frequency (radians per sample)')
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)

xlabel('Frequency (radians per sample)')
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)

xlabel('Frequency (radians per sample)')
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)

xlabel('Frequency (radians per sample)')
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

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.

Categories: Uncategorized

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;

Categories: Uncategorized

runLTspice demo

October 29, 2014 2 comments
%% 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

October 29, 2014 6 comments
% 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)) ltspice_path = fileread(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 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 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; res2 = regexp(tline{k}, '.tran (?.*?)(?(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
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);

varargout = grab_specified_data(vars2read, raw_data);
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

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