Discrete FIR low-pass filter (unfinished)

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

$$ y[n]=x[n]+x[n-1] \label{eq:difeq} $$

where $x[n]$ is the input and $y[n]$ is the output at sample $n$.

The signal flow graph for this filter is

FIR Filter Signal Flow — Signal flow for simple FIR Low Pass Filter
($z^{-1}$ means a delay of one sample)

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]$

$$ x[n]=A\mathrm{e}^{j(\omega nT+\phi)}=A\Big(\cos(\omega nT+\phi)+j\sin(\omega nT+\phi)\Big) \label{eq:xn} $$

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.

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

$$ x[n]=\mathrm{e}^{j\omega nT} \label{eq:xn2} $$

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

$$ \begin{align} y[n] &= \mathrm{e}^{j\omega nT}+\mathrm{e}^{j\omega (n-1)T} \nonumber \\ &= \mathrm{e}^{j\omega nT}+\mathrm{e}^{j\omega nT}\mathrm{e}^{-j\omega T} \nonumber \\ &= \underbrace{(1+\mathrm{e}^{-j\omega T})}_{H(\omega)}\ \ \underbrace{\mathrm{e}^{j\omega nT}}_{x[n]} \label{eq:yn} \end{align} $$

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)$.

$$ \shaded{ H(\omega)\triangleq 1+\mathrm{e}^{-j\omega T} } $$

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

$$ \begin{align} H\left(\mathrm{e}^{j\omega T}\right) &=1+\mathrm{e}^{-j\omega T}\nonumber\\ &=\left( \mathrm{e}^{j\omega T/2} + \mathrm{e}^{-j\omega T/2} \right)\mathrm{e}^{-j\omega T/2}\nonumber\\ &=2\cos\left(\tfrac{\omega T}{2}\right)\ \mathrm{e}^{-j\omega T/2} \end{align} $$

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

$$ H(\omega)=|H(\omega)|\,e^{j\varphi(\omega)} $$

Gain and Phase

The gain $G(\omega)$ follows

$$ \begin{align} G(\omega) &=\left|2\cos\left(\tfrac{\omega T}{2}\right)\ \mathrm{e}^{-j\omega T/2}\right|&\overset{|\mathrm{e}^{-j\alpha}|=1}{\implies}\nonumber\\ &=2\left|\cos\left(\tfrac{\omega T}{2}\right)\right|&\overset{\cos\alpha=-\cos\alpha}{\implies}\nonumber\\ &=2\cos\left(\tfrac{\omega T}{2}\right) \end{align} $$

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

$$ \begin{align} \varphi(\omega)&=\angle\Big(2\cos\left(\tfrac{\omega T}{2}\right)\ \mathrm{e}^{-j\omega T/2}\Big)\nonumber\\ &=\angle{2}+\angle{\cos\left(\tfrac{\omega T}{2}\right)}+\angle\mathrm{e}^{-j\omega T/2}\\ &=0 + 0 -\frac{\omega T}{2}\nonumber\\ &=-\frac{\cancel{2}\pi f}{\cancel{2}}\ \frac{1}{f_s}\nonumber\\ &=-\pi\frac{f}{f_s},&|f|\leq\frac{f_s}{2} \end{align} $$

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

$$ \left\{ \shaded{ \begin{align} H(\omega)&=|H(\omega)|\ \mathrm{e}^{-\varphi(\omega)}\nonumber\\ |H(\omega)|&=2\cos\left(\tfrac{\omega T}{2}\right)\nonumber\\ \varphi(\omega)&=-\pi\frac{f}{f_s},&|f|\leq\frac{f_s}{2}\nonumber \end{align} } \right. $$

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}$

$$ \begin{align} y[n] &=\overbrace{|H(\omega)|\ \mathrm{e}^{j\varphi(\omega)}}^{H(\omega)}\ \overbrace{A\mathrm{e}^{j(\omega nT+\phi)}}^{x[n]}\nonumber\\ &=\underbrace{|H(\omega)|A}_\text{amplitude}\ \mathrm{e}^{j(\omega nT+\overbrace{\phi+\varphi(\omega)}^\text{phase})}\nonumber\\ \end{align} $$

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

$$ H(z)=\frac{b_0+b_1z^{-1}+b_2z^{-2}+\ldots}{1 – a_1z^{-1}+a_2z^{-2}+\ldots}=\frac{\sum_{k=0}^\infty b_k\,z^{-k}}{1-\sum_{k=1}^\infty a_k\,z^{-k}} $$

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

$$ H(z)=\frac{bz}{z-a}=\frac{b}{1-az^{-1}} $$

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

$$ y[n]=b\,x[n]+a\,y[n-1] $$

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

