Accelerating the pace of engineering and science

# Documentation

## Frequency Response

### Digital Domain

freqz uses an FFT-based algorithm to calculate the z-transform frequency response of a digital filter. Specifically, the statement

```[h,w] = freqz(b,a,p)
```

returns the p-point complex frequency response, H(ejω), of the digital filter.

$H\left({e}^{jw}\right)=\frac{b\left(1\right)+b\left(2\right){e}^{-jw}+...+b\left(n+1\right){e}^{-jwn}}{a\left(1\right)+a\left(2\right){e}^{-jw}+...+a\left(m+1\right){e}^{-jwm}}$

In its simplest form, freqz accepts the filter coefficient vectors b and a, and an integer p specifying the number of points at which to calculate the frequency response. freqz returns the complex frequency response in vector h, and the actual frequency points in vector w in rad/s.

freqz can accept other parameters, such as a sampling frequency or a vector of arbitrary frequency points. The example below finds the 256-point frequency response for a 12th-order Chebyshev Type I filter. The call to freqz specifies a sampling frequency fs of 1000 Hz:

```[b,a] = cheby1(12,0.5,200/500);
[h,f] = freqz(b,a,256,1000);
```

Because the parameter list includes a sampling frequency, freqz returns a vector f that contains the 256 frequency points between 0 and fs/2 used in the frequency response calculation.

 Note   This toolbox uses the convention that unit frequency is the Nyquist frequency, defined as half the sampling frequency. The cutoff frequency parameter for all basic filter design functions is normalized by the Nyquist frequency. For a system with a 1000 Hz sampling frequency, for example, 300 Hz is 300/500 = 0.6. To convert normalized frequency to angular frequency around the unit circle, multiply by π. To convert normalized frequency back to hertz, multiply by half the sample frequency.

If you call freqz with no output arguments, it plots both magnitude versus frequency and phase versus frequency. For example, a ninth-order Butterworth lowpass filter with a cutoff frequency of 400 Hz, based on a 2000 Hz sampling frequency, is

```[b,a] = butter(9,400/1000);
```

To calculate the 256-point complex frequency response for this filter, and plot the magnitude and phase with freqz, use

```freqz(b,a,256,2000)
```

freqz can also accept a vector of arbitrary frequency points for use in the frequency response calculation. For example,

```w = linspace(0,pi);
h = freqz(b,a,w);
```

calculates the complex frequency response at the frequency points in w for the filter defined by vectors b and a. The frequency points can range from 0 to 2π. To specify a frequency vector that ranges from zero to your sampling frequency, include both the frequency vector and the sampling frequency value in the parameter list.

### Analog Domain

freqs evaluates frequency response for an analog filter defined by two input coefficient vectors, b and a. Its operation is similar to that of freqz; you can specify a number of frequency points to use, supply a vector of arbitrary frequency points, and plot the magnitude and phase response of the filter.

### Magnitude and Phase

MATLAB® functions are available to extract magnitude and phase from a frequency response vector h. The function abs returns the magnitude of the response; angle returns the phase angle in radians. To extract the magnitude and phase of a Butterworth filter:

```[z,p,k] = butter(9,400/1000);
fvtool(zp2sos(z,p,k))
```

and click the Magnitude and Phase Response button on the toolbar or select Analysis > Magnitude and Phase Response to display the plot.

The unwrap function is also useful in frequency analysis. unwrap unwraps the phase to make it continuous across 360º phase discontinuities by adding multiples of ±360°, as needed. To see how unwrap is useful, design a 25th-order lowpass FIR filter:

```h = fir1(25,0.4);
```

Obtain the filter's frequency response with freqz, and plot the phase in degrees:

```[H,f] = freqz(h,1,512,2);
plot(f,angle(H)*180/pi)
grid```

It is difficult to distinguish the 360° jumps (an artifact of the arctangent function inside angle) from the 180° jumps that signify zeros in the frequency response.

unwrap eliminates the 360° jumps:

```plot(f,unwrap(angle(H))*180/pi)
```

Alternatively, you can use phasez to see the unwrapped phase:

`phasez(h,1)`

### Delay

The group delay of a filter is a measure of the average time delay of the filter as a function of frequency. It is defined as the negative first derivative of a filter's phase response. If the complex frequency response of a filter is H(e), then the group delay is

${\tau }_{g}\left(\omega \right)=-\frac{d\theta \left(\omega \right)}{d\omega }$

where θ(ω) is the phase, or argument of H(e). Compute group delay with

```[gd,w] = grpdelay(b,a,n)
```

which returns the n-point group delay, ,τg(ω) of the digital filter specified by b and a, evaluated at the frequencies in vector w.

The phase delay of a filter is the negative of phase divided by frequency:

${\tau }_{p}\left(\omega \right)=-\frac{\theta \left(\omega \right)}{\omega }$

To plot both the group and phase delays of a system on the same FVTool graph, type

```[z,p,k] = butter(10,200/1000);
fvtool(zp2sos(z,p,k),'Analysis','grpdelay', ...
'OverlayedAnalysis','phasedelay','Legend','on')```