Code covered by the BSD License  

Highlights from
GUI for Rational Fitting

image thumbnail

GUI for Rational Fitting

by

 

Use this GUI for fitting 2-port S-parameters using the functions from RF Toolbox

RationalFittingGUI(varargin)
%%Copyright 2013 The MathWorks, Inc
function varargout = RationalFittingGUI(varargin)
%RATIONALFITTINGGUI M-file for RationalFittingGUI.fig
%      RATIONALFITTINGGUI, by itself, creates a new RATIONALFITTINGGUI or raises the existing
%      singleton*.
%
%      H = RATIONALFITTINGGUI returns the handle to a new RATIONALFITTINGGUI or the handle to
%      the existing singleton*.
%
%      RATIONALFITTINGGUI('Property','Value',...) creates a new RATIONALFITTINGGUI using the
%      given property value pairs. Unrecognized properties are passed via
%      varargin to RationalFittingGUI_OpeningFcn.  This calling syntax produces a
%      warning when there is an existing singleton*.
%
%      RATIONALFITTINGGUI('CALLBACK') and RATIONALFITTINGGUI('CALLBACK',hObject,...) call the
%      local function named CALLBACK in RATIONALFITTINGGUI.M with the given input
%      arguments.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help RationalFittingGUI
% Last Modified by GUIDE v2.5 05-Mar-2013 11:28:48
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @RationalFittingGUI_OpeningFcn, ...
    'gui_OutputFcn',  @RationalFittingGUI_OutputFcn, ...
    'gui_LayoutFcn',  [], ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

function Rational_Fitting_GUI_Callback(hObject, eventdata, handles)
guidata(hObject,handles);

% Executed once when the GUi is made visible
function RationalFittingGUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
ch = get(handles.uipanel10,'Children');
for ii = 1:numel(ch),
    x = ch(ii);
    if(strcmp(get(x,'Tag'),'freq_resp'))
        handles.freq_resp = x;
    end
    if(strcmp(get(x,'Tag'),'ang_resp'))
        handles.ang_resp = x;
    end
end
% make the individual fitting button not visible
set(handles.SparamIndFit,'Visible','off');
guidata(hObject, handles);

% Execute when closing the GUI
function GUI_CloseRequestFcn(hObject, eventdata, handles)
delete(hObject);

function varargout = RationalFittingGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% Import new data in the GUI
function ImportData_Callback(hObject, eventdata, handles)
% When new data is imported, everything get reset.
% Section Import
set(handles.file_name, 'String', []);
set(handles.f1a, 'String', []);
set(handles.f2a, 'String', []);
set(handles.makePassive, 'Value', 0);
% Section plot
handles.plotParams = [0 0 0 0];
set(handles.radiobutton_s11, 'Value', 0);
set(handles.radiobutton_s12, 'Value', 0);
set(handles.radiobutton_s21, 'Value', 0);
set(handles.radiobutton_s22, 'Value', 0);
handles.PlotBw = 30;
set(handles.PlotBW, 'String', num2str(30));
% Section fit response
handles.FitInd = 0;
set(handles.FitIndividually, 'Value', 0);
set(handles.SparamIndFit, 'Visible', 'off');
handles.SIndFit = 'S11';
handles.IRF = [1, 1];
set(handles.SparamIndFit, 'Value', 1);
handles.tendsToZero = logical(0);
set(handles.tendstozero, 'Value', 0);
handles.MaxNpoles = 40;
set(handles.maxnpoles, 'String', num2str(40));
handles.Tol = -40;
set(handles.tolerance, 'String', num2str(-40));
handles.delay_factor = 0;
set(handles.del_fact, 'String', num2str(0));
handles.weightType = 1;
handles.weight_vec = [];
set(handles.weight_vec_menu, 'Value', 1);
weight_fig = findobj('Tag','WeightGUI');
if ~isempty(weight_fig)
    delete(weight_fig)
