gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\freqs.m

    function [h,ww] = freqs(b,a,w)
%FREQS Laplace-transform (s-domain) frequency response.  
%   H = FREQS(B,A,W) returns the complex frequency response vector H 
%   of the filter B/A:
%                        nb-1      nb-2
%   	         B(s)   b(1)s  +  b(2)s   + ... +  b(nb)
%     H(s) = ---- = --------------------------------
%                        na-1      na-2
%   	         A(s)   a(1)s  +  a(2)s   + ... +  a(na)
%
%   given the numerator and denominator coefficients in vectors B and A. 
%   The frequency response is evaluated at the points specified in 
%   vector W.  The magnitude and phase can be graphed by calling
%   FREQS(B,A,W) with no output arguments.
%
%   [H,W] = FREQS(B,A) automatically picks a set of 200 frequencies W on 
%   which the frequency response is computed.  FREQS(B,A,N) picks N 
%   frequencies. 
%
%   See also LOGSPACE, POLYVAL, INVFREQS, and FREQZ.

% 	Author(s): J.N. Little, 6-26-86
%   	   T. Krauss, 3-19-93, default plots and frequency vector
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 1998/07/20 18:53:02 $

if nargin == 2,
    w = 200;
end
if length(w) == 1,
    wlen = w;
    w_long = freqint(b,a,wlen);
    % need to interpolate long frequency vector:
    xi = linspace(1,length(w_long),wlen).';
    w = 10.^interp1(1:length(w_long),log10(w_long),xi,'linear');
end

s = j*w;
hh = polyval(b,s) ./ polyval(a,s);

if nargout == 0,
    newplot;
    mag = abs(hh);   phase = angle(hh)*180/pi;
    subplot(211),loglog(w,mag),set(gca,'xgrid','on','ygrid','on') 
    set(gca,'xlim',[w(1) w(length(w))])
    xlabel('Frequency (radians)')
    ylabel('Magnitude')
    ax = gca;
    subplot(212), semilogx(w,phase),set(gca,'xgrid','on','ygrid','on') 
    set(gca,'xlim',[w(1) w(length(w))])
    xlabel('Frequency (radians)')
    ylabel('Phase (degrees)')
    axes(ax)
elseif nargout == 1,
    h = hh;
elseif nargout == 2,
    h = hh;
    ww = w;
end
% end freqs

function w=freqint(a,b,c,d,npts)
%FREQINT Auto-ranging algorithm for Bode frequency response
%   W=FREQINT(A,B,C,D,Npts)
%   W=FREQINT(NUM,DEN,Npts)

%   Andy Grace  7-6-90
%   Was Revision: 1.9,  Date: 1996/07/25 16:43:37

% Generate more points where graph is changing rapidly.
% Calculate points based on eigenvalues and transmission zeros. 

[na,ma] = size(a);

if (nargin==3)&(na==1),		% Transfer function form.
  npts=c;
  ep=roots(b);
  tz=roots(a);
else				% State space form
  if nargin==3, npts=c; [a,b,c,d] = tf2ss(a,b); end
  ep=[eig(a)];
  tz=tzero(a,b,c,d);
end

if isempty(ep), ep=-1000; end

% Note: this algorithm does not handle zeros greater than 1e5
ez=[ep(find(imag(ep)>=0));tz(find(abs(tz)<1e5&imag(tz)>=0))];

% Round first and last frequencies to nearest decade
integ = abs(ez)<1e-10; % Cater for systems with pure integrators
highfreq=round(log10(max(3*abs(real(ez)+integ)+1.5*imag(ez)))+0.5);
lowfreq=round(log10(0.1*min(abs(real(ez+integ))+2*imag(ez)))-0.5);

% Define a base range of frequencies
diffzp=length(ep)-length(tz);
w=logspace(lowfreq,highfreq,npts+diffzp+10*(sum(abs(imag(tz))<abs(real(tz)))>0));
ez=ez(find(imag(ez)>abs(real(ez))));

% Oscillatory poles and zeros
if ~isempty(ez)
  f=w;
  npts2=2+8/ceil((diffzp+eps)/10);
  [dum,ind]=sort(-abs(real(ez)));
  z=[];
  for i=ind'
    r1=max([0.8*imag(ez(i))-3*abs(real(ez(i))),10^lowfreq]);
    r2=1.2*imag(ez(i))+4*abs(real(ez(i)));
    z=z(find(z>r2|z<r1));
    indr=find(w<=r2&w>=r1);
    f=f(find(f>r2|f<r1));
    z=[z,logspace(log10(r1),log10(r2),sum(w<=r2&w>=r1)+npts2)];
  end
  w=sort([f,z]);
end
w = w(:);

% end freqint