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