end
% Every time new data is imported, the data for the rational fitting is
% reset
temp = rfmodel.rational;
handles.rationalfunc = temp(ones(2,2));
handles.np = zeros(2,2);
handles.acc = nan(2,2);
handles.delayFit = 0;
set(handles.n_poles_h,'String', '');
set(handles.Acc,'String', '');
set(handles.del,'String', '');
% Section Export
handles.exportType = 1;
set(handles.SelectExport, 'Value', 1);
% Reset the plotting
set(gcf,'CurrentAxes',handles.freq_resp);
hold off;
plot(0,0);
axis([0 1 0 1]);
set(gcf,'CurrentAxes',handles.ang_resp);
hold off;
plot(0,0);
axis([0 1 0 1]);
% Save the reset data
guidata(hObject, handles);

% Get the file name
[FileName, PathName] = uigetfile({'*.s2p'},'2 Port S-parameter Touchstone file');
if FileName~=0
    [~, ~, ext] = fileparts(FileName);
    % make sure that the file is a 2 port (s2p)
    if strcmp(ext, '.s2p') || strcmp(ext, '.S2P')
        % Print the file name in the GUI
        file_name = strcat(PathName,FileName);
        set(handles.file_name, 'String', FileName);
        % Import the data from the file
        handles.FileName = file_name;                % name
        Sparam = sparameters(file_name);    % S-parameters
        freq = Sparam.Frequencies;   % frequency range
        f1 = Sparam.Frequencies(1,1)/1e6;            % F min
        f2 = Sparam.Frequencies(numel(freq),1)/1e6;  % Fmax
        z0 = Sparam.Impedance;       % Impedance: not used
        % Store the data in the structure handles
        handles.z0 = z0;
        handles.f1 = f1;
        handles.f2 = f2;
        handles.freq = freq;
        handles.freqPlot = freq;        % Default initialization for plotting
        handles.Sparam = Sparam;
        handles.SparamOrig = Sparam; % This is used to turn on and off "make passive"
        handles.resp = zeros(2,2,0);    % Default initialization
        % Print the frequency range of the data in the GUI
        set(handles.f1a, 'String', handles.f1);
        set(handles.f2a, 'String', handles.f2);
        handles.weight_vec = ones(1,length(freq));
        guidata(hObject, handles);
        % This is required to update the custom weights in case new data is
        % read in, keeping non-default settings for the weights.
        weight_vec_menu_Callback(handles.weight_vec_menu, eventdata, handles);
    else
        errordlg('File must be .s2p - 2 ports ONLY','Error');
    end
end
% Make the data passive (or return to the original data)
function makePassive_Callback(hObject, eventdata, handles)
if isfield(handles, 'Sparam')
    if get(hObject, 'Value')==1
        Sparam = handles.Sparam;
        handles.SparamOrig = Sparam;
        Sparam = makepassive(Sparam);
        handles.Sparam = Sparam;
    else
        handles.Sparam = handles.SparamOrig;
    end
    guidata(hObject,handles);
end

%% Four radio buttons used for plotting
function radiobutton_s11_Callback(hObject, eventdata, handles)
if get(hObject,'Value')
    handles.plotParams(1) = 1;
else
    handles.plotParams(1) = 0;
end
guidata(hObject,handles);

function radiobutton_s21_Callback(hObject, eventdata, handles)
if get(hObject,'Value')
    handles.plotParams(3) = 1;
else
    handles.plotParams(3) = 0;
end
guidata(hObject,handles);

function radiobutton_s12_Callback(hObject, eventdata, handles)
if get(hObject,'Value')
    handles.plotParams(2) = 1;
else
    handles.plotParams(2) = 0;
end
guidata(hObject,handles);

function radiobutton_s22_Callback(hObject, eventdata, handles)
if get(hObject,'Value')
    handles.plotParams(4) = 1;
else
    handles.plotParams(4) = 0;
end
guidata(hObject,handles);

% Read the plot excess bandwidth
function PlotBW_Callback(hObject, eventdata, handles)
bw = str2double(get(hObject,'String'));
% Error checking
if isnan(bw) || bw <0
    handles.PlotBw = 0;
    errordlg('Plot excess bandwidth must be a positive number','Error');
else
    handles.PlotBw = bw;
end
guidata(hObject, handles);

