\(\)Similar to the Analog RLC Filters, discrete filters can remove low and/or high frequencies to suppress unwanted frequencies such as background noise. This article shows the math behind these filters and visualizes their system response.

General .. FIR filters

See ccrma.stanford.edu for more.

For a start we will start with an abbreviated description of the filter discussed in [Smith, J.O. Introduction to Digital Filters with Audio Applications, http://ccrma.stanford.edu/~jos/filters/, online book, 2007 edition, accessed 2017-11-24].

The simplest low-pass filter is given by the difference equation

The signal flow graph for this filter is

## Euler to the rescue

We could study the system response to a general sinusoid input \(x[n]=A\cos(\omega nT+\phi)\). To derive the gain and phase of this simple filter, we need over a page of equations using various trig identities.

We can avoid the trig identities by using Euler’s formula and its trigonometry connection

$$ \begin{align} \mathrm{e}^{j\varphi}&=\cos\varphi+j\sin\varphi\nonumber\\ \cos\varphi&=\frac{\mathrm{e}^{j\varphi}+\mathrm{e}^{-j\varphi}}{2}\nonumber \end{align} \nonumber $$

Consider a complex sinusoid input \(x[n]\)

We can simplify equation \(\eqref{eq:xn}\):

- Because of linearity, input \(\mathrm{e}^{j\omega nT}\) is equivalent to putting \(\cos(\omega nT)\) and \(j\sin(\omega nT)\) in different copies of the filter and adding the outputs together. To study the filter, we may input the imaginary part, for as long as we ignore the associated imaginary output.
- The filter gain doesn’t depend on the amplitude \(A\) due to linearity.
- The filter phase response does not depend on \(\phi\) due to time-invariance.
- choose frequencies \(\omega_1\ldots\omega_L\) at which the frequency response should contain a null, a dip or a peak.
- Choose zeros \(r_i=|z|\mathrm{j\omega_i}\) at those frequencies \(\neq0\) at which a null or dip should occur. With \(|z|=1\), or close to \(1\) as desired. For each such \(\omega_i\), choose also a zero \(r_j\) that is the complex conjugate of \(r_i\). Let \(M\) be the total number of zeros chosen.
- Choose poles \(p_i=|z|\mathrm{j\omega_i}\) at those frequencies \(\neq0\) at which a null or dip should occur. With \(|z|=1\), or close to \(1\) as desired. For each such \(\omega_i\), choose also a zero \(r_j\) that is the complex conjugate of \(r_i\). Let \(N\) be the total number of poles chosen.
- Form the system function \(H(z)=K\frac{(1-r_1z^{-1})(1-r_2z^{-1})\ldots(1-r_Mz^{-1})}{(1-p_1z^{-1})(1-p_2z^{-1})\ldots(1-p_Nz^{-1})}\), where \(K\) is the desired gain.
- Cross multiply these factors of \(H(z)\), to express \(H(z)\) as the ratio of two polynomials whose terms are powers of \(z^{-1}\).
- Identify the IIR filter coefficients {a_0,\ldots,a_N,b_0,\ldots,b_M}, which are simply the coefficients of the polynomials.
- src: http://www.eecs.umich.edu/courses/eecs206/archive/spring02/lab.dir/Lab9/lab9_v3_0_release.pdf

With letting \(A=1\) and \(\phi=0\), and promising to ignore the imaginary part of the output, equation \(\eqref{eq:xn}\) can simplified to

Substituting \(\eqref{eq:xn2}\) in \(\eqref{eq:difeq}\) gives the output in response to the complex sinusoid input

### Transfer function

In equation \(\eqref{eq:yn}\), the output \(y[n]\) is the input signal \(x[n]\) multiplied with what we will define as the filter’s transfer function \(H(\omega)\).

Balance the exponents, so that we can apply Euler’s trig connection to \(\cos()\)

Based on Euler’s formula, we can express \(H(\omega)\) in polar coordinates, where \(G(\omega)\) is the gain (amplitude response), and \(\varphi(\omega)\) is the phase response

#### Gain and Phase

The gain \(G(\omega)\) follows

The phase response \(\varphi(\omega)\) follows as

Combining these we get the filter’s amplitude and phase response

The gain goes sinusoidally from 2 to 0 as \(\omega T\) goes from \(0\) to \(\pi\). It is somewhat of a low-pass filter, but the slope of the filter is more gently as we usually would like to see in a low-pass filter. Then again, this filter was very simple and would be fast to implement.

The output signal y[n] follows from \(\eqref{eq:xn}\)

Like any LTI filter, it can change the amplitude and phase of an incoming sinusoid signa.

Continue at ccrma.stanford.edu.

## IIR

Here we describe a single-pole low-pass infinite impulse response (IIR) filter. In IIR filters the output depends on the past and current value of the input signal and past values of the output signal

See ccrma.stanford.edu.

The general form for IIR filters is

A single-pole low-pass infinite impulse response (IIR) filter is given by the Z-Transform

A low-pass single-pole IIR filter as a single design parameter, the decay value \(d\). It is customary to define parameters \(a=d\) and \(b=1-d\). For a typical value of \(d=0.99\), we get \(a=0.99\) and \(b=0.01\). The relation is given by

This relation directly shows the effect of the filter. The previous output value of the filter \(y[n-1]\) is decreased with the decay factor \(a\), and a small fraction \(b\) of the current input \(x[n]\) is taken into account.

Substituting \(b=1-a\) and rewriting gives the expression

In the programming code, this is expressed as \(y += b * (x – y)\)

The response of this filter is completely analogous to the response of an analog low-pass filter containing a resistor and a capacitor.

Decay is related to the time constant \(\tau\) of the filter.

The cutoff frequency \(f_c\) is also related to \(d\) as

Further look at IIR filters:

Note decay and cutoff frequency formulas.

One of the simplest forms of passive filters consists of a resistor and capacitor in series The output is the voltage over the capacitor \(y(t)\) as shown in the schematic below.

This type of filter is called an Infinite-Impulse Response (IIR) filter, because if you give it an impulse input, the output takes an infinite time to go down to *exactly* zero.

Even though this article shows a low pass filter, the same principles apply to a high pass filter where the output is taken over the resistor.

We will derive the transfer function for this filter and determine the step and frequency response functions. Required prior reading includes Laplace Transforms, Impedance and Transfer Functions.

## Transfer Function

In the RC circuit, shown above, the current is the input voltage divided by the sum of the impedance of the resistor \(Z_R=R\) and that of the capacitor \(Z_C\). The output is the voltage over the capacitor and equals the current through the system multiplied with the capacitor impedance. The transfer function follows as the quotient of the output and input signals.

The denominator of \(\eqref{eq:voltagedivider}\) is a first-order polynomial. The root of this polynomial is called the system’s pole

This first order system \(\eqref{eq:transferpolynomial}\) has no zeros and one stable pole \(p\) on the left real axis \(p<0\) as visualized in the \(s\)-plane.

## Unit Step Response

The step response gives an impression of the system behavior when the input signal going from \(0\) to \(1\) volt at time \(t=0\). This input is called the Unit Step Function, here represented by \(u(t)=\gamma(t)\).

Combining the Laplace transform of \(\gamma(t)\), \(\Gamma(s)\), with the transfer function \(\eqref{eq:transferpolynomial}\), gives the unit step response \(Y(s)\)

Split up this complicated fraction into forms that are in the Laplace Transform table. According to Heaviside, this can be expressed as partial fractions. [swarthmore, MIT-cu]

Substitute \(K=-p\) from \(\eqref{eq:transferpolynomial}\) and find expressions for the constants \(c_{0,1}\), by multiplying with respectively \(s\) and \((s-p)\)

Given that these equations are true for any value of \(s\), choose two convenient values of \(s\) that help us find \(c_0\) and \(c_1\).

The unit step response \(y(t)\) follows from the inverse Laplace transform of \(\eqref{eq:heaviside}\)

Substituting the constants \(\eqref{eq:constants1a}, \eqref{eq:constants1b}\) gives the unit step response

As shown in the graph, the unit step response is a relatively slow decaying exponential curve (with \(p\lt 0\)).

The GNU Octave code used to print this figure is listed in the Appendix.

## Frequency Response

The frequency response \(y_{ss}(t)\) is defined as the steady state response to a sinusoidal input signal \(u(t)=\sin(\omega t)\,\gamma(t)\). It describes how well the filter can distinguish between different frequencies.

In Evaluating Transfer Functions, we have proven that

The transfer function \(H(s)\) for this RC Filter is given by \(\eqref{eq:transferpolynomial}\)

The system behavior at \(\omega=0\) and at \(\omega\implies \infty\) indicates that this is a low pass filter.

Based on Euler’s formula, we can express \(H(s)\) in polar coordinates

This transfer function with pole \(p\), evaluated for \(s=j\omega\) can be visualized with a vector from the pole to \(j\omega\).

The length of the vector corresponds to \(|(H(j\omega)|\), and minus the angle with the real axis corresponds to phase shift \(\angle H(j\omega)\).

Substitute \(p=-\frac{1}{RC}\) and \(K=\frac{1}{RC}\)

The output signal \(y_{ss}(t)\) for a sinusoidal input signal \(\sin(\omega t)\,\gamma(t)\) $$ \shaded{ \begin{aligned} y_{ss}(t)&=|H(j\omega)|\,\sin(\omega t+\angle H(j\omega))\,\gamma(t)\nonumber \\ |H(j\omega)| &=\frac{1}{\sqrt{(1+\omega RC)^2}}\\ \angle{H(j\omega)}&=-\arctan(\omega RC)\nonumber \end{aligned} } \label{eq:frequencyresponse} $$

This frequency response for different frequencies can be visualized in a Bode plot or a Nyquist diagram. Each of these are a topic of the remaining sections.

### Effect on Input with Harmonics

As a side step, we examine the effect of the filter on a square wave input signal. The Fourier series of the square wave shows that it consists of a base frequency and odd harmonics.

Substituting \(C=470\ \mathrm{nF}\), \(R=100\ \Omega\) and 20 kHz in \(\eqref{eq:frequencyresponse}\), gives the output signal \(y_{ss}(t)\)

UNFINISHED

### Bode plot

A Bode plots frequency as the horizontal axis and usually consists of two separate plots to that show the magnitude and phase of the frequency response \(y_{tt}\). Since the range of magnitudes may also be large, the amplitude scale is usually expressed in decibels \(20\log_{10}\left|H(j\omega)\right|\). The frequency axis uses a logarithmic scale as well.

The magnitude of the frequency response has a relatively shallow drop-off.

The phase shift depends on the frequency, causing signals composed of multiple frequencies to be distorted.

Angular frequency \(\omega_c=|p|\), is is known as the *cutoff, break, -3dB* or *half-power* frequency because the magnitude of the transfer function \(\eqref{eq:polar2}\) equals \(1/\sqrt{2}\)

The attenuation slope is calculated by first expressing the magnitude relative to the cutoff angular frequency \(\omega_c\)

The rate of change of attenuation is usually expressed in dB/decade, where an decade is a factor of 10 in frequency, *\(\omega_2=10\omega_1\)*

This single-pole filter gives has a relatively shallow -20 dB/decade drop-off.

In general, the cutoff frequency is equal to the radial distance of the poles or zeros from the origin of the \(s\)-plane. For information on sketching the Bode magnitude plot from the poles and zeros, refer to Understanding Poles and Zeros [MIT 3.1].

### Nyquist plot

The Nyquist plots display both amplitude and phase angle on a single plot, using the angular frequency as the parameter. It helps visualize if a system is stable or unstable.

Starting with the transfer function \(\eqref{eq:transferpolynomial}\)

Evaluate at \(s=j\omega\) and split into real and imaginary parts

Plot the frequency transfer function for \(-\infty\lt\omega\lt\infty\), indicating an increase of frequency using an arrow. A dashed line is used for negative frequencies. (The plot was generated using the GNU/Octave as shown in the appendix.)

From the plot we see that for \(\omega=0\) the gain is 1, and for \(\omega\to\infty\) the gain becomes 0. The high frequency portion of the plot approaches the origin at an angle of -90 degrees. For more information on Nyquist refer to Determining Stability using the Nyquist Plot [swarthmore].

## Appendix

### Unit Step Response in GNU/Octave

clc; close all; clear all; format short eng R=100; # 100Ohm C=470e-9; # 470nF w=logspace(3,5,200); f=w/(2*pi); t=linspace(0,2e-3,200); p=-1/(R*C); u=1-e.^(p*t); h=plot(t,u); axis([min(t) max(t) 0 2]); #1.75 xlabel('time [s]'); ylabel('|h(t)|'); t=['Step Response(t), C=' num2str(C*1e9) 'nF, R=' num2str(R) '\Omega'] title(t, &amp;amp;amp;amp;amp;amp;amp;amp;quot;fontsize&amp;amp;amp;amp;amp;amp;amp;amp;quot;, 15);

### Frequency Response in GNU/Octave

clc; close all; clear all; format short eng R=100; # 100Ohm C=470e-9; # 470nF f=logspace(1,6,200); w=2*pi*f; p=-1/(R*C); u=-20*log10(sqrt(1+(w*R*C).^2)); h=semilogx(f,u); hold on; wn=-p; plot([wn/(2*pi) wn/(2*pi)], get(gca,'YLim'),'k--'); text(wn/(2*pi),5,'|p|/2\pi'); f1=-p/(2*pi); fmax=max(f); asymp=-20*log10((fmax-f1)/f1); plot([min(f) f1 fmax],[0 0 asymp ],'k--'); hold off poles=[-p -p]; figure(1); grid off; axis([min(f) max(f) -80 40]); xlabel('frequency [Hz]'); ylabel('20log| H(t)|'); leg=[strread(num2str(R,1),'%s');'asymptote']; t=['Bode Magnitude in dB(f), C=' num2str(C*1e9) 'nF, R=' num2str(R) '\Omega']; title(t, &amp;amp;amp;amp;amp;amp;amp;amp;quot;fontsize&amp;amp;amp;amp;amp;amp;amp;amp;quot;, 15); hold off;

### Nyquist Diagram in GNU/Octave

clc; close all; clear all; format short eng pkg load control R=100; # 100Ohm C=470e-9; # 470nF p=-1/(R*C); H = tf([-p], [ 1 -p ]); [mag, phi, w] = bode(H); nyquist(H); h=gcf; axis ([-0.2, 1.2, -.7, .7], &amp;amp;amp;amp;amp;amp;amp;amp;quot;square&amp;amp;amp;amp;amp;amp;amp;amp;quot;);

### Square Wave in GNU/Octave

clf; t=linspace(0,2e-3,1e4); # t from 0 to 2 msec, with 10,000 steps f=2e3; # input frequency [Hz] R=100; # 100 Ohms C=470e-9; # 470 nF M=1e5; # number of harmonics ut=0; # input signal (square wave) yt=0; # output signal w=2*pi*f; # omega fc=1/(2*pi*R*C) # cutoff (-3dB) frequency for n=1:2:M, nwt = n*w*t; nwRC = n*w*R*C; ut = ut + 4/pi * sin(nwt) / n; argH = 1 / sqrt( 1 + nwRC^2 ); angH = -atan(nwRC); yt = yt + 4/pi * argH * sin(nwt + angH) / n; end plot(t*1e3,ut, 'b-',t*1e3,yt) title('Square Wave input to RC filter') xlabel('t [msec]') ylabel('[Volt]') grid on; legend('u(t)','y(t)') saveas(1,&amp;amp;amp;amp;amp;amp;amp;amp;quot;square.svg&amp;amp;amp;amp;amp;amp;amp;amp;quot;)