function varargout = pole_zero_place(varargin)
% function pole_zero_place(linear_gains, fund_frq, B, A)
%
% POLE_ZERO_PLACE  -- A GUI for filter design using pole-zero placement.
%
% Input Parameters:
%   linear_gains  -- desired transfer function gains at harmonically spaced points
%   fund_frq      -- fundamental frequency of the above gains
%   B             -- feedforward filter coefficients (see help filter) -- optional
%   A             -- feedback filter coefficients (see help filter) -- optional
%         NOTE: If B and A do not define filers with only conjugate pairs,
%               these inputs will be ignored.
%

% Bad programming practice warning: this code has not been properly commented.
%
% Some useful tidbits of info:
%   handles.Mode defines several "modes" of operation
%     Mode 0:  Normal operation
%     Mode 1:  Ready to place a zero
%     Mode 2:  Placed/selected zero or pole, waiting for update
%     Mode 4:  Ready to place a pole
%
% Also note that poles and zeros are treated identically in most places and 
%   are stored in the same arrays.  They can be distinguished using handles.pz
%   0 means zero, 1 means pole.


if nargin == 0  
    
    error('Not enough parameters');
    
elseif ~ischar(varargin{1})  % Create figure
 
    set(0,'ShowHiddenHandles','On')
    h = findobj('Tag','PZ_Place_Figure');
    if ~isempty(h)
        close(h);
    end
    set(0,'ShowHiddenHandles','Off');
    
    fig = openfig(mfilename,'reuse');
    
    % Use system color scheme for figure:
    set(fig,'Color',get(0,'defaultUicontrolBackgroundColor'));
    
    % Generate a structure of handles to pass to callbacks, and store it. 
    handles = guihandles(fig);
    
    handles.match_vals = varargin{1};
    if size(handles.match_vals,1) > size(handles.match_vals,2)
        handles.match_vals = handles.match_vals';
    end
    
    handles.match_frq = varargin{2};
    
    if length(varargin) > 2
        handles.B = varargin{3};
        if size(handles.B,1) > size(handles.B,2)
            handles.B = handles.B';
        end
        if length(varargin) > 3
            handles.A = varargin{4};
            if size(handles.A,1) > size(handles.A,2)
                handles.A = handles.A';
            end
        else
            handles.A = 1;
        end
    else
        handles.B = 1;
        handles.A = 1;
    end
    
    handles.fs = 8192;
    
    % Update handles structure
    guidata(fig, handles);
    
    if nargout > 0
        varargout{1} = fig;
    end
    
    set(fig,'HandleVisibility','On');
    
    % Perform a some initialization
    init(fig,[],handles);
    
    set(fig,'HandleVisibility','Callback');
    
    