% Plotting callback
function Plot_Callback(hObject, eventdata, handles)
Sparam = handles.Sparam;
plotParams = handles.plotParams;
freq = handles.freq;

% Plot over excess bandwidth
if handles.PlotBw ~= 0
    step = handles.freq(2) - handles.freq(1);
    range = handles.freq(end) - handles.freq(1);
    minF = handles.freq(1) - range/100*handles.PlotBw;
    if minF < 0
        % make sure that frequency is only positive
        minF =  0;
    end
    maxF = handles.freq(end) + range/100*handles.PlotBw;
    freqPlot1 = fliplr((handles.freq(1)-step:-step:minF));
    freqPlot2 = (handles.freq(end)+step:step:maxF);
    freqPlot = [freqPlot1'; freq; freqPlot2'];
else
    minF = handles.freq(1);
    maxF = handles.freq(end);
    freqPlot=handles.freq;
end
handles.freqPlot = freqPlot;

% If present, determine the frequency response of the fitting
% This is done everytime you plot, in case the plotting excess BW changes
% or in case the data changes
if isempty(handles.resp)
    fit = 0;
else
    fit = 1;
    rationalfunc = handles.rationalfunc;
    for i = 1:2
        for j = 1:2
            resp(i,j,:) = freqresp(rationalfunc(i,j), freqPlot);
        end
    end
    handles.resp = resp;
    np = handles.np;
end

% Create and rfckt used for plotting the data
if ispassive(Sparam.Parameters)
    Data = rfckt.passive('NetworkData', ...
        rfdata.network('Data',  Sparam.Parameters, 'Freq', freq));
else
    warning('off','rf:rfdata:data:checkproperty:PoutNotConsistentWithNetworkData')
    Data = rfckt.amplifier('NetworkData', ...
        rfdata.network('Data',  Sparam.Parameters, 'Freq', freq));
end
% Plotting without fitting data
if fit == 0
    for i = (1:4)
        switch i
            case 1
                str = 'S11';
                c = 'b';
            case 2
                str = 'S12';
                c = 'r';
            case 3
                str = 'S21';
                c = 'g';
            case 4
                str = 'S22';
                c = 'k';
        end
        % Plot only selected data
        if plotParams(i) == 1
            % Plot magnitude response
            set(gcf,'CurrentAxes',handles.freq_resp);
            hfreq_resp = plot(Data, str);
            set(hfreq_resp, 'Color', c);
            set(hfreq_resp, 'LineWidth', 2);
            tmp = get(hfreq_resp, 'UserData');
            if strcmp(tmp.XUnit,'[MHz]')
                f = 1e6;
            elseif strcmp(tmp.XUnit,'[GHz]')
                f = 1e9;
            elseif strcmp(tmp.XUnit,'[kHz]')
                f = 1e3;
            end
            title(sprintf('Magnitude response'),'fonts',10);
            set(gca, 'XLimMode', 'manual', 'XLim', [minF maxF]./f);
            hold on;
            % Plot phase response
            set(gcf,'CurrentAxes',handles.ang_resp);
            hang_resp = plot(Data, str, 'angle');
            set(hang_resp, 'Color', c);
            set(hang_resp, 'LineWidth', 2);
            title(sprintf('Phase response'),'fonts',10);
            set(gca, 'XLimMode', 'manual', 'XLim', [minF maxF]./f);
            hold on;
        end
    end
