www.gusucode.com > 利用MATLAB GUI设计滤波器界面,可以设计IIR滤波器 > AFD/GuiBuildStageInformation.m

    function varargout = GuiBuildStageInformation(varargin)
% GUIBUILDSTAGEINFORMATION M-file for guiBuildStageInformation.fig
%      GUIBUILDSTAGEINFORMATION, by itself, creates a new GUIBUILDSTAGEINFORMATION or raises the existing
%      singleton*.
%
%      H = GUIBUILDSTAGEINFORMATION returns the handle to a new GUIBUILDSTAGEINFORMATION or the handle to
%      the existing singleton*.
%
%      GUIBUILDSTAGEINFORMATION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUIBUILDSTAGEINFORMATION.M with the given input arguments.
%
%      GUIBUILDSTAGEINFORMATION('Property','Value',...) creates a new GUIBUILDSTAGEINFORMATION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before guiBuildStageInformation_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to guiBuildStageInformation_OpeningFcn via varargin.
%
%      *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 guiBuildStageInformation

% Last Modified by GUIDE v2.5 27-Dec-2003 17:45:11

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GuiBuildStageInformation_OpeningFcn, ...
                   'gui_OutputFcn',  @GuiBuildStageInformation_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(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


% --- Executes just before guiBuildStageInformation is made visible.
function GuiBuildStageInformation_OpeningFcn(hObject, eventdata, handles, varargin)
% Choose default command line output for guiBuildStageInformation
handles.output = hObject;

% Load data
if length(varargin) == 2
    strStage = varargin{1};
    curStage = varargin{2};
else
   temp=load('curStage.mat');
   disp([mfilename ' called in debug mode using StageInformation.mat datafile'])
   strStage = temp.curStage;
   curStage = 1;
end
handles.curStage = curStage;
handles.strStage = strStage;

% create the figure title and status bar
set(handles.uiBuildStageInformation,'Name',sprintf('Stage %g',curStage))
set(handles.uitxTitle,'String',['Detailed Information: Stage ' num2str(curStage)])
set(handles.uitxStatus,'String','Select the information to view')

% Figure out if there is a toleranced R or C and if so display it by default
global strCircuit
if strCircuit.nRTol + strCircuit.nCTol == 0
    set(handles.uirbComponentsExact,'Value',1)
    set(handles.uirbComponentsToleranced,'Value',0)
else
    set(handles.uirbComponentsExact,'Value',0)
    set(handles.uirbComponentsToleranced,'Value',1)
end
% create upper plot and table, and lower plot and table
RecalculateRightSide(handles)

% fill the component selection popup
nCSelect = length(strStage.vfCSelect);
nRSelect = length(strStage.vfRSelect);
nCCalc = length(strStage.vfCCalc);
nRCalc = length(strStage.vfRCalc);
cComponents = {};
for i=1:nCSelect
    cComponents{end+1} = ['C' char(96+i)];
end
for i=1:nRSelect
    cComponents{end+1} = ['R' char(96+i)];
end
for i=1:nCCalc
    cComponents{end+1} = ['C' num2str(i)];
end
for i=1:nRCalc
    cComponents{end+1} = ['R' num2str(i)];
end
set(handles.uipmSensitivityIn,'String',cComponents)
CalculateSensitivity(handles)

% Update handles structure
guidata(hObject, handles);

% --- Outputs from this function are returned to the command line.
function varargout = GuiBuildStageInformation_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Components, Frequency, Gain                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uirbComponentsExact_Callback(hObject, eventdata, handles)
set(hObject,'Value',1);
set(handles.uirbComponentsToleranced,'Value',0)
RecalculateRightSide(handles)

function uirbComponentsToleranced_Callback(hObject, eventdata, handles)
set(hObject,'Value',1);
set(handles.uirbComponentsExact,'Value',0)
RecalculateRightSide(handles)

function uirbFrequencyRads_Callback(hObject, eventdata, handles)
set(hObject,'Value',1);
set(handles.uirbFrequencyHz,'Value',0)
RecalculateRightSide(handles)

function uirbFrequencyHz_Callback(hObject, eventdata, handles)
set(hObject,'Value',1);
set(handles.uirbFrequencyRads,'Value',0)
RecalculateRightSide(handles)

function uirbGainLinear_Callback(hObject, eventdata, handles)
set(hObject,'Value',1);
set(handles.uirbGainDb,'Value',0)
RecalculateRightSide(handles)

function uirbGainDb_Callback(hObject, eventdata, handles)
set(hObject,'Value',1);
set(handles.uirbGainLinear,'Value',0)
RecalculateRightSide(handles)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              CustomFrequency                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uitxCustomFreq_Callback(hObject, eventdata, handles)
RecalculateRightSide(handles)

function uipmCustomFreq_Callback(hObject, eventdata, handles)
RecalculateRightSide(handles)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Sensitivity                                 %
%          uipmSensitivityIn, uipmSensitivityOut, uitxSensitivity           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uipmSensitivityIn_Callback(hObject, eventdata, handles)
CalculateSensitivity(handles)

function uipmSensitivityOut_Callback(hObject, eventdata, handles)
CalculateSensitivity(handles)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           RecalculateRightSide                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function RecalculateRightSide(handles)
% setup
strStage = handles.strStage;
nStage = handles.curStage;
if get(handles.uirbComponentsExact,'Value') % exact components
    % Temporarily, refigure entire circuit out as if RTol=CTol=0 and use the 
    % calculated z1,p1,k1 from this for the exact stage.
    % The issue is that the ideal stage z,p,k is not a function of the circuit type,
    % so, e.g., a MFB biquad section that cannot select a gain may have a RTol=CTol=0
    % yielding a max gain of 0.1, whereas an "ideal" biquad that lets one set the
    % gain may have a gain of 1.  We don't want the user getting confused by seeing
    % a 0 toleranced implementation with a gain of 0.1 and an ideal stage having a
    % gain of 1, so, we do this.
    global strCircuit
    oldstrCircuit=strCircuit;
    strCircuit.RTol=0;
    strCircuit.CTol=0;
    for i=1:strCircuit.nBiquads
        strCircuit.vStage(i) = BuildCircuit_UpdateComponents(strCircuit.vStage(i),0,0);
    end
    strCircuit = BuildCircuit_LastStage(strCircuit);
    z = strCircuit.vStage(nStage).z1; % what they would be with 0 toleranced components
    p = strCircuit.vStage(nStage).p1;
    k = strCircuit.vStage(nStage).k1;
    % now account for rounding by making small values = 0
    z = RemoveSmallValues(z);
    p = RemoveSmallValues(p);
    strCircuit = oldstrCircuit;
else
    z = strStage.z1; zIdeal=strStage.z;
    p = strStage.p1; pIdeal=strStage.p;
    k = strStage.k1; kIdeal=strStage.k;
end
isHz = get(handles.uirbFrequencyHz,'Value');
isDb = get(handles.uirbGainDb,'Value');

% begin by loading the schematic
persistent strGraphic
if isempty(strGraphic)
    t=load('StageInformation.mat');
    strGraphic = t.strGraphic;
end
nPoles=length(strStage.p);
nZeros=length(strStage.z);
switch nPoles
    case 0
        schName = 'ZO';
    case 1
        if nZeros == 1
            schName = 'FO_HP';
        else
            schName = 'FO_LP';
        end
    case 2
        if nZeros == 0
            schName = 'SO_LP';
        elseif strStage.z(1) == 0
            schName = 'SO_HP';
        else
            schName = 'SO_Notch';
        end
    otherwise
        error('Unknown strStage type')
end
axes(handles.uiaxEqn)
hImage = image(strGraphic.(schName), 'Parent', handles.uiaxEqn);
set(handles.uiaxEqn, ...
    'Visible', 'off', ...
    'YDir'   , 'reverse'       , ...
    'XLim'   , get(hImage,'XData'), ...
    'YLim'   , get(hImage,'YData')  ...
    );
set(handles.uiBuildStageInformation,'Colormap',gray(255))

% Next, load the table values under the schematic
sK = sprintf('%g',k);
switch schName
    case 'ZO'
        sZ = 'No zeros';
        sP = 'No poles';
        sQ = 'Not second order';
        sWp = 'No poles';
        sWz = 'No zeros';
        LFGain = abs(k);
        HFGain = abs(k);
        minGainFreq = 'All';
        minGain = abs(k);
        maxGainFreq = 'All';
        maxGain = abs(k);
        halfPowerFreq = 'None';
        halfPowerGain = 1/sqrt(2);
    case 'FO_LP'
        sZ = 'No zeros';
        sP = sprintf('%0.5g',p);
        sQ = 'Not second order';
        sWp = abs(p);
        sWz = 'No zeros';
        LFGain = abs(k/p);
        HFGain = 0;
        minGainFreq = 'Infinity';
        minGain = 0;
        maxGainFreq = 0;
        maxGain = abs(k/p);
        halfPowerFreq = -p;
        halfPowerGain = abs(k/(p*sqrt(2)));
    case 'FO_HP'
        sZ = '0';
        sP = sprintf('%0.5g',p);
        sQ = 'Not second order';
        sWp = abs(p);
        sWz = '0';
        LFGain = 0;
        HFGain = abs(k);
        minGainFreq = 0;
        minGain = 0;
        maxGainFreq = 'Infinity';
        maxGain = abs(k);
        halfPowerFreq = -p;
        halfPowerGain = 1/sqrt(2);
    case 'SO_LP'
        sZ = 'No zeros';
        sP = sprintf('%0.5g + j%0.5g',real(p(1)),abs(imag(p(1))));
        sQ = sprintf('%0.5g',-abs(p(1))/(2*real(p(1))));
        sWp = sprintf('%0.5g',abs(p(1)));
        sWz = 'No zeros';
        LFGain = abs(k)/abs(p(1))^2;
        HFGain = 0;
        minGainFreq = 'Infinity';
        minGain = 0;
        alpha2 = real(p(1))^2; sigma2 = imag(p(1))^2; alpha=-sqrt(alpha2); sigma=sqrt(sigma2);
        if sigma2>alpha2  % peaks
            maxGain = abs(k/(2*alpha*sigma));
            maxGainFreq = sqrt(sigma2-alpha2);
            halfPowerFreq = sqrt(sigma2-2*alpha*sigma-alpha2);
        else
            maxGain = abs(k)/(alpha2+sigma2);
            maxGainFreq = 0;
            halfPowerFreq = sqrt(sigma2-alpha2+sqrt(2)*sqrt(alpha2^2+sigma2^2));
        end
        halfPowerGain = maxGain/sqrt(2);
    case 'SO_HP'
        sZ = '0 + j0';
        sP = sprintf('%0.5g + j%0.5g',real(p(1)),abs(imag(p(1))));
        sQ = sprintf('%0.5g',-abs(p(1))/(2*real(p(1))));
        sWp = sprintf('%0.5g',abs(p(1)));
        sWz = '0';
        LFGain = 0;
        HFGain = abs(k);
        minGainFreq = 0;
        minGain = 0;
        alpha2 = real(p(1))^2; sigma2 = imag(p(1))^2; sigma=sqrt(sigma2); alpha=-sqrt(alpha2);
        if sigma2>alpha2 % peaks
            maxGainFreq = (sigma2+alpha2)/sqrt(sigma2-alpha2);
            maxGain = -abs(k)*(alpha2+sigma2)/(2*alpha*sigma);
            halfPowerFreq = (alpha2+sigma2)/sqrt(sigma2-2*alpha*sigma-alpha2);
        else
            maxGainFreq = 'Infinity';
            maxGain = abs(k);
            halfPowerFreq = sqrt(alpha2-sigma2+sqrt(2)*sqrt(sigma2^2+alpha2^2));
        end
        halfPowerGain = maxGain/sqrt(2);
    case 'SO_Notch'
        sZ = sprintf('%0.5g+j%0.5g',real(z(1)),abs(imag(z(1))));
        sP = sprintf('%0.5g+j%0.5g',real(p(1)),abs(imag(p(1))));
        sQ = sprintf('%0.5g',-abs(p(1))/(2*real(p(1))));
        sWp = sprintf('%0.5g',abs(p(1)));
        sWz = sprintf('%0.5g',abs(z(1)));
        LFGain = abs(k)*abs(z(1))^2/abs(p(1))^2;
        HFGain = abs(k);
        minGainFreq = abs(z(1));
        minGain = 0;
        alpha2 = real(p(1))^2; sigma2 = imag(p(1))^2; z2 = abs(z(1))^2; alpha = -sqrt(alpha2); sigma = sqrt(sigma2);
        if abs(p(1)) < abs(z(1)) % LP
            if sigma  <= -alpha % Q<1/2, so no hump
                maxGainFreq = 0;
                maxGain = abs(k*z2/(alpha2+sigma2));
            else % Q>1/2, so resonant hump
                maxGainFreq = sqrt((z2*(sigma2-alpha2)-(alpha2+sigma2)^2)/(z2+alpha2-sigma2));
                maxGain = sqrt(k^2*(z2^2+2*z2*(alpha-sigma)*(alpha+sigma)+(alpha2+sigma2)^2)/(4*alpha2*sigma2));
            end  
            a = alpha2-sigma2;
            b = alpha2+sigma2;
            c = alpha2^2+sigma2^2;
            d=-2*alpha2^2-4*alpha2*sigma2-2*sigma2^2-z2*a;
            e=2*z2*a*b^2+b^4+z2^2*c;
            halfPowerFreq = sqrt(z2)*sqrt((d+sqrt(2*e))/(z2^2-2*b^2));
            halfPowerGain = abs(k*z2/(alpha2+sigma2))/sqrt(2);     
        else % HP
            if sigma  <= -alpha % Q<1/2, so no hump
                maxGainFreq = 'Infinity';
                maxGain = abs(k);
            else % Q>1/2, so resonant hump
                maxGainFreq = sqrt((z2*(sigma2-alpha2)-(alpha2+sigma2)^2)/(z2+alpha2-sigma2));
                maxGain = sqrt(k^2*(z2^2+2*z2*(alpha-sigma)*(alpha+sigma)+(alpha2+sigma2)^2)/(4*alpha2*sigma2));
            end
            alpha4 = alpha2^2;
            z4 = z2^2;
            sigma4 = sigma2^2;
            halfPowerFreq = sqrt(2*z2+alpha2-sigma2+sqrt(2)*sqrt(z4+alpha4+sigma4+2*z2*(alpha2-sigma2)));
            halfPowerGain = abs(k)/sqrt(2);     
        end
    otherwise
        error('Unrecognized schName')
end
set(handles.uitxZ,'String',sZ)
set(handles.uitxP,'String',sP)
set(handles.uitxK,'String',sK)
set(handles.uitxQ,'String',sQ)
set(handles.uitxWo,'String',sWp)
set(handles.uitxWz,'String',sWz)

% fill in the table in the upper right
sLFGain = ConvertGain(LFGain,isDb);
sHFGain = ConvertGain(HFGain,isDb);
sMinGainFreq = ConvertFreq(minGainFreq,isHz);
sMinGain = ConvertGain(minGain,isDb);
sMaxGainFreq = ConvertFreq(maxGainFreq,isHz);
sMaxGain = ConvertGain(maxGain,isDb);
sHalfPowerFreq = ConvertFreq(halfPowerFreq,isHz);
sHalfPowerGain = ConvertGain(halfPowerGain,isDb);

set(handles.uitxLFGain,'String',sLFGain)
set(handles.uitxHFGain,'String',sHFGain)
set(handles.uitxMinGainFreq ,'String',sMinGainFreq)
set(handles.uitxMinGain ,'String',sMinGain)
set(handles.uitxMaxGainFreq ,'String',sMaxGainFreq)
set(handles.uitxMaxGain ,'String',sMaxGain)
set(handles.uitxHalfPowerFreq ,'String',sHalfPowerFreq)
set(handles.uitxHalfPowerGain ,'String',sHalfPowerGain)

% do the custom field
customMan = str2num(get(handles.uitxCustomFreq,'String'));
customExp = 3*(get(handles.uipmCustomFreq,'Value')-2);
customFreq = customMan*10^customExp*2*pi;
customGain = CalcGains(customFreq,z,p,k,isDb,0);
if isDb
    sCustomGain = sprintf('%g dB',customGain);
else
    sCustomGain = sprintf('%g',customGain);
end
set(handles.uitxCustomGain,'String',sCustomGain)

% create plot
if get(handles.uirbComponentsExact,'Value')==0 
    % if toleranced, use ideal values for finding FreqRange
    % so can compare to ideal without changes in axis
    [fMin,fMax]=CalcFreqRange(zIdeal,pIdeal,kIdeal,isHz); 
else
    [fMin,fMax]=CalcFreqRange(z,p,k,isHz); 
end
fRange = linspace(fMin,fMax,400);
gains = CalcGains(fRange,z,p,k,isDb,isHz);
% if toleranced components exist, figure out both tol and exact
% and set scale to the exact so that it won't flip between two
% nearly identical but different scalings
axes(handles.uiaxFreq)
plot(fRange,gains)
v = axis;
xMin = v(1); xMax = v(2); 
yMin = v(3); yMax = v(4);
v=[fMin fMax v(3) v(4)];
hold off
plot(fRange,gains,'b')
ylabel('|Gain|')
hold on
% draw the minGain lines
if ischar(minGain)
    minGain = str2num(minGain(1:4));
end
if ischar(minGainFreq)
    minGainFreq = str2num(minGainFreq(1:3));
end
if ~ischar(minGain) && ~ischar(minGainFreq)  && ~isempty(minGain) && ~isempty(minGainFreq)
    if isHz, minGainFreq = minGainFreq/(2*pi); end 
    if isDb 
        minGain = abs(minGain);
        minGain(minGain==0)=eps;
        minGain = 20*log10(minGain); 
    end
    if isfinite(minGain) && isfinite(minGainFreq) % draw a corner
        xMinGain = [xMin minGainFreq minGainFreq];
        yMinGain = [minGain minGain yMin];
    elseif isfinite(minGainFreq) % draw vertical line
        xMinGain = [minGainFreq minGainFreq];
        yMinGain = [yMin yMax];
    elseif isfinite(minGain) % draw horizontal line
        xMinGain = [xMin xMax];
        yMinGain = [minGain minGain];
    end
    plot(xMinGain,yMinGain,'g')
end
% draw the maxGain lines
if ischar(maxGain)
    maxGain = str2num(maxGain(1:4));
end
if ischar(maxGainFreq)
    maxGainFreq = str2num(maxGainFreq(1:3));
end
if ~ischar(maxGain) && ~ischar(maxGainFreq)  && ~isempty(maxGain) && ~isempty(maxGainFreq)
    if isHz, maxGainFreq = maxGainFreq/(2*pi); end
    if isDb
        if maxGain ==0
            maxGain = -inf;
        else
            maxGain = 20*log10(abs(maxGain)); 
        end
    end
    if isfinite(maxGain) && isfinite(maxGainFreq) % draw a corner
        xMaxGain = [xMin maxGainFreq maxGainFreq];
        yMaxGain = [maxGain maxGain yMin];
    elseif isfinite(maxGainFreq) % draw vertical line
        xMaxGain = [maxGainFreq maxGainFreq];
        yMaxGain = [yMin yMax];
    elseif isfinite(maxGain) % draw horizontal line
        xMaxGain = [xMin xMax];
        yMaxGain = [maxGain maxGain];
    end
    plot(xMaxGain,yMaxGain,'r')
end
% draw the half power lines
if ~ischar(halfPowerFreq)
    if isHz, halfPowerFreq = halfPowerFreq/(2*pi); end
    if isDb, halfPowerGain = 20*log10(halfPowerGain); end
    xHalfPower = [xMin halfPowerFreq halfPowerFreq];
    yHalfPower = [halfPowerGain halfPowerGain yMin];
    plot(xHalfPower,yHalfPower,'c')
end
% draw the custom frequency lines
if isHz, customFreq = customFreq/(2*pi); end
xCustomGain = [xMin customFreq customFreq];
yCustomGain = [customGain customGain yMin];
plot(xCustomGain,yCustomGain,'k')
% reset the axis
hold off
axis(v)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            CalculateSensitivity                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CalculateSensitivity(handles)
curStage = handles.strStage;
curStage = BuildCircuit_UpdateComponents(curStage,0,0,'New Tolerance');
newStage = curStage;
cComponents = get(handles.uipmSensitivityIn,'String');
sComponent = cComponents{get(handles.uipmSensitivityIn,'Value')};
isCap = NaN;
switch sComponent(1)
    case 'C', isCap = 1;
    case 'R', isCap = 0;
    otherwise, error('Unknown component')
end
isCalc = ~isempty(str2num(sComponent(2)));
if isCalc
    nComponent = str2num(sComponent(2));
else
    nComponent = double(sComponent(2))-96;
end

if ~isCalc % Selected components
    if isCap
        curCap = curStage.vfCSelect(nComponent);
        newStage.vfCSelect(nComponent) = curCap*1.01;
    else % isResistor
        curRes = curStage.vfRSelect(nComponent);
        newStage.vfRSelect(nComponent) = curRes*1.01;
    end
else % Calculated components
    if isCap
        curCap = curStage.vfCCalc(nComponent);
        newStage.vfCCalc(nComponent) = curCap*1.01;      
    else % isResistor
        curRes = curStage.vfRCalc(nComponent);
        newStage.vfRCalc(nComponent) = curRes*1.01;
    end
end
newStage = BuildCircuit_FindZ1P1K1(newStage);

switch get(handles.uipmSensitivityOut,'Value')
    case 1 % |z|
        if isempty(curStage.z) || curStage.z(1)==0
            S = 0;
        else
            oldZ = abs(curStage.z1(1));
            newZ = abs(newStage.z1(1));
            S = (newZ-oldZ)*100/oldZ;
        end
    case 2 % |p|
        if isempty(curStage.p)
            S = 0;
        else
            oldP = abs(curStage.p1(1));
            newP = abs(newStage.p1(1));
            S = (newP-oldP)*100/oldP;
        end
    case 3 %  k
        oldK = abs(curStage.k1);
        newK = abs(newStage.k1);
        S = (newK-oldK)*100/oldK;
    case 4 %  Q
        if isempty(curStage.Q1) || isnan(curStage.Q1) || ~isfinite(curStage.Q1) || curStage.Q1==0
            S = 0;
        else
            oldQ = abs(curStage.Q1);
            newQ = abs(newStage.Q1);
            S = (newQ-oldQ)*100/oldQ;
        end
    case 5 % wo
        if isempty(curStage.wp1) || isnan(curStage.wp1) || ~isfinite(curStage.wp1) || curStage.wp1==0  
            S = 0;
        else
            oldW = abs(curStage.wp1);
            newW = abs(newStage.wp1);
            S = (newW-oldW)*100/oldW;
        end
    case 6 % wz
        if isempty(curStage.wz1) || isnan(curStage.wz1) || ~isfinite(curStage.wz1) || curStage.wz1==0  
            S = 0;
        else
            oldW = abs(curStage.wz1);
            newW = abs(newStage.wz1);
            S = (newW-oldW)*100/oldW;
        end        
end
set(handles.uitxSensitivity,'String',num2str(S))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               ConvertGain                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sGain = ConvertGain(gain,isDb)
if ischar(gain)
    sGain = gain;
    return
end
if isDb
    if gain==0
        sGain = '-Infinity dB';
    else
        sGain = sprintf('%0.5g dB',20*log10(abs(gain)));
    end
else
    sGain = sprintf('%0.5g',gain);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               ConvertFreq                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sFreq = ConvertFreq(freq,isHz)
if ischar(freq)
    sFreq = freq;
    return
end
sSuffix = 'rad/s';
if isHz
    freq = freq / (2*pi);
    sSuffix = 'Hz';
end
sFreq = Utility_EngOutput(freq, sSuffix,5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               CalcFreqRange                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fMin,fMax]=CalcFreqRange(z,p,k,isHz)
if length(z) && z(1) == 0 % If HP
    fMin = abs(p(1))/4;
else
    fMin=0;
end
if length(p)==0 % ZO
    fMax = 10;
    return
end
if length(z)==0 % LP
    fMax = 2.5*abs(p(1));
elseif z(1)==0 % HP
    fMax = 3*max([abs(p(1)) abs(z(1))]);
elseif abs(p(1))< abs(z(1)) % LP Notch
    fMax = 2.5*abs(p(1));
else               % HP Notch
    fMax = 3*max([abs(p(1)) abs(z(1))]);
end
if isHz
    fMin = fMin/(2*pi);
    fMax = fMax/(2*pi);
end
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 CalcGains                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gains = CalcGains(f,z,p,k,isDb,isHz)
if isHz, f=f*2*pi; end
gains = zeros(size(f));
for i=1:length(f)
    w = f(i);
    gains(i) = k*prod(j*w-z)/prod(j*w-p);
end
gains = abs(gains);
gains(gains==0)=eps;
if isDb
    gains = 20*log10(gains);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             RemoveSmallValues                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = RemoveSmallValues(in)
TINY = 1e-10;
if isempty(in)
    out=[];
    return
end
for index=1:length(in)
    num=in(index);
    if abs(real(num))<TINY
        num=i*imag(num);
    end
    if abs(imag(num))<TINY
        num=real(num);
    end
    if real(num)~=0 && abs(imag(num)/real(num))<TINY
        num=real(num);
    end
    out(index)=num;
end