$$ y[n]=y[n-1]+b(\,x[n]-y[n-1]\,) $$

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.

$$ \begin{align} d&=\mathrm{e}^{-\frac{1}{\tau}}\ \ \implies \tau=-\frac{1}{\ln d} \end{align} $$

The cutoff frequency $f_c$ is also related to $d$ as

$$ \begin{align} d&=\mathrm{e}^{-2\pi f_c}\ \ \implies f_c=-\frac{\ln d}{2\pi} \end{align} $$

Further look at IIR filters:

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

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.

Schematic RC filter

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.

$$ H(s) = \frac{Y(s)}{U(s)}=\frac{Z_C}{Z_C+Z_R} = \frac{1/sC}{1/sC+R} = \frac{1}{sRC+1} \label{eq:voltagedivider} $$

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

$$ \shaded{ \begin{array}{ccc} H(s) = K\frac{1}{s-p}, &K = \frac{1}{RC}, &p = -\frac{1}{RC} \end{array} } \label{eq:transferpolynomial} $$

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)$.

$$ \begin{align} \gamma(t)&=\begin{cases} 0 & t<0 \\ 1 & t\geq0 \end{cases}\nonumber\\ \implies\ \Gamma(s)&=\mathcal{L}\left\{\gamma(t)\right\}=\frac{1}{s} \end{align} \label{eq:unitstep} $$

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

$$ \begin{align} Y(s)&=\Gamma(s)\cdot H(s)\nonumber \\ &= \frac{1}{s}\cdot K\frac{1}{(s-p)}\nonumber \\ &= K\frac{1}{s(s-p)} \end{align} \label{eq:multiplication} $$

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]

$$ Y(s)=K\frac{1}{s(s+a)}\equiv\frac{c_0}{s}+\frac{c_1}{s+a} \label{eq:heaviside} $$

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