else % Plot including fitted data
    for i = (1:2)
        for j = (1:2)
            k = 2*(i-1)+j;
            switch k
                case 1
                    str = 'S11';
                    c1 = 'b';
                    c2 = 'b:';
                case 2
                    str = 'S12';
                    c1 = 'r';
                    c2 = 'r:';
                case 3
                    str = 'S21';
                    c1 = 'g';
                    c2 = 'g:';
                case 4
                    str = 'S22';
                    c1 = 'k';
                    c2 = 'k:';
            end
            % plot only selected data
            if plotParams(k) == 1
                plotOnlyOne = sum(plotParams);
                % Plot magnite response and fitting
                set(gcf, 'CurrentAxes',handles.freq_resp);
                S_param = Data.AnalyzedResult.S_Parameters;
                freqPlot = handles.freqPlot;
                resp = handles.resp;
                npoles = handles.npoles;
                hfreq_resp = plot(freq*1.e-9, 20*log10(abs(squeeze(S_param(i,j,:)))), c1, ...
                    freqPlot*1.e-9, 20*log10(abs(squeeze(resp(i,j,:)))), c2);
                set(hfreq_resp, 'LineWidth', 2);
                if plotOnlyOne == 1
                    title([str,' Rational Fitting with ',num2str(np(i,j)), ' poles'],'fonts',10);
                else
                    title(['Rational Fitting with ',num2str(npoles), ' poles max'],'fonts',10);
                end
                ylabel('Magnitude (decibels)'); xlabel('Frequency (GHz)');
                % Update the legend
                l = legend;
                if isempty(l)
                    legend([str ' Original data'], [str ' Fitting result']);
                else
                    string = get(l, 'String');
                    exist_legend = numel(string);
                    string{exist_legend+1} = [str ' Original data'];
                    string{exist_legend+2} = [str ' Fitting result'];
                    legend(string{1:end}, 'Location', 'NorthEast');
                end
                grid on;
                hold on;
                
                % Plot phase response
                set(gcf,'CurrentAxes',handles.ang_resp);
                hang_resp = plot(freq*1.e-9, unwrap(angle(squeeze(S_param(i,j,:))))*180/pi, c1, ...
                    freqPlot*1.e-9, unwrap(angle(squeeze(resp(i,j,:))))*180/pi, c2);
                set(hang_resp, 'LineWidth', 2);
                if plotOnlyOne == 1
                    title([str,' Rational Fitting with ',num2str(np(i,j)), ' poles'],'fonts',10);
                else
                    title(['Rational Fitting with ',num2str(npoles), ' poles max'],'fonts',10);
                end
                ylabel('Absolute phase (deg)'); xlabel('Frequency (GHz)');
                % Update legend
                l = legend;
                if isempty(l)
                    legend([str ' Original data'], [str ' Fitting result']);
                else
                    string = get(l, 'String');
                    exist_legend = numel(string);
                    string{exist_legend+1} = [str ' Original data'];
                    string{exist_legend+2} = [str ' Fitting result'];
                    legend(string{1:end}, 'Location', 'NorthEast');
                end
                grid on;
                hold on;
            end
        end
    end
end

if sum(plotParams) == 0
    set(gcf,'CurrentAxes',handles.freq_resp);
    hold off;
    plot(0,0);
    axis([0 1 0 1]);
    set(gcf,'CurrentAxes',handles.ang_resp);
    hold off;
    plot(0,0);
    axis([0 1 0 1]);
end

set(gcf,'CurrentAxes',handles.freq_resp);
hold off;
set(gcf,'CurrentAxes',handles.ang_resp);
hold off;
guidata(hObject, handles);

% Fit individually checkbox
function FitIndividually_Callback(hObject, eventdata, handles)
FitInd = get(hObject,'Value');
% Enable / disble the pulldown menu to select which S-parameter to fit
if FitInd == 1
    set(handles.SparamIndFit,'Visible','on');
else
    set(handles.SparamIndFit,'Visible','off');
end
handles.FitInd = FitInd;
guidata(hObject,handles);
% This is required in case you toggle on and off the individual fitting to
% reset the correct values for the weights
weight_vec_menu_Callback(handles.weight_vec_menu, eventdata, handles);

% Reads which S-parameter should be fitted individually
function SparamIndFit_Callback(hObject, eventdata, handles)
SIndFit = cellstr(get(hObject,'String'));
handles.SIndFit = SIndFit{get(hObject,'Value')};
switch handles.SIndFit
    case 'S11'
        i = 1; j = 1;
    case 'S12'
        i = 1; j = 2;
    case 'S21'
        i = 2; j = 1;
    case 'S22'
        i = 2; j = 2;
end
handles.IRF = [i, j];
guidata(hObject,handles);
weight_vec_menu_Callback(handles.weight_vec_menu, eventdata, handles);