else % INVOKE NAMED SUBFUNCTION OR CALLBACK
    
    try
        if (nargout)
            [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
        else
            feval(varargin{:}); % FEVAL switchyard
        end
    catch
        disp(lasterr);
    end
    
end


%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and 
%| sets objects' callback properties to call them through the FEVAL 
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the 
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.



% --------------------------------------------------------------------
function init(h, eventdata, handles, varargin)

draw_axes(h,[],handles);

% Store Data in handles

dif = length(handles.B) - length(handles.A);
[z,p,k] = tf2zp(handles.B,conv(handles.A,[1 zeros(1,dif)]));
handles.zero_list = reshape(conj([z;p]),[2,length([z;p])/2])';
handles.pz = reshape([zeros(size(z)); ones(size(p))],[2,length([z;p])/2])';

if any(handles.zero_list(:,1) ~= conj(handles.zero_list(:,2)))
    disp('Input filter is not composed of conjugate pairs.  Reinitializing...');
    handles.zero_list = zeros(0,1);
    handles.pz = zeros(0,2);
end

handles.zero_handle_list = [];
axes(handles.ZeroPlot);
for i = 1:size(handles.zero_list,1)
    hdl = line(real(handles.zero_list(i,1))*[1 1],imag(handles.zero_list(i,1))*[1 -1]);
    set(hdl,'LineStyle','none','LineWidth',1,'MarkerSize',8,'Color',[0 1 0]);
    set(hdl,'ButtonDownFcn',sprintf('pole_zero_place(''zero_select'',gcbo,[],guidata(gcbo));'));
    set(hdl,'EraseMode','xor');
    handles.zero_handle_list = [handles.zero_handle_list; hdl];
    if handles.pz(i,1) == 0
        set(hdl,'Marker','o');
    else
        set(hdl,'Marker','x');
    end
end
    
set(handles.Gain,'string',num2str(k));

handles.selected_zero = 0;

handles.Mode = 0;

guidata(h,handles);

slow_update_frq_resp(h,[],handles);


% --------------------------------------------------------------------
function slow_update_frq_resp(h, eventdata, handles, varargin);

axes(handles.MagResponse);

zlist = handles.zero_list(find(handles.pz == 0));
zlist = zlist(:);
set(handles.NumZeros,'string',num2str(size(zlist,1)));
if isempty(zlist);
    zlist = 0;
end
plist = handles.zero_list(find(handles.pz == 1));
plist = plist(:);
set(handles.NumPoles,'string',num2str(size(plist,1)));
if isempty(plist);
    plist = 0;
end
[handles.B,handles.A] = zp2tf(zlist,plist,str2double(get(handles.Gain,'string')));


[fr,w] = freqz(handles.B,handles.A,1024);
w = w/2/pi*handles.fs;
if get(handles.LinearButton,'value')
    handles.fr = plot(w,abs(fr),'b');
    hold on;
    stem(handles.match_frq*(1:length(handles.match_vals)),handles.match_vals,'r');
    hold off;
    ylabel('|H(\omega)|');
else
    handles.fr = plot(w,20*log10(abs(fr)),'b');        
    hold on;
    stem(handles.match_frq*(1:length(handles.match_vals)),20*log10(handles.match_vals),'r');
    hold off;
    ylabel('|H(\omega)| (dB)');
end
set(handles.fr,'EraseMode','xor');

xlabel('Frequency (Hz)');

axis tight; a1 = axis;
axis auto; a2 = axis;
axis([a1(1:2),a2(3:4)]);

if length(handles.match_vals) == 1
    H = freqz(handles.B,handles.A,[handles.match_frq 1]/handles.fs*2*pi);
    H = H(1);
else
    H = freqz(handles.B,handles.A,handles.match_frq*(1:length(handles.match_vals))*2*pi/handles.fs);
end
lin_err = sqrt(mean((handles.match_vals - abs(H)).^2));
db_err = sqrt(mean((20*log10(handles.match_vals) - 20*log10(abs(H))).^2));

set(handles.LinError,'string',num2str(lin_err));
set(handles.dBError,'string',num2str(db_err));

guidata(h,handles);

% --------------------------------------------------------------------

function draw_axes(h, eventdata, handles, varargin)

axes(handles.ZeroPlot);
ax_ext = 1.2;

ce = exp(j*2*pi*(0:.01:1));

h = line(real(ce),imag(ce));
set(h,'linestyle',':','color',[0 0 0])

h = line([0 0],[-ax_ext ax_ext]);
set(h,'linestyle',':','color',[0 0 0])
h = line([-ax_ext ax_ext],[0 0]);
set(h,'linestyle',':','color',[0 0 0])
axis([-ax_ext ax_ext -ax_ext ax_ext]);

xlabel('Real(z)');
ylabel('Imag(z)');

% --------------------------------------------------------------------
function varargout = ZeroPlot_CreateFcn(h, eventdata, handles, varargin)

% --------------------------------------------------------------------
function varargout = ZeroMag_Callback(h, eventdata, handles, varargin)

if handles.selected_zero > 0
    try
        newabs = eval(get(h,'string'));
        set(h,'string',num2str(newabs));
        zeros = handles.zero_list(handles.selected_zero,:);
        zeros = zeros./abs(zeros)*newabs;
        handles.zero_list(handles.selected_zero,:) = zeros;
        set(handles.zero_handle_list(handles.selected_zero),'XData',[real(zeros)]);
        set(handles.zero_handle_list(handles.selected_zero),'YData',[imag(zeros)]);
        
        slow_update_frq_resp(h,[],handles);
    catch
        set(h,'string',num2str(abs(handles.zero_list(handles.selected_zero,1))));
    end
else
    set(h,'string','');
end

% --------------------------------------------------------------------
function varargout = ZeroAngle_Callback(h, eventdata, handles, varargin)

if handles.selected_zero > 0
    try
        newang = eval(get(h,'string'));
        newang = mod(newang,2*pi);
        if newang > pi
            newang = newang - 2*pi;
        end
        set(h,'string',num2str(newang));
        zeros = handles.zero_list(handles.selected_zero,:);
        zeros = abs(zeros).*exp(j*[newang, -newang]);
        handles.zero_list(handles.selected_zero,:) = zeros;
        set(handles.zero_handle_list(handles.selected_zero),'XData',[real(zeros)]);
        set(handles.zero_handle_list(handles.selected_zero),'YData',[imag(zeros)]);
        
        slow_update_frq_resp(h,[],handles);
        
    catch
        set(h,'string',num2str(abs(handles.zero_list(handles.selected_zero,1))));
    end
else
    set(h,'string','');
end

% --------------------------------------------------------------------
function varargout = DeleteZero_Callback(h, eventdata, handles, varargin)

if handles.selected_zero > 0
    handles.zero_list(handles.selected_zero,:) = [];
    handles.pz(handles.selected_zero,:) = [];
    delete(handles.zero_handle_list(handles.selected_zero));
    handles.zero_handle_list(handles.selected_zero) = [];
    handles.selected_zero = 0;
    
    slow_update_frq_resp(h,[],handles);

end


% --------------------------------------------------------------------
function varargout = MagResponse_CreateFcn(h, eventdata, handles, varargin)


% --------------------------------------------------------------------
function varargout = AddZero_Callback(h, eventdata, handles, varargin)

if eventdata == 0
    if get(handles.AddZero,'value')
        handles.Mode = 1;
        set(handles.AddPole,'value',0)
    else
        handles.Mode = 0;
    end
else
    if get(handles.AddPole,'value')
        handles.Mode = 4;
        set(handles.AddZero,'value',0)
    else
        handles.Mode = 0;
    end
end

guidata(h,handles);


% --------------------------------------------------------------------
function varargout = ZeroPlot_ButtonDownFcn(h, eventdata, handles, varargin)

if handles.Mode == 1 | handles.Mode == 4
    pt = get(handles.ZeroPlot,'currentpoint');
    pt = pt(1,1:2);
    pt(2) = abs(pt(2));
    handles.zero_list = [handles.zero_list; pt(1) + j*pt(2), pt(1) - j*pt(2)];
    handles.selected_zero = size(handles.zero_list,1);
    
    axes(handles.ZeroPlot);
    hdl = line([pt(1) pt(1)],[pt(2) -pt(2)]);
    set(hdl,'LineStyle','none','LineWidth',2,'MarkerSize',8,'Color',[0 1 0]);
    set(hdl,'ButtonDownFcn',sprintf('pole_zero_place(''zero_select'',gcbo,[],guidata(gcbo));'));
    set(hdl,'EraseMode','xor');
    for ind = 1:length(handles.zero_handle_list)
        set(handles.zero_handle_list(ind),'LineWidth',1);
    end
    
    if handles.Mode == 1
        set(hdl,'Marker','o');
        handles.pz = [handles.pz; 0 0];
    else
        set(hdl,'Marker','x');
        handles.pz = [handles.pz; 1 1];
    end
    
    handles.Mode = 2;
    set(handles.AddZero,'value',0);
    set(handles.AddPole,'value',0);
    
    handles.zero_handle_list(end+1) = hdl; 
    
    set(handles.ZeroMag,'string',num2str(abs(pt(1) + i*pt(2))));
    set(handles.ZeroAngle,'string',num2str(angle(pt(1) + i*pt(2))));   
    
end

guidata(h,handles);


% --------------------------------------------------------------------
function varargout = zero_select(h, eventdata, handles, varargin)

if handles.Mode == 0
    
    for ind = 1:length(handles.zero_handle_list)
        set(handles.zero_handle_list(ind),'LineWidth',1);
    end
    
    handles.selected_zero = find(handles.zero_handle_list == h);
    set(h,'LineWidth',2);
    
    handles.Mode = 2;
    
    zero = handles.zero_list(handles.selected_zero,:);
    pt = get(handles.ZeroPlot,'currentpoint');
    if pt(1,2) > 0
        zero = zero(1);
    else
        zero = zero(2);
    end
    
    set(handles.ZeroMag,'string',num2str(abs(zero)));
    set(handles.ZeroAngle,'string',num2str(angle(zero)));   
    
end

guidata(h,handles);


% --------------------------------------------------------------------
function varargout = figure1_WindowButtonUpFcn(h, eventdata, handles, varargin)

if handles.Mode == 2
    handles.Mode = 0;
   
    zero = handles.zero_list(handles.selected_zero,:);
    pt = get(handles.ZeroPlot,'currentpoint');
    if pt(1,2) > 0
        zero = zero(1);
    else
        zero = zero(2);
    end
    
    set(handles.ZeroMag,'string',num2str(abs(zero)));
    set(handles.ZeroAngle,'string',num2str(angle(zero)));   
    
    slow_update_frq_resp(h,[],handles);
    
end



% --------------------------------------------------------------------
function varargout = Gain_Callback(h, eventdata, handles, varargin)

if isnan(str2double(get(h,'string')))
    set(h,'string','1.0');
end

handles.B = zp2tf(handles.zero_list(:),1,str2double(get(handles.Gain,'string')));
slow_update_frq_resp(h,[],handles);


% --------------------------------------------------------------------
function varargout = figure1_WindowButtonMotionFcn(h, eventdata, handles, varargin)

if handles.Mode == 2

    pt = get(handles.ZeroPlot,'currentpoint');

    zero = handles.zero_list(handles.selected_zero,:);
    if pt(1,2) > 0
        zero = zero(1);
    else
        zero = zero(2);
    end
    
    set(handles.ZeroMag,'string',num2str(abs(zero)));
    set(handles.ZeroAngle,'string',num2str(angle(zero)));   
    
    pt = pt(1,1:2);
    pt(2) = abs(pt(2));

    set(handles.zero_handle_list(handles.selected_zero),'XData',[pt(1) pt(1)]);
    set(handles.zero_handle_list(handles.selected_zero),'YData',[pt(2) -pt(2)]);

    handles.zero_list(handles.selected_zero,:) = [pt(1) + j*pt(2), pt(1) - j*pt(2)];
    
    zlist = handles.zero_list(find(handles.pz == 0));
    zlist = zlist(:);
    if isempty(zlist);
        zlist = 0;
    end
    plist = handles.zero_list(find(handles.pz == 1));
    plist = plist(:);
    if isempty(plist);
        plist = 0;
    end
    [handles.B,handles.A] = zp2tf(zlist,plist,str2double(get(handles.Gain,'string')));
    
    [fr,w] = freqz(handles.B,handles.A,1024);
    
    if get(handles.LinearButton,'value')
        set(handles.fr,'YData',abs(fr));
    else
        set(handles.fr,'YData',20*log10(abs(fr)));
    end        
    
    if length(handles.match_vals) == 1
        H = freqz(handles.B,handles.A,[handles.match_frq 1]/handles.fs*2*pi);
        H = H(1);
    else
        H = freqz(handles.B,handles.A,handles.match_frq*(1:length(handles.match_vals))*2*pi/handles.fs);
    end
    lin_err = sqrt(mean((handles.match_vals - abs(H)).^2));
    db_err = sqrt(mean((20*log10(handles.match_vals) - 20*log10(abs(H))).^2));
    
    set(handles.LinError,'string',num2str(lin_err));
    set(handles.dBError,'string',num2str(db_err));
    
  
    guidata(h,handles);
    
end

% --------------------------------------------------------------------
function varargout = LinearButton_Callback(h, eventdata, handles, varargin)

if get(h,'value') == 1
    set(handles.DecibelButton,'value',0);
else
    set(handles.DecibelButton,'value',1);
end

slow_update_frq_resp(h,[],handles);



% --------------------------------------------------------------------
function varargout = DecibelButton_Callback(h, eventdata, handles, varargin)

if get(h,'value') == 1
    set(handles.LinearButton,'value',0);
else
    set(handles.LinearButton,'value',1);
end

slow_update_frq_resp(h,[],handles);



% --------------------------------------------------------------------
function varargout = FIRCoefs_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = LinError_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = dBError_Callback(h, eventdata, handles, varargin)




% --------------------------------------------------------------------
function varargout = pushbutton2_Callback(h, eventdata, handles, varargin)

cstring = ['B = [ ' num2str(handles.B,'%7g, ')];
cstring(end-1:end) = ' ]';

disp(cstring);

cstring = ['A = [ ' num2str(handles.A,'%7g, ')];
cstring(end-1:end) = ' ]';

disp(cstring);

disp(' ');


% --------------------------------------------------------------------
function varargout = pushbutton3_Callback(h, eventdata, handles, varargin)


disp('Filter Coefficients:');

cstring = ['B = [ ' num2str(handles.B,'%7g, ')];
cstring(end-1:end) = ' ]';

disp(cstring);

cstring = ['A = [ ' num2str(handles.A,'%7g, ')];
cstring(end-1:end) = ' ]';

disp(cstring);

disp(' ');

delete(gcf);


% --------------------------------------------------------------------
function varargout = PlaySound_Callback(h, eventdata, handles, varargin)

g = glottal_source(handles.match_frq,2);

snd = filter(handles.B,handles.A,g);

soundsc(snd);


% --------------------------------------------------------------------
function varargout = printgui_Callback(h, eventdata, handles, varargin)

printdlg -crossplatform


% --------------------------------------------------------------------
function varargout = ToClipboard_Callback(h, eventdata, handles, varargin)

print -dbitmap

