Archive
Archive for the ‘matlab code’ Category
Negative Frequency GUI Code
March 16, 2022
Leave a comment
function negfgui(varargin)
% An interactive GUI which attempts to help students visualise negative frequencies
% by David Dorran, TU Dublin, March 2022
% reuses some code to track mouse movement from the zpgui function by Tom Krauss from Perdue University
% Youtube video explanation provided at https://youtu.be/Rxc_ypVdruw
global zero_locations diff_equ_text op_text_samps op_freq_pts op_points draw_plt disp_plt op_disp_plt freq_plt xyz_pts op_xyz_pts xy_pts ptr negfqui_fh freq_pts points text_samps b a y gain pole_vals zero_vals freq_resp_pts
buf_len = 50;
tail_len = 10;
max_points = 100000;
fft_len = 1024;
win_func = hamming(buf_len);
if nargin == 0
action = 'init';
else
action = varargin{1};
end
if(strcmp(class(action), 'matlab.ui.control.UIControl'))
zero_vals = str2num(zero_locations.String);
negfgui('update_freq_response');
end
switch action
case 'init'
zero_vals = [ -j j];
pole_vals = [0 0];
if(strcmp(class(negfqui_fh),'matlab.ui.Figure'))
if(~isvalid(negfqui_fh))
negfqui_fh = figure;
end
if(nargin ==0)
close(negfqui_fh);
negfqui_fh = figure;
end
else
negfqui_fh = figure;
end
set(negfqui_fh,'position',[ 1 41 1366 651]);
ptr = 1;
points = ones(1,max_points)*NaN;
op_points = ones(1,max_points)*NaN;
freq_plt = subplot(2,3,1);
fax = (([0:fft_len-1])/fft_len*2*pi)-pi;
yyaxis left
freq_pts = plot(fax,points(1:fft_len));
hold on
op_freq_pts = plot(fax,points(1:fft_len),'r');
xlim([-pi pi])
ylim([0 buf_len/2])
ylabel('Magnitude Content')
yyaxis right
freq_resp_pts = plot(fax,points(1:fft_len));
xlim([-pi pi])
ylim([0 1])
title('Double sided-spectrum');
xlabel('frequency (radians per sample)')
ylabel('Magnitude Response')
draw_plt = subplot(2,3,2);
xy_pts = plot([]);
for m =1:tail_len
hold on
xy_pts(m) = plot(points(m)+points(m)*j,'r.','MarkerSize',tail_len-m+1);
end
ylabel('Y')
xlabel('X')
xlim([-1 1]);
ylim([-1 1]);
disp_plt = subplot(2,3,3);
xyz_pts = plot3(real(points),imag(points),[0:length(points)-1]);
hold on
surf([0 0 ],[-1.1 1.1],[-1 length(points); -1 length(points)],'FaceAlpha',0.2)
surf([-1.1 1.1], [ 0 0 ],[ length(points) length(points); -1 -1],'FaceAlpha',0.2)
xlabel('X')
ylabel('Y')
zlabel('Samples')
xlim([-1 1]);
ylim([-1 1]);
zlim([0 buf_len])
set(gca,'CameraUpVector', [0 1 0 ])
grid on
set(draw_plt,'buttondownfcn', 'negfgui(''draw'')')
op_disp_plt = subplot(2,3,6);
op_xyz_pts = plot3(real(op_points),imag(op_points),[0:length(op_points)-1]);
hold on
surf([0 0 ],[-1.1 1.1],[-1 length(op_points); -1 length(op_points)],'FaceAlpha',0.2)
surf([-1.1 1.1], [ 0 0 ],[ length(op_points) length(op_points); -1 -1],'FaceAlpha',0.2)
xlabel('X')
ylabel('Y')
zlabel('Samples')
xlim([-1 1]);
ylim([-1 1]);
zlim([0 buf_len])
set(gca,'CameraUpVector', [0 1 0 ])
grid on
text_samps = annotation('textbox','position',[.0500 0.48 0.85 0.05],'LineStyle','none','FontSize',10);
op_text_samps = annotation('textbox','position',[.0500 0 0.85 0.05],'LineStyle','none','FontSize',10);
diff_equ_text = annotation('textbox',...
[0.34 0.337941628264209 0.285237188872621 0.0844854070660522],...
'String',[''],...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none');
annotation('textbox',...
[0.34 0.2 0.285237188872621 0.0844854070660522],...
'String',['Zero Locations:'],...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none');
zero_locations = uicontrol('style','edit','units','normalized','FontSize',14,'position', [0.34 0.17 0.2 0.05]);
update_button = uicontrol('style','pushbutton', 'string','update locations','units','normalized','position',[0.34 0.12 0.2 0.05],'callback',{@negfgui});
negfgui('update_freq_response')
case 'update_freq_response'
b = poly(zero_vals);
a = poly(pole_vals);
H = freqz(b,a,fft_len,'whole');
gain = max(abs(H));
freq_resp_pts.YData= fftshift(abs(H)/gain); %ensure max gain is 1
b = b/gain;
subplot(2,3,4)
zplane(b,a);
title('Pole-Zero Plot')
try
str = create_difference_equation(b,1,'tex');
diff_equ_text.String = {['Difference Equation:'], str};
catch
diff_equ_text.String = {['b coefficients: ' num2str(b)], ['a coefficients:' num2str(a)]};
end
zero_locations.String = num2str(zero_vals);
case 'draw'
set(gcf,'userdata','')
set(gcf,'windowbuttonmotionfcn','set(gcf,''userdata'',''motion'')')
set(gcf,'windowbuttonupfcn','set(gcf,''userdata'',''up'')')
done = 0;
while ~done
waitfor(gcf,'userdata')
switch get(gcf,'userdata')
case 'motion'
pt = get(draw_plt,'currentpoint');
if(abs(pt(1,1))>1)
pt(1,1) = 1*sign(pt(1,1));
end
if(abs(pt(1,2))>1)
pt(1,2) = 1*sign(pt(1,2));
end
points(ptr) = pt(1,1)+pt(1,2)*j;
xyz_pts.XData(ptr) = real(points(ptr));
xyz_pts.YData(ptr) = imag(points(ptr));
%filter the signal (it's not perfect!)
try
op_points(ptr) = 0;
for k= 1:length(b)
op_points(ptr) = op_points(ptr)+ b(k)*points(ptr-k+1);
end
catch
end
op_xyz_pts.XData(ptr) = real(op_points(ptr));
op_xyz_pts.YData(ptr) = imag(op_points(ptr));
if(ptr>buf_len)
disp_plt.ZLim = [ptr-buf_len ptr];
op_disp_plt.ZLim = [ptr-buf_len ptr];
freq_pts.YData = abs(fftshift(fft(points(ptr-buf_len+1:ptr).*win_func', fft_len)));
op_freq_pts.YData = abs(fftshift(fft(op_points(ptr-buf_len+1:ptr).*win_func', fft_len)));
end
for m =1:tail_len
try
xy_pts(m).XData(1) = real(points(ptr-m+1));
xy_pts(m).YData(1) = imag(points(ptr-m+1));
text_samps.String = num2str(fliplr(round(points(ptr-m+1:ptr)*100)/100));
op_text_samps.String = num2str(fliplr(round(op_points(ptr-m+1:ptr)*100)/100));
catch
text_samps.String = num2str(fliplr(round(points(1:ptr)/100)*100));
end
end
ptr=ptr+1;
if(ptr> max_points)
negfgui('init')
end
case 'up'
done = 1;
end
set(gcf,'userdata','')
end
set(gcf,'windowbuttonmotionfcn','')
set(gcf,'windowbuttonupfcn','')
end
Categories: matlab code, Uncategorized, youtube demo code
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))
Categories: matlab code, youtube demo code
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
Categories: matlab code
LTspice Matlab, run LTspice matlab
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
Categories: matlab code, youtube demo code
pitch/period tracking using autocorrelation
September 24, 2014
Leave a comment
%% Using Autocorrelation to track the local period of a signal
% This code is used as part of a youtube video demonstration
% See http://youtube.com/ddorran
%
% Code available at https://dadorran.wordpress.com
%
% The following wav file can be downloaded from
% https://www.dropbox.com/s/3y25abf1xuqpizj/speech_demo.wav
%% speech analysis example
[ip fs] = wavread('speech_demo.wav');
max_expected_period = round(1/50*fs);
min_expected_period = round(1/200*fs);
frame_len = 2*max_expected_period;
for k = 1 : length(ip)/frame_len -1;
range = (k-1)*frame_len + 1:k*frame_len;
frame = ip(range);
%show the input in blue and the selected frame in red
plot(ip);
set(gca, 'xtick',[],'position',[ 0.05 0.82 0.91 0.13])
hold on;
temp_sig = ones(size(ip))*NaN;
temp_sig(range) = frame;
plot(temp_sig,'r');
hold off
%use xcorr to determine the local period of the frame
[rxx lag] = xcorr(frame, frame);
subplot(3,1,3)
plot(lag, rxx,'r')
rxx(find(rxx < 0)) = 0; %set any negative correlation values to zero
center_peak_width = find(rxx(frame_len:end) == 0 ,1); %find first zero after center
%center of rxx is located at length(frame)+1
rxx(frame_len-center_peak_width : frame_len+center_peak_width ) = min(rxx);
% hold on
% plot(lag, rxx,'g');
% hold off
[max_val loc] = max(rxx);
period = abs(loc - length(frame)+1);
title(['Period estimate = ' num2str(period) 'samples (' num2str(fs/period) 'Hz)']);
set(gca, 'position', [ 0.05 0.07 0.91 0.25])
[max_val max_loc] = max(frame);
num_cycles_in_frame = ceil(frame_len/period);
test_start_positions = max_loc-(period*[-num_cycles_in_frame:num_cycles_in_frame]);
index = find(test_start_positions > 0,1, 'last');
start_position = test_start_positions(index);
colours = 'rg';
subplot(3,1,2)
plot(frame);
set(gca, 'position',[ 0.05 0.47 0.91 0.33])
pause
for g = 1 : num_cycles_in_frame
if(start_position+period*(g) <= frame_len && period > min_expected_period)
cycle_seg = ones(1, frame_len)*NaN;
cycle_seg(start_position+period*(g-1):start_position+period*(g)) =...
frame(start_position+period*(g-1):start_position+period*(g));
hold on
plot(cycle_seg,colours(mod(g, length(colours))+1)) %plot one of the available colors
hold off
end
end
pause
end
%% synthesise a periodic signal to use as a basic demo
fs = 500;
T = 1/fs;
N = 250; % desired length of signal
t = [0:N-1]*T; %time vector
f1 = 8; f2=f1*2;
x = sin(2*pi*f1*t-pi/2) + sin(2*pi*f2*t);
plot(t, x)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Synthesised Signal');
%% Determine the autocorrelation function
[rxx lags] = xcorr(x,x);
figure
plot(lags, rxx)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')
%% Illustrate the auto correlation process
%function available from https://dadorran.wordpress.com
illustrate_xcorr(x,x)
%% Identify most prominent peaks
% Most prominent peak will be at the center of the correlation function
first_peak_loc = length(x) + 1;
% Lots of possible ways to identify second prominent peak. Am going to use a crude approach
% relying on some assumed prior knowledge of the signal. Am going to assume
% that the signal has a minimum possible period of .06 seconds = 30 samples;
min_period_in_samples = 30;
half_min = min_period_in_samples/2 ;
seq = rxx;
seq(first_peak_loc-half_min: first_peak_loc+half_min) = min(seq);
plot(rxx,'rx');
hold on
plot(seq)
[max_val second_peak_loc] = max(seq);
period_in_samples = abs(second_peak_loc -first_peak_loc)
period = period_in_samples*T
fundamental_frequency = 1/period
%% Autocorrelation of a noisy signal
x2 = x + randn(1, length(x))*0.2;
plot(x2)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Noisy Synthesised Signal');
[rxx2 lags] = xcorr(x2,x2);
figure
plot(lags, rxx2)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')
%% Autocorrelation technique can be problematic!
% Consider the following signal
f1 = 8; f2=f1*2;
x3 = sin(2*pi*f1*t) + 5*sin(2*pi*f2*t);
plot(t, x3)
ylabel('Amplitude')
xlabel('Time (seconds)')
title('Synthesised Signal');
[rxx3 lags] = xcorr(x3,x3,'unbiased');
figure
plot(lags, rxx3)
xlabel('Lag')
ylabel('Correlation Measure')
title('Auto-correlation Function')
seq = rxx3;
seq(first_peak_loc-half_min: first_peak_loc+half_min) = min(seq);
plot(seq)
[max_val second_peak_loc] = max(seq);
period_in_samples = abs(second_peak_loc -first_peak_loc)
Categories: matlab code, youtube demo code
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
Categories: matlab code, youtube demo code
requantise
September 17, 2014
Leave a comment
% 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
June 2, 2014
Leave a comment
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 &gt; 1
signal = signal';
end;
[r, c] = size(signal);
if r &gt; 1
disp('Note :only works on mono signals');
synthSignal = [];
return
end;
% add a few zeros to stop the algorithm failing
zpad = zeros(1, 44100/4);
signal = [signal, zpad];
if nargin &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;
Categories: audio processing, matlab code
Audio time-scale modification – VSOLA algorithm in matlab
June 2, 2014
Leave a comment
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));
ip = [ip zpad];
%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_seg = linear_cross_fade(seg_1, ip_frame, opt_overlap_length);
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-|
function op = linear_cross_fade(seg_1, seg_2, overlap_length, cross_fade_duration)
error(nargchk(3,4,nargin));
if nargin < 4
cross_fade_duration = overlap_length;
end
if cross_fade_duration > overlap_length
cross_fade_duration = overlap_length;
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)
linear_ramp1 = (cross_fade_duration-1:-1:0) ./(cross_fade_duration-1);
linear_ramp2 = (0:1:cross_fade_duration-1) ./(cross_fade_duration-1);
overlapping_region = (seg_1(length(seg_1)-cross_fade_duration+1:length(seg_1)).*linear_ramp1) + (seg_2(1: cross_fade_duration).*linear_ramp2);
end;
op = [seg_1(1:length(seg_1)-cross_fade_duration) ,overlapping_region , seg_2(cross_fade_duration+1:length(seg_2))];
% END of linear_cross_fade function
Categories: audio processing, matlab code, Uncategorized
dadorran 