% Tends to zero checkbox
function tendstozero_Callback(hObject, eventdata, handles)
tendsToZero = get(hObject,'Value');
handles.tendsToZero = logical(tendsToZero);
guidata(hObject,handles);

% Read the maximum number of poles used for the fitting
function maxnpoles_Callback(hObject, eventdata, handles)
MaxNpoles = str2double(get(hObject,'String'));
% Error checking
if isnan(MaxNpoles) || MaxNpoles<=0
    handles.MaxNpoles = 40;
    set(hObject, 'String', num2str(40));
    errordlg('Max number of poles must be a positive number','Error');
else
    handles.MaxNpoles = MaxNpoles;
end
guidata(hObject,handles);

% Read the tolerance for the fitting
function tolerance_Callback(hObject, eventdata, handles)
Tol = str2double(get(hObject, 'String'));
% Error checking
if isnan(Tol) || Tol>0
    set(hObject, 'String', '-40');
    errordlg('Tolerance must be a negative number','Error');
    Tol = -40;
end
% Save the new Rational fit relative tolerance
handles.Tol = Tol;
guidata(hObject,handles)

% Read the delay factor for the fitting
function del_fact_Callback(hObject, eventdata, handles)
delay_factor = str2double(get(hObject, 'String'));
% Error checking
if isnan(delay_factor) || (delay_factor<0) || (delay_factor>1)
    set(hObject, 'String', '0');
    errordlg('Delay factor must be a number between 0 and 1','Error');
    delay_factor = 0;
end
% Save the new delay factor
handles.delay_factor = delay_factor;
guidata(hObject,handles)

% Nightmare callback to determine and update weights
function weight_vec_menu_Callback(hObject, eventdata, handles)
% If the object does not exists then the weights cannot be computed
if isfield(handles, 'Sparam')
    Data = handles.Sparam.Parameters;
    freq = handles.freq;
    weightType  = get(hObject,'value');
    % Detrmine the Weight vector. Is is 4x4xNumFreqs is fit individually is 
    % used, otherwise it is just a vecotr with as many elements as the
    % number of frequencies
    switch(weightType)
        case 1 % Uniform weights
            if handles.FitInd
                weight_vec = ones(2,2,length(freq));
            else
                weight_vec = ones(1,length(freq));
            end
        case 2 % Non uniform weights work only for individual fitting
            if handles.FitInd
                % Weight vector initialization if the vector exists for
                % non-individual fitting
                weight_vec = handles.weight_vec;
                if ~((size(weight_vec,1) == 2 ) && size(weight_vec,2) == 2 )
                    weight_vec = ones(2,2,length(freq));
                end
                d = squeeze(abs(Data(handles.IRF(1),handles.IRF(2),:)));
                M = repmat(max(d),[length(freq) 1]);
                weight_vec(handles.IRF(1),handles.IRF(2),:) = squeeze(ones(size(d))./abs(d+eps).*M);
            else
                errordlg('Non-uniform weights only for individual fitting','Error');
                set(handles.weight_vec_menu,'Value',1);
                weight_vec = ones(1,length(freq));
            end
        case 3 % Custom weights work only for individual fitting
            if handles.FitInd
                weight_vec = handles.weight_vec;
                % Weight vector initialization if the vector exists for
                % non-individual fitting
                if ~((size(weight_vec,1) == 2 ) && size(weight_vec,2) == 2 )
                    weight_vec = ones(2,2,length(freq));
                end
                weight_fig = findobj('Tag','WeightGUI');
                if isempty(weight_fig) % New weight GUI
                    simple_weight([], freq, squeeze(Data(handles.IRF(1),handles.IRF(2),:)));
                else % reconnect to the existing GUI
                    set(0, 'CurrentFigure', weight_fig);
                    set(weight_fig, 'Selected', 'on');
                    simple_weight('set_FandSxy', weight_fig, freq, squeeze(Data(handles.IRF(1),handles.IRF(2),:)));
                end
                % This is needed when the GUI is first opened if weight
                % is not get written in the workspace fast enough
                weight = simple_weight('get_weight');
                weight_vec(handles.IRF(1),handles.IRF(2),:) = weight;
            else % error out when not using it for individual fitting
                % Gray out the custom option?
                errordlg('Custom weights only for individual fitting','Error');
                set(handles.weight_vec_menu,'Value',1);
                weight_vec = ones(1,length(freq));
            end
    end
