freqz2.m

来自「有关matlab的电子书籍有一定的帮助希望有用」· M 代码 · 共 122 行

function [hout,w1,w2] = freqz2(a,n1,n2,dx)
%FREQZ2 Compute two-dimensional frequency response.
%   [H,F1,F2] = FREQZ2(h,N1,N2) returns H, the N2-by-N1
%   frequency response of h, and the frequency vectors F1 (of
%   length N1) and f2 (of length N2). h is a two-dimensional FIR
%   filter, in the form of a computational molecule. F1 and F2
%   are returned as normalized frequencies in the range -1.0 to
%   1.0, where 1.0 corresponds to half the sampling frequency, or
%   pi radians.
%
%   [H,F1,F2] = FREQZ2(h,[N2 N1]) returns the same result as
%   [H,F1,F2] = FREQZ2(h,N1,N2).
%
%   [H,F1,F2] = FREQZ2(h) uses [N2 N1] = [64 64].
%
%   [H,F1,F2] = freqz2(h,F1,F2) returns the frequency response
%   for the FIR filter h at frequency values in F1 and F2. These
%   frequency values must be in the range -1.0 to 1.0, where 1.0
%   corresponds to half the sampling frequency, or pi radians.
%
%   [...] = FREQZ2(h,...,[DX DY]) uses [DX DY] to override the
%   intersample spacing in h. DX determines the spacing for the
%   x-dimension and DY determines the spacing for the
%   y-dimension. The default spacing is 0.5, which corresponds to
%   a sampling frequency of 2.0.
%
%   [...] = FREQZ2(h,...,DX) uses DX to determine the intersample
%   spacing in both dimensions.
%
%   With no output arguments, FREQZ2(...) produces a mesh plot of
%   the two-dimensional frequency response.
%
%   Class Support
%   -------------
%   The input matrix h can be of class uint8 or double. All other
%   inputs to freqz2 must be of class double. All outputs are of
%   class double.
%
%   Example
%   -------
%   Use the window method to create a 16-by-16 filter, then view
%   its frequency response using FREQZ2.
%
%       Hd = zeros(16,16);
%       Hd(5:12,5:12) = 1;
%       Hd(7:10,7:10) = 0;
%       h = fwind1(Hd,bartlett(16));
%       freqz2(h,[32 32]); axis([-1 1 -1 1 0 1]); colormap(jet(64))
%
%   See also FREQZ.

%   Clay M. Thompson 1-7-93
%   Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved.
%   $Revision: 5.8 $  $Date: 1997/11/24 15:34:39 $

if nargin<2, n1 = 64; n2 = 64; end     % Default number of frequency points
if nargin==2, 
  if all(size(n1)==[1 2]), 
    n2 = n1(1); n1 = n1(2);
  elseif prod(size(n1)==1)
    n2 = n1;
  else
    error('Not enough input arguments.');
  end
end
if nargin<4, dx = [0.5 0.5]; end % Default intersample spacing
if prod(size(dx))==1, dx = [dx dx]; end % Allow scalar dx

% Allow other numeric types for first input argument
if (~strcmp(class(a),'double'))
    a = double(a);
end

% Check for n1 and n2 < size(a).  If so, switch to w1,w2 based calc.
if (prod(size(n1))==1) & (prod(size(n2))==1), % N equally spaced points
  if any(size(a)>[n2 n1]),
    n1 =  ((0:n1-1)-floor(n1/2))*(1/(n1*dx(1)));
    n2 = (((0:n2-1)-floor(n2/2))*(1/(n2*dx(2))))';
  end
end

a = rot90(a,-2); % Unrotate filter since FIR filters are rotated.
if (prod(size(n1))==1) & (prod(size(n2))==1), % N equally spaced points
  % Pad A if necessary
  if any(size(a)~=[n2 n1]),
    [ma,na] = size(a);
    aa = zeros(n2,n1);
    aa(floor(n2/2)-floor(ma/2)+(1:ma),floor(n1/2)-floor(na/2)+(1:na)) = a;
    a = aa;
  end
  a = rot90(fftshift(rot90(a,2)),2); % Inverse fftshift
  w1 =  ((0:n1-1)-floor(n1/2))*(1/(n1*dx(1)));
  w2 = (((0:n2-1)-floor(n2/2))*(1/(n2*dx(2))))';
  h = fftshift(fft2(a));

else % evaluate frequency response at points in n1.
  % Expand w1 and w2 if they are vectors
  if min(size(n1))==1 | min(size(n2))==1, 
    [w1,w2] = meshgrid(n1,n2);
  else
    w1 = n1; w2 = n2;
  end
  [t1,t2] = freqspace(size(a)); 
  t1 = t1*size(a,2)/2; 
  t2 = t2*size(a,1)/2;
  [t1,t2] = meshgrid(t1,t2);
  h = zeros(size(w1));
  h(:) = (exp(-sqrt(-1)*pi*w1(:)*t1(:)').*exp(-sqrt(-1)*pi*w2(:)*t2(:)'))*a(:);
end

% Convert to real if possible
if all(max(abs(imag(h)))<sqrt(eps)), h = real(h); end

if nargout==0, % Plot data
  mesh(w1,w2,abs(h))
  xlabel('Frequency'), ylabel('Frequency'), zlabel('Magnitude')
  return
end
hout = h;
w1 = w1(1,:)';
w2 = w2(:,1);