$$ \left\{ \begin{eqnarray} -p\frac{\cancel{s}}{\cancel{s}(s-p)} &\equiv& \frac{\cancel{s}c_0}{\cancel{s}}+\frac{sc_1}{s-p}\nonumber \\ -p\frac{\cancel{s-p}}{s\cancel{(s-p)}} &\equiv& \frac{c_0(s-p)}{s}+\frac{c_1\cancel{(s-p)}}{\cancel{s-p}}\nonumber \end{eqnarray} \right. $$

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$.

$$ \begin{eqnarray} c_0 &=& \left.\frac{-p}{s-p}\right|_{s=0}=1 \label{eq:constants1a} \\ c_1 &=& \left.\frac{-p}{s}\right|_{s=p}=-1 \label{eq:constants1b} \end{eqnarray} $$

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

$$ \begin{align} y(t)&= \mathcal{L}^{-1}\left\{\frac{c_0}{s}\right\}+\mathcal{L}^{-1}\left\{\frac{c_1}{s-p}\right\}, & t\geq0\nonumber \\ &=c_0+c_1e^{pt}, & t\geq0 \end{align} $$

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

$$ \shaded{ \begin{array}{ccr} y(t)=1-e^{pt}, & p=-\frac{1}{RC}, & t\geq0 \end{array} } $$

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

$$ y_{ss}(t)=|H(j\omega)|\,\sin(\omega t+\angle H(j\omega))\,\gamma(t) $$

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

$$ H(s)=K\frac{1}{s-p},\ K = \frac{1}{RC},\ p = -\frac{1}{RC} \label{eq:splane} $$

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

$$ \begin{align} \lim_{s \implies j0} |H(s)| &=-p\frac{1}{0-p} =1 \\ \lim_{s \implies j\infty} |H(s)| &=-p\frac{1}{\infty-p} =0 \end{align} $$

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

$$ \left\{ \begin{align} H(s) &= |H(s)|\ e^{j\angle{H(s)}}\nonumber \\ |H(s)| &= K \frac{\prod_{i=1}^m\left|(s-z_i)\right|}{\prod_{i=1}^n\left|(s-p_i)\right|} =K \frac{1}{\left|s-p\right|}\nonumber \\ \angle{H(s)}&=\sum_{i=1}^m\angle(s-z_i)-\sum_{i=1}^n\angle(s-p_i) =-\angle(s-p)\nonumber \end{align} \right. $$

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)$.

$$ \left\{ \begin{align} |H(j\omega)| &=K \frac{1}{\left|j\omega-p\right|}= K \frac{1}{\sqrt{\omega^2+p^2}}\nonumber \\ \angle{H(j\omega)}&=-\angle(j\omega-p)=\mathrm{atan2}(\omega,-p)\nonumber\\ &=-\arctan\frac{\omega}{-p},\ p\lt 0\land p\in\mathbb{R}\nonumber \end{align} \right. \label{eq:polar1} $$

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

$$ \left\{ \begin{eqnarray} |H(j\omega)| &=&\frac{1}{RC} \frac{1}{\sqrt{\omega^2+\left(\frac{1}{RC}\right)^2}}\nonumber \\ \angle{H(j\omega)}&=&-\arctan\frac{\omega}{\frac{1}{RC}}\nonumber \end{eqnarray} \right. \label{eq:eq101} $$

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.

$$ x(t)=\frac{4}{\pi}\sum_{n=1,3,\dots }^\infty \frac{\sin\left(n\omega t\right) }{n} $$

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

$$ \left\{ \begin{align} y_{ss}(t)&=\frac{4}{\pi}\sum_{n=1,3,\dots }^\infty |H(jn\omega)| \frac{\sin\left(n\omega t+\angle H(jn\omega)\right) }{n}\nonumber\\ |H(jn\omega)| &=\frac{1}{\sqrt{(1+n\omega RC)^2}}\nonumber\\ \angle{H(jn\omega)}&=-\arctan(n\omega RC)\nonumber\\ RC&=47\cdot 10^{-6}\nonumber\\ \omega&=2\pi20\cdot 10^{3}\nonumber \end{align} \right. $$

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.

$$ \begin{align} |H(j\omega)|\ &= \frac{1}{\sqrt{1+(\omega RC)^2}} \nonumber\\ |H_{dB}(j\omega)|\ &= -20\log\sqrt{1+ (\omega RC)^2} \nonumber\\ \angle{H(j\omega)}\ &=-\arctan\left(\omega RC\right) \label{eq:polar2} \end{align} $$

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}$

$$ \begin{align} |H(j\omega_c)|\ &=\dfrac{1}{\sqrt{1+(\omega_c RC)^2}}=\frac{1}{\sqrt{2}},\ \omega_c=\left|p\right|=\frac{1}{RC} \\ |H_{dB}(j\omega_c)|\ &= 20\log\frac{1}{\sqrt{2}} \approx-3\rm{\ dB} \end{align} $$

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

$$ \begin{array}{llr} &|H(j\omega)| = \frac{1}{\sqrt{1+\left(\frac{\omega}{\omega_c}\right)^2}}\\ \implies&|H(j\omega)| =\frac{\omega}{\omega_c} & \forall\ {\omega\gg\omega_c} \end{array} $$

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$

$$ \begin{gather} \begin{aligned} |H_{dB}(j\omega_2)|-|H_{dB}(\omega_1)|\ &= -20\log\left(\frac{10\omega_1}{\omega_c}\right)+20\log\left(\frac{\omega_1}{\omega_c}\right) \\ &=20\log\left(\frac{\omega_1}{\omega_c}\frac{\omega_c}{10\omega_1}\right) \\ &=-20 \mathrm{\ dB/decade} \end{aligned} \end{gather} $$

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}$

$$ \begin{gather} \begin{array}{ccc} H(s) = K\frac{1}{s-p}, &K = \frac{1}{RC}, &p = -\frac{1}{RC} \end{array} \end{gather} $$

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

$$ \begin{gather} \begin{aligned} H(j\omega)\ &=K\frac{1}{j\omega-p}\\ &=K\frac{1}{j\omega-p}\frac{j\omega+p}{j\omega+p}\\ &=K\frac{j\omega+p}{-\omega^2-p^2} \\ &=K\frac{-j\omega-p}{\omega^2+p^2} \end{aligned}\\ \implies \left\{ \begin{aligned} \Re\left\{{H(j\omega)}\right\}=&K\frac{-p}{p^2+\omega^2}\\ \Im\left\{{H(j\omega)}\right\}=&K\frac{-\omega}{p^2+\omega^2}\\ \end{aligned} \right. \end{gather} $$

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;amp;quot;fontsize&amp;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;amp;quot;fontsize&amp;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;amp;quot;square&amp;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;amp;quot;square.svg&amp;amp;amp;amp;amp;amp;amp;amp;amp;quot;)

Cookie	Duration	Description
cookielawinfo-checkbox-analytics	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Analytics".
cookielawinfo-checkbox-functional	11 months	The cookie is set by GDPR cookie consent to record the user consent for the cookies in the category "Functional".
cookielawinfo-checkbox-necessary	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookies is used to store the user consent for the cookies in the category "Necessary".
cookielawinfo-checkbox-others	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Other.
cookielawinfo-checkbox-performance	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Performance".
viewed_cookie_policy	11 months	The cookie is set by the GDPR Cookie Consent plugin and is used to store whether or not user has consented to the use of cookies. It does not store any personal data.