else
    errordlg('S-parameters must be loaded to compute weights','Error');
end
h_fig=findobj('Tag','GUI');
handles=guidata(h_fig);
handles.weightType = weightType;
handles.weight_vec = weight_vec;
guidata(hObject,handles)

% Rational fitting callback
function calc_rat_fit_Callback(hObject, eventdata, handles)
freq = handles.freq;
Data = handles.Sparam.Parameters;
warning('off','rf:rationalfit:ErrorToleranceNotMet');
if handles.FitInd == 0
    % Fit individually
    % Make sure that you get the right data, from the right GUI before
    % fitting
    h_fig = findobj('Tag','GUI');
    handles = guidata(h_fig);
    Tol = handles.Tol;
    MaxNpoles = handles.MaxNpoles;
    weight_vec = handles.weight_vec;
    freqPlot = handles.freqPlot;
    tendsToZero = handles.tendsToZero;
    
    % Rational fitting for all 4 S-parameters
    [rationalfunc Acc]= rationalfit(freq, Data, ...
        'Tolerance', Tol, ...
        'Weight', weight_vec, ...
        'DelayFactor', handles.delay_factor, ...
        'NPoles',  [ 1 handles.MaxNpoles], ...
        'TendsToZero', tendsToZero, ...
        'WaitBar', true);
    
    % Make the rational fitting GUI active
    set(h_fig, 'Selected', 'on');
    Acc = ceil(Acc);
    acc = Acc.*ones(2,2);
    npoles = length(rationalfunc(1,1).A);
    np = npoles*ones(2,2);
    resp = zeros(2,2,length(freqPlot));
    for i = 1:2
        for j = 1:2
            resp(i,j,:) = freqresp(rationalfunc(i,j), freqPlot);
            np(i,j) = npoles;
            %% Plot results rational fit
            k = 2*(i-1)+j;
            switch k
                case 1
                    str = 'S11';
                    c1 = 'b';
                    c2 = 'b:';
                case 2
                    str = 'S12';
                    c1 = 'r';
                    c2 = 'r:';
                case 3
                    str = 'S21';
                    c1 = 'g';
                    c2 = 'g:';
                case 4
                    str = 'S22';
                    c1 = 'k';
                    c2 = 'k:';
            end
            if handles.plotParams(k) == 1
                % Plot magniture response
                set(0,'CurrentFigure',handles.GUI)
                set(gcf, 'CurrentAxes',handles.freq_resp);
                hfreq_resp = plot(freq*1.e-9, 20*log10(abs(squeeze(Data(i,j,:)))), c1, ...
                    freqPlot*1.e-9, 20*log10(abs(squeeze(resp(i,j,:)))), c2);
                set(hfreq_resp, 'LineWidth', 2);
                title(['Rational Fitting with ',num2str(np(i,j)), ' poles shared by all S-parameters'],'fonts',10);
                ylabel('Magnitude (decibels)'); xlabel('Frequency (GHz)');
                l = legend;
                if isempty(l)
                    legend([str ' Original data'], [str ' Fitting result']);
                else
                    string = get(l, 'String');
                    exist_legend = numel(string);
                    string{exist_legend+1} = [str ' Original data'];
                    string{exist_legend+2} = [str ' Fitting result'];
                    legend(string{1:end}, 'Location', 'NorthEast');
                end
                grid on;
                hold on;
                % Plot phase response
                set(gcf,'CurrentAxes',handles.ang_resp);
                hang_resp = plot(freq*1.e-9, unwrap(angle(squeeze(Data(i,j,:))))*180/pi, c1, ...
                    freqPlot*1.e-9, unwrap(angle(squeeze(resp(i,j,:))))*180/pi, c2);
                set(hang_resp, 'LineWidth', 2);
                title(['Rational Fitting with ',num2str(np(i,j)), ' poles shared by all S-parameters'],'fonts',10);
                ylabel('Absolute phase (deg)'); xlabel('Frequency (GHz)');
                l = legend;
                if isempty(l)
                    legend([str ' Original data'], [str ' Fitting result']);
                else
                    string = get(l, 'String');
                    exist_legend = numel(string);
                    string{exist_legend+1} = [str ' Original data'];
                    string{exist_legend+2} = [str ' Fitting result'];
                    legend(string{1:end}, 'Location', 'NorthEast');
                end
                grid on;
                hold on;
            end            
        end
    end
else
    % Individual fitting
    % Which S-parameter to fit
    SIndFit = handles.SIndFit;
    switch SIndFit
        case 'S11'
            i = 1; j = 1;
        case 'S12'
            i = 1; j = 2;
        case 'S21'
            i = 2; j = 1;
        case 'S22'
            i = 2; j = 2;
    end
    
    % Initialize the response of the fitting (if it does not exist)
    resp = handles.resp;
    if isempty(resp)
        resp = zeros(2,2,length(handles.freqPlot));
    end
    acc = handles.acc;
    np = handles.np;
    rationalfunc = handles.rationalfunc;
    
    handles.IRF = [i, j];
    guidata(hObject, handles);
    
    % Make sure that you get the right data, from the right GUI before
    % fitting
    h_fig = findobj('Tag','GUI');
    handles = guidata(h_fig);
    Tol = handles.Tol;
    MaxNpoles = handles.MaxNpoles;
    weight_vec = handles.weight_vec;
    freqPlot = handles.freqPlot;
    tendsToZero = handles.tendsToZero;
    delay_factor = handles.delay_factor;
    if (handles.weightType == 3)
        weight_fig = findobj('Tag','WeightGUI');
        if isempty(weight_fig) % New weight GUI
            simple_weight([], freq, squeeze(Data(handles.IRF(1),handles.IRF(2),:)));
        end
        weight = simple_weight('get_weight');
        weight_vec(i,j,:) = weight;
    end
    
    % rational fitting
    [rationalfunc(i,j) acc(i,j)]= rationalfit(freq, Data(i,j,:), ...
        'Tolerance', Tol, ...
        'Weight', squeeze(weight_vec(i,j,:)), ...
        'DelayFactor', delay_factor, ...
        'NPoles',  [ 1 handles.MaxNpoles], ...
        'TendsToZero', tendsToZero, ...
        'WaitBar',true);
    
    np(i,j) = length(rationalfunc(i,j).A);
    resp(i,j,:) = freqresp(rationalfunc(i,j), freqPlot);
    if acc(i,j) > Tol
        h = msgbox(['Accuracy achieved = ', num2str(acc(i,j))],'Warning');
        uiwait(h);
    end
    set(0,'CurrentFigure',handles.GUI);
    set(h_fig, 'Selected', 'on');
    %% Plot results rational fit
    k = 2*(i-1)+j;
    switch k
        case 1
            str = 'S11';
            c1 = 'b';
            c2 = 'b:';
        case 2
            str = 'S12';
            c1 = 'r';
            c2 = 'r:';
        case 3
            str = 'S21';
            c1 = 'g';
            c2 = 'g:';
        case 4
            str = 'S22';
            c1 = 'k';
            c2 = 'k:';
    end
    % Plot magnitude response
    set(gcf, 'CurrentAxes',handles.freq_resp);
    hfreq_resp = plot(freq*1.e-9, 20*log10(abs(squeeze(Data(i,j,:)))), c1, ...
        freqPlot*1.e-9, 20*log10(abs(squeeze(resp(i,j,:)))), c2);
    set(hfreq_resp, 'LineWidth', 2);
    title([str,' Rational Fitting with ',num2str(np(i,j)), ' poles'],'fonts',10);
    ylabel('Magnitude (decibels)'); xlabel('Frequency (GHz)');
    l = legend;
    if isempty(l)
        legend([str ' Original data'], [str ' Fitting result']);
    else
        string = get(l, 'String');
        exist_legend = numel(string);
        string{exist_legend+1} = [str ' Original data'];
        string{exist_legend+2} = [str ' Fitting result'];
        legend(string{1:end}, 'Location', 'NorthEast');
    end
    grid on;
    
    % plot phase response
    set(gcf,'CurrentAxes',handles.ang_resp);
    hang_resp = plot(freq*1.e-9, unwrap(angle(squeeze(Data(i,j,:))))*180/pi, c1, ...
        freqPlot*1.e-9, unwrap(angle(squeeze(resp(i,j,:))))*180/pi, c2);
    set(hang_resp, 'LineWidth', 2);
    title([str, ' Rational Fitting with ',num2str(np(i,j)), ' poles'],'fonts',10);
    ylabel('Absolute phase (deg)'); xlabel('Frequency (GHz)');
    l = legend;
    if isempty(l)
        legend([str ' Original data'], [str ' Fitting result']);
    else
        string = get(l, 'String');
        exist_legend = numel(string);
        string{exist_legend+1} = [str ' Original data'];
        string{exist_legend+2} = [str ' Fitting result'];
        legend(string{1:end}, 'Location', 'NorthEast');
    end
    grid on;
    npoles = max(max(np));
    Acc = ceil(max(max(acc)));
end

handles.resp = resp;
handles.rationalfunc = rationalfunc;
handles.npoles = npoles;
handles.np = np;
handles.acc = acc;
handles.delayFit = rationalfunc(2,1).Delay;
set(handles.n_poles_h,'String', npoles);
set(handles.Acc,'String', Acc);
set(handles.del,'String', handles.delayFit./1e-9);

set(gcf,'CurrentAxes',handles.freq_resp);
hold off;
set(gcf,'CurrentAxes',handles.ang_resp);
hold off;
guidata(hObject, handles);

% Get the type of export
function SelectExport_Callback(hObject, eventdata, handles)
handles.exportType = get(hObject,'value');
guidata(hObject,handles)

% Export the fitting
function Export_Callback(hObject, eventdata, handles)
rationalfunc = handles.rationalfunc;
num_lines = 1;
% make sure that all 4 fittings are present. It might happen that some is
% missing after using "individual fitting"
if isempty(rationalfunc(1,1).A) || isempty(rationalfunc(1,2).A) || ...
        isempty(rationalfunc(2,1).A) || isempty(rationalfunc(2,2).A)
    h = msgbox('At least one fitting is missing','Warning');
    uiwait(h);
end

% Three possible exports
switch(handles.exportType)
    case 1 % Export to MATLAB Workspace rational object
        prompt = {'Enter name for rational model:'};
        dlg_title = 'Save to MATLAB workspace';
        num_lines = 1;
        def = {'rationalmodel'};
        ratmodel = inputdlg(prompt, dlg_title, num_lines, def);
        try 
            assignin('base', ratmodel{1}, rationalfunc);
        catch
            errordlg('Error: Export failed');
        end
        guidata(hObject, handles)
    case 2 % Export to MATLAB Workspace A C D and delay
        A = get(handles.rationalfunc,'A');
        C = get(handles.rationalfunc,'C');
        D = get(handles.rationalfunc,'D');
        try 
            assignin('base', 'A', A);
            assignin('base', 'C', C);
            assignin('base', 'D', D);
            assignin('base', 'Delay', handles.rationalfunc(2,1).Delay);
        catch
            errodlg('Error: Export failed');
        end
    case 3 % Export to .mat File
        prompt = {'Enter mat output file name with extension (*.mat):'};
        dlg_title = 'Generate Rational fit data file';
        rationalfunc = handles.rationalfunc;
        def = {'rationalmodel.mat'};
        outfile = inputdlg(prompt, dlg_title, num_lines, def);
        outfile1 = char(outfile);
        try 
            save(outfile1, 'rationalfunc');
        catch
            errordlg('Error: Export failed');
        end
        guidata(hObject,handles)
end
guidata(hObject,handles)

%% Create Fcn - nothing to say about them
function GUI_CreateFcn(hObject, eventdata, handles)

function SparamIndFit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function PlotBW_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function SelectExport_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function weight_vec_menu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function maxnpoles_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function tolerance_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), ...
        get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function del_fact_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

Contact us