View Raw SPL
/*****************************************************************************
* *
* FREQZ.SPL Copyright (C) 2009 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Evaluates the frequency response of a Z transform *
* *
* Revisions: 17 Aug 2009 RRR Creation *
* *
*****************************************************************************/
#include
#if @HELP_FREQZ
FREQZ
Purpose: Evaluates the frequency response of a Z transform
Syntax: FREQZ(b, a, N, "whole", Fs)
(h, w) = FREQZ(b, a, N, "whole", Fs)
b - a series, numerator (i.e. zero) coefficients
a - a series, denominator (i.e. pole) coefficients,
if the first coefficient is not 1.0, the coefficents
are assumed to be in difference equation form
N - an optional integer, number of output samples,
defaults to 512
"whole" - an optional string, if specified, the transform is
evaluated over the entire unit circle.
Fs - an optional real, sample rate of data. If not specified,
the transform is evaluated in normalized angular frequency
with values of pi * radians / s.
Alternative Syntax:
FREQZ(b, a, w, "whole", Fs)
(h, w) = FREQZ(b, a, w, "whole", Fs)
b - a series, numerator (i.e. zero) coefficients
a - a series, denominator (i.e. pole) coefficients,
if the first coefficient is not 1.0, the coefficents
are assumed to be in difference equation form
w - a series, the frequencies at which to evaluate the
transform.
"whole" - an optional string, if specified, the transform is
evaluated over the entire unit circle.
Fs - an optional real, sample rate of data. If the rates of
the numerator and denominator coefficients are equal,
the rate defaults to the coefficient rate, else the rate
defaults to 1.0.
Returns: Displays the magnitude and phase response in two windows.
h = FREQZ(b, a) returns the complex frequency response
as one interval or XY series.
(h, w) = FREQZ(b, a) returns the complex frequency response
and frequency values as two separate series.
Example:
h = freqz({1}, {1, -0.5, 0.8})
h contains 512 uniformly spaced samples of the frequency
response of the Z transform:
1
H(z) = ------------------
-1 -2
1 - 0.5z + 0.8z
The angular frequency samples range from 0 to pi radians/s.
Example:
freqz({1}, {1, -0.5, 0.8})
Same as the first example except the magnitude and phase
responses are displayed in two separate Windows.
The displayed normalized frequency samples range from
0 to 1.0 in units of pi radians/sec.
Example:
freqz({1}, {1, 0.5, -0.8}, 1024, 500)
Same as the previous example except 1024 samples are computed
and the frequency values range from 0 to 250 Hertz. The
magnitude and phase responses are automatically displayed in
two separate Windows.
Remarks:
If no explicit frequencies are specified, FREQZ uses the
FFT method to evaluate N uniformly spaced samples over the
unit circle of a Z transform in direct form:
-1 -Q
Y(z) b[1] + b[2]*z + ... + b[Q+1]*z
H(z) = ---- = -------------------------------------
-1 -P
X(z) 1 + a[2]*z + ... + a[P+1]*z
If no sample rate is specified, the angular frequencies are
used in units of radians/s.
If a sample rate is specified, the frequencies are
determined in Hertz.
If no output arguments are specified, the magnitude
response in dB and the phase response in degrees are
displayed in two Windows. If no sample rate is specified,
the displayed frequencies are in normalized angular units
of pi * radians/s.
h = FREQZ returns a complex series. Use MAG or PHASE to obtain
the magnitude and/or phase components separately.
See ZFREQ to specify the input coefficients in combined direct
form or cascaded bi-quad form.
See Also:
Filteq
Freqs
Invfreqs
Invfreqz
Mag
Phase
Residue
Residuez
Sfreq
Zfreq
#endif
/* evaluate the frequency response of a discrete system */
freqz(b, a, N, whole, Fs)
{
local w, srate;
/* parse args */
(b, a, N, whole, Fs) = freqz_parse_args(b, a, N, whole, Fs);
if (length(N) > 1)
{
/* frequency vector */
(h, w) = freqz_vector(b, a, N, whole, Fs);
}
else
{
/* use FFT method */
srate = isempty(Fs) ? 1 : Fs;
/* let ZFREQ do the calculation without assuming difference eq form */
(h, w) = zfreq(b, a, N, srate, whole, 0);
}
/* convert to radian frequency */
if (isempty(Fs))
{
if (outargc == 0)
{
w *= 2;
setvunits(w, "$\pi$ radians/s");
}
else
{
w *= 2 * pi;
setvunits(w, "radians/s");
sethunits(h, "radians/s");
setdeltax(h, deltax(h) * 2 * pi);
}
}
else
{
setvunits(w, "Hertz");
sethunits(h, "Hertz");
}
if (outargc == 0)
{
freqzplot(h, w, Fs);
}
else if (outargc > 1)
{
return(h, w);
}
else
{
if (length(N) > 1)
{
return(xy(w, h));
}
else
{
return(h);
}
}
}
freqz_vector(b, a, f, whole, Fs)
{
local w, h, c = 1;
if (not(isempty(Fs)))
{
w = 2 * pi * f / Fs;
}
else
{
w = f;
c = 2 * pi;
}
w = exp(-i*w);
h = polyval(rev(b), w) / polyval(rev(a), w);
return(h, f / c);
}
/* parse input arguments to freqz */
freqz_parse_args(b, a, N, whole, Fs)
{
local tmp;
switch (argc)
{
case 0:
error("freqz - input coefficients required");
break;
case 1:
a = {};
case 2:
N = {};
case 3:
whole = {};
case 4:
Fs = {};
case 5:
break;
default:
error("freqz - too many arguments");
break;
}
if (isstring(N))
{
whole = 1;
n = {};
}
else if (isstring(Fs))
{
tmp = whole;
whole = 1;
Fs = tmp;
}
else if (isstring(whole))
{
whole = 1;
}
else
{
Fs = whole;
whole = 0;
}
if (isempty(N))
{
N = 512;
}
return(b, a, N, whole, Fs);
}