freqz_matlab求频响_ahead6yh_whole_freqz函数matlab_

function [hh,ww] = freqz(varargin) %FREQZ Z-transform digital filter frequency response. % When N is an integer, [H,W] = FREQZ(B,A,N) returns the N-point frequency % vector W in radians and the N-point complex frequency response vector H % of the filter B/A: % jw -jw -jnbw % jw B(e) b(1) + b(2)e + .... + b(nb+1)e % H(e) = ---- = ---------------------------- % jw -jw -jnaw % A(e) a(1) + a(2)e + .... + a(na+1)e % given numerator and denominator coefficients in vectors B and A. The % frequency response is evaluated at N points equally spaced around the % upper half of the unit circle. If N isn't specified, it defaults to 512. % % [H,W] = FREQZ(B,A,N,'whole') uses N points around the whole unit circle. % % H = FREQZ(B,A,W) returns the frequency response at frequencies % designated in vector W, in [radians/sample] (normally between 0 and pi). % % [H,F] = FREQZ(B,A,N,Fs) and [H,F] = FREQZ(B,A,N,'whole',Fs) given a % sampling freq Fs in Hz return a frequency vector F in Hz. % % H = FREQZ(B,A,F,Fs) given sampling frequency Fs in Hz returns the % complex frequency response at the frequencies designated in vector F, % also in Hz. % % FREQZ(B,A,...) with no output arguments plots the magnitude and % unwrapped phase of B/A in the current figure window. % % See also FILTER, FFT, INVFREQZ, FREQS and GRPDELAY. % CAUTION, parameter-value pairs can be used, but this feature will be % removed in the next release. % Author(s): J.N. Little, 6-26-86 % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.25.1.2 $ $Date: 1999/01/22 03:42:34 $ % ------------------------- % % Parse the input arguments % % ------------------------- % [params,fvflag] = freqzparse(varargin{:}); % fvflag will be set to one if a freqvector was specified by the user. b = params.numerator(:).'; % Make it a row vector a = params.denominator(:).'; % fs = params.fs; nb = length(b); na = length(a); % ------------------------------------- % % Actual Frequency Response Computation % % ------------------------------------- % if fvflag, % Frequency vector specified. Use Horner's method of polynomial % evaluation at the frequency points and divide the numerator % by the denominator. % % Note: we use positive i here because of the relationship % polyval(a,exp(i*w)) = fft(a).*exp(i*w*(length(a)-1)) % ( assuming w = 2*pi*(0:length(a)-1)/length(a) ) % a = [a zeros(1,nb-na)]; % Make sure a and b have the same length b = [b zeros(1,na-nb)]; [digw,xlab,xtickflag,xlim] = freqconv(params.freqvector,fs,... params.frequency,params.normfreq,params.range,'freq2dig'); % Scale w correctly s = exp(i*digw); % Digital frequency must be used for this calculation h = polyval(b,s) ./ polyval(a,s); w = params.freqvector; % Assign the correct freqvector to w else % freqvector not specified, use nfft and RANGE in calculation s = strmatch(lower(params.range),{'whole','half'}); n = params.nfft; if s*n < na | s*n < nb nfft = lcm(n,max(na,nb)); % dividenowarn temporarily shuts off warnings to avoid "Divide by zero" h = dividenowarn(fft([b zeros(1,s*nfft-nb)]),... fft([a zeros(1,s*nfft-na)])).'; h = h(1+(0:n-1)*nfft/n); else % dividenowarn temporarily shuts off warnings to avoid "Divide by zero" h = dividenowarn(fft([b zeros(1,s*n-nb)]),... fft([a zeros(1,s*n-na)])).'; h = h(1:n); end h = h(:); % Make it a column only when nfft is given (backwards comp.) [w,xlab,xtickflag,xlim] = freqconv(n,fs,params.frequency,... params.normfreq,params.range,'dig2freq',params.return_nyquist); % Scale w correctly w = w(:); % Make it a column only when nfft is given (backwards comp.) end % ------------------------------ % % Plot or Compute the Output % % ------------------------------ % if nargout == 0, % Plot when no output arguments are given % If it's normalized angular frequency divide by pi % because the plot's x domain is 0 to 1. if strmatch(params.normfreq,'yes') & strmatch(params.frequency,'angular'), w = w./pi; end freqzplot(h,w,params,fvflag,xlim,xlab,xtickflag); else % Don't plot, just return the (complex) frequency response. hh = h; ww = w; end % ------------------------------ % % Local Functions % % ------------------------------ % %------------------------------------------------------------------------------ function freqzplot(h,w,params,fvflag,xlim,xlab,xtickflag) % FREQZPLOT Plots the Magnitude and Phase response. % % Inputs: % h - complex frequency response vector % w - frequency vector % params - structure containing the following fields % normmag - YES/NO flag indicating if magnitude is normalized % magnitude - magnitude units string: 'DB', 'LINEAR' or 'SQUARED' % phase - YES/NO flag indicating if phase plot is required % range - frequency range, 'WHOLE' or 'HALF', which corresponds to % [0,Fs) or [0,Fs/2) respectively % fvflag - flag indicating if a frequency vector was specified by user % xlim - x-limits which corresponds to correct x-axis units % xlab - x-axis label with correct units % xtickFlag - specifies the frequency units (angular, normalized angular, % linear, etc.) newplot; if ishold, holdflag = 1; else holdflag = 0; end % Check if we should scale the magnitude, this is for plotting only [magh,ylab] = correctmag(h,params.normmag,params.magnitude); if strmatch(lower(params.phase),'yes'), subplot(212); % We plot the phase first to retain the functionality of freqz when hold is on plot(w,unwrap(angle(h))*180/pi); ax = gca; xlabel(xlab); ylabel('Phase (degrees)'); subplot(211); plot(w,magh); xlabel(xlab); ylabel(ylab); ax = [ax gca]; if ~holdflag, hfig = get(ax(1),'parent'); set(hfig,'nextplot','replace'); % Resets the figure so that next plot does not subplot end else, plot(w,magh); xlabel(xlab); ylabel(ylab); ax = gca; end axes(ax(1)); % Always bring the plot to the top set(ax,'xgrid','on','ygrid','on'); if fvflag, xlim = [w(1) w(end)]; % If a freqvector was given, overwrite the xlim returned by freqconv end % If hold is on, check if xlim in current plot is greater than the new calculated xlim if holdflag, currxlim = get(ax(end),'xlim'); xlim(1) = min([currxlim(1) xlim(1)]); xlim(2) = max([currxlim(2) xlim(2)]); end set(ax,'xlim',xlim); if ~fvflag, % If xtickflag='normang' sets the ticks of x-axis in units of pi. %setxticks(ax,params.range,xtickflag); end %------------------------------------------------------------------------------ function [magh,ylab] = correctmag(h,normmag,magnitude) % CORRECTMAG scale the magnitude correctly according to given parameters. % CORRECTMAG will take the precomputed frequency response vector, h, % and calculate its magnitude according to the desired parameters. % MAGNITUDE can be either 'DB','LINEAR' or 'SQUARED' and NORMMAG can be % either 'YES' (divide everything by maximum value) or 'NO'. % [MAGH,YLAB] = correctmag(...) returns the correct y-axis label % YLAB. magh = abs(h); ylab = 'Magnitude'; if strmatch(lower(normmag),'yes'), magh = magh./max(magh); ylab = ['Normalized ' ylab]; end if strmatch(lower(magnitude),'db'), % These next few lines of code are here to avoid "Log of zero" warnings zerosIndx = find(magh==0); % Find the indicies where magnitude is zero magh(zerosIndx) = 1; % Place hold