FFT based FIR filter design

July 19, 2019 Leave a comment
%% 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
%start with the DC prototype filter
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))

Advertisements

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 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
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\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_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