Category: Tips & Tricks

Design FIR filter to reject unwanted frequencies

Let’s see how to design a simple digital FIR filter to reject unwanted frequencies in an incoming signal. As a per-requisite, I urge you to read through this post: Introduction to digital filter design

Background on transfer function

The transfer function of a system provides the underlying support for ascertaining vital system response characteristics without solving the difference equations of the system. As mentioned earlier, the transfer function of a generic filter in Z-domain is given by ratio of polynomials in z

\[H(z) = \frac{ \displaystyle{\sum_{i=0}^{M} b_k z^{-1}}}{ 1 + \displaystyle{\sum_{i=1}^{N} a_k z^{-1}} } \quad\quad (1)\]

The values of z when H(z) =0 are called zeros of H(z). The values of z when H(z) = ∞ are called poles of H(z).

It is often easy to solve for zeros {zi} and poles {pj}, when the polynomials in the numerator and denominator are expressed as resolvable factors.

\[H(z) = \frac{ \displaystyle{\sum_{i=0}^{M} b_k z^{-1}}}{ 1 + \displaystyle{\sum_{i=1}^{N} a_k z^{-1}} } = \frac{N(z)}{D(z)} = \frac{b_M}{a_N} \frac{(z – z_1)(z – z_2)\cdots (z – z_M)}{(z – p_1)(z – p_2)\cdots (z – p_N)} \quad\quad (2) \]

The zeros {zi} are obtained by finding the roots of the equation

\[N(z) = 0 \quad\quad (3)\]

The poles {pj} are obtained by finding the roots of the equation

\[D(z) = 0 \quad\quad (4) \]

Pole-zero plots are suited for visualizing the relationship between the Z-domain and the frequency response characteristics. As mentioned before, the frequency response of the system H(e jω) can be computed by evaluating the transfer function H(z) at specific values of z = e jω. Because the frequency response is periodic with period 2π, it is sufficient to evaluate the frequency response for the range -π <= ω < π (that is one loop around the unit circle on the z-plane starting from z=-1 and ending at the same point).

FIR filter design

FIR filters contain only zeros and no poles in the pole-zero plot (in fact all the poles sit at the origin for a causal FIR). For an FIR filter, the location of zeros of H(z) on the unit-circle nullifies specific frequencies. So, to design an FIR filter to nullify specific frequency ω, we just have to place the zeros at corresponding locations on the unit circle (z=ejω) where the gain of the filter needs to be 0. Let’s illustrate this using an example.

For this illustration, we would use this audio file as an input to the filtering system. As a first step in the filter design process, we should understand the nature of the input signal. Discrete-time Fourier transform (DTFT) is a tool for analyzing the frequency domain characteristics of the given signal.

The following function is used to compute the DTFT of the sequence read from the audio file.

import numpy as np
import matplotlib.pyplot as plt

from math import ceil,log,pi,cos
from scipy.fftpack import fft,fftfreq,fftshift

def compute_DTFT(x,M=0):
    """
    Compute DTFT of the given sequence x
    M is the desired length for computing DTFT (optional).
    Returns the DTFT X[k] and corresponding frequencies w (omega) arranged as -pi to pi
    """
    N = max(M,len(x))
    N = 2**(ceil(log(N)/log(2)))
    
    X = fftshift(fft(x,N))
    w = 2*pi*fftshift(fftfreq(N))    
    return (X,w)

Let’s read the audio file, load the samples as a signal sequence , and plot the sequence in time-domain/frequency domain (using DTFT).

from scipy.io.wavfile import read

samplerate, x = read('speechWithNoise.wav')
duration = len(x)/samplerate
time = np.arange(0,duration,1/samplerate)

fig1, (ax1,ax2) = plt.subplots(nrows=2,ncols=1)
ax1.plot(time,x)
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Amplitude')
ax1.set_title('speechWithNoise.wav - x[n]')

(X,w)= compute_DTFT(x)
ax2.plot(w,abs(X))
ax2.set_xlabel('Normalized frequency (radians/sample)')
ax2.set_ylabel('|X[k]|')
ax2.set_title('Magnitude vs. Frequency')
Time-domain and frequency domain characteristics of the given audio sample
Figure 1: Time-domain and frequency domain characteristics of the given audio sample

The magnitude vs. frequency plot simply shows huge spikes at θ=±1.32344 radians. The location of the spikes are captured by using the numpy.argmax↗ function.

maxIndex = np.argmax(X)
theta = w[maxIndex]
print(theta)

Since a sinusoid can be mathematically represented as

\[x[n] = cos (\theta n) = \frac{1}{2}\left( e^{j \theta n} + e^{-j \theta n }\right) \quad\quad (5)\]

The two spikes at θ=±1.32344 radians in the frequency domain, will manifest as a sinusoidal signal in the time domain.

Zooming in the area between θ= ±0.4 radians, the frequency domain plot reveals a hidden signal.

Figure 2: Hidden signal revealed in frequency domain

Now, our goal is to design an FIR filter that should reject the sinusoid at θ=±1.32344 radians, so that only the hidden signal gets filtered in.

Since the sinusoid that we want to reject is occurring at some θ radians in the frequency domain, for the FIR filter design, we place two zeros at

\[z_1 = e^{j \theta} \quad\quad z_2 = e^{-j \theta}\quad\quad (6)\]

Therefore, the transfer function of the filter system is given by

\[\begin{aligned} H_f(z) &= \left( 1 – z_1 z^{-1}\right)\left(1 – z_2 z^{-1} \right) \\ &= \left( 1 – e^{j \theta} z^{-1}\right)\left(1 -e^{-j \theta}z^{-1} \right) \\ & = 1 – 2\;cos\left(\theta\right)z^{-1} + z^{-2} \end{aligned}\]

This is a second order FIR filter. The coefficients of the filter are

\[b_0 = 1,\; b_1 = – 2 cos (\theta),\; b_2 = 1 \text{ and } a_0=1\]

For the given problem, to should reject the sinusoid at θ=±1.32344 radians, we set θ=1.32344 in the filter coefficients above.

Filter in action

Filter the input audio signal through the designed filter and plot the filtered output in time-domain and frequency domain. The lfilter function from scipy.signal package↗ is utilized for the filtering operation.

from scipy.signal import lfilter
y_signal = lfilter(b, a, x)
fig3, (ax3,ax4) = plt.subplots(nrows=1,ncols=2)
ax3.plot(time,y_signal,'g')
ax3.set(title='Noise removed speech - y[n]',xlabel='Time (s)',ylabel='Amplitude')

(Y,w)= compute_DTFT(y_signal)
ax4.plot(w,abs(Y)/max(abs(Y)),'r')
ax4.set(title='Frequency content of Y',xlabel='Normalized frequency (radians/sample)',ylabel='Magnitude - |Y[k]|')

The filter has effectively removed the sinusoidal noise, as evident from both time-domain and frequency domain plots.

Figure 4: Extracted speech and its frequency content

Save the filtered output signal as .wav file for audio playback

from scipy.io.wavfile import write
output_data = np.asarray(y_signal, dtype=np.int16)#convert y to int16 format
write("noise_removed_output.wav", samplerate, output_data)

Characteristics of the designed filter

Let’s compute the double sided frequency response of the designed FIR filter. The frequency response of the designed FIR digital filter is computed using freqz function from the scipy.signal package↗.

from scipy.signal import freqz
b = [1,-2*cos(theta),1] #filter coefficients
a = [1]
w, h = freqz(b,a,whole=True)#frequency response h[e^(jw)]
#whole = True returns output for whole range 0 to 2*pi
#To plot double sided response, use fftshift
w = w - 2*np.pi*(w>=np.pi) #convert to range -pi to pi
w = fftshift(w)
h = fftshift(h)

Plot the magnitude response, phase response, pole-zero plot and the impulse response of the designed filter.

#Plot Magnitude-frequency response
fig2, (ax) = plt.subplots(nrows=2,ncols=2)
ax[0,0].plot(w,abs(h),'b')
ax[0,0].set(title='Magnitude response',xlabel='Frequency [radians/sample]',ylabel='Magnitude [dB]')
ax[0,0].grid();ax[0,0].axis('tight');

#Plot phase response
angles = np.unwrap(np.angle(h))
ax[0,1].plot(w,angles,'r')
ax[0,1].set(title='Phase response',xlabel='Frequency [radians/sample]',ylabel='Angles [radians]')
ax[0,1].grid();ax[0,1].axis('tight');

#Transfer function to pole-zero representation
from scipy.signal import tf2zpk
z, p, k = tf2zpk(b, a)

#Plot pole-zeros on a z-plane
from  matplotlib import patches
patch = patches.Circle((0,0), radius=1, fill=False,
                    color='black', ls='dashed')
ax[1,0].add_patch(patch)
ax[1,0].plot(np.real(z), np.imag(z), 'xb',label='Zeros')
ax[1,0].plot(np.real(p), np.imag(p), 'or',label='Poles')
ax[1,0].legend(loc=2)
ax[1,0].set(title='Pole-Zero Plot',ylabel='Real',xlabel='Imaginary')
ax[1,0].grid()

#Impulse response
#create an impulse signal
imp = np.zeros(20)
imp[0] = 1

from scipy.signal import lfilter
y_imp = lfilter(b, a, imp) #drive the impulse through the filter
ax[1,1].stem(y_imp,linefmt='-.')
ax[1,1].set(title='Impulse response',xlabel='Sample index [n]',ylabel='Amplitude')
ax[1,1].axis('tight')
Figure 3: Magnitude response, phase response, pole-zero plot and impulse response of the designed second order FIR filter

Questions

Use the comment form below to answer the following questions

1) Is the filter that we designed a lowpass, highpass, bandpass or a bandstop filter ?
2) To reject a single sinusoidal signal, why two zeros are needed in the above filter design ?
3) What do you understand from the phase response plotted above ?

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Similar topics

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Digital filter design – Introduction

Key focus: Develop basic understanding of digital filter design. Learn about fundamentals of FIR and IIR filters and the design choices.

Analog filters and digital filters are the two major classification of filters, depending on the type of signal signal they process. An analog filter, processes continuous-time signal analog signals. Whereas, digital filters process sampled, discrete-time signals.

Design of digital filters, involve the use of both frequency domain and time domain techniques.This is because, the filter specifications are often specified in frequency domain and the implementation is done in time-domain in the form of difference equations. Z-transform and discrete-time frequency transform (DTFT) are typical tools used for frequency domain analysis of filters.

Linear Time-Invariant Causal filter

The choice of filter and the design process depends on design specification, application and the performance issues associates with them. In this section, the focus is on the principles used for designing linear time-invariant causal filters.

A linear system obeys the principle of superposition. This means that an arbitrary signal can be represented as the weighted sums of shifted unit impulse functions. In a linear time-invariant filtering system, the filter coefficients do not change with time. Therefore, if the input signal is time-shifted, there will be a corresponding time-shift in the output signal. The term causal implies that the output of the system depends only on the present and past samples of input or output, and not on the future samples. Causality is required for real-time applications.

Solution for LTI causal filtering

The input-output relationship in a linear time-invariant causal filter system (shown above) is mathematically described using the following difference equation

\[ \begin{align} \sum_{k=0}^{N} a_k y[n-k] &= \sum_{i=0}^{M} b_i x[n-i] \\ \Rightarrow y[n] &= \sum_{i=0}^{M} b_i x[n-i] – \sum_{k=1}^{N} a_k y[n-k] \end{align} \quad \quad (1) \]

where, x[n] and y[n] are the input and output (sampled discrete-time) signals of the filtering system, ak and bi are the coefficients of the filter that is programmed to certain values for achieving the given filtering task. The values of M and N determine the number of such coefficients, in other words, the filter has M zeros and N poles.

Linear time-invariant systems are completely characterized by it response to an impulse signal δ[n]. Therefore, the analysis and solution for a given filtering task, is easily achieved by using the impulse response h[n], which is the response of the LTI system to an impulse signal δ[n] .

Therefore, the solution for the recursive equation in [1] is found by exciting the input with an impulse , that is x[n]= δ[n].

\[h[n] = \sum_{i=0}^{M} b_i \delta[n-i] – \sum_{k=1}^{N} a_k h[n-k] \quad \quad (2) \]

Z-transform and DTFT

Z-transform and discrete-time frequency transform (DTFT) are important tools for analyzing difference equations and frequency response of filters. Z-transform converts a discrete-time signal into a complex frequency domain representation.

The impulse response of a discrete-time causal system is analyzed using the unilateral or one-sided Z-transform. The unilateral Z-transform of a discrete-time signal x(n) is given by

where, n is an integer and z is a complex number with magnitude r and complex argument ω in radians.

\[ z = re^{j \omega} = r \cdot \left[ cos(\omega) + j\; sin (\omega) \right] \quad \quad (4)\]
Figure 1: Z-transform and Z-plane

If we let |z| = r = 1, then the equation for Z-transform transforms into discrete-time Fourier transform.

The Z-plane contains all values of z, whereas, the unit circle (defined by |z|=r=1) contains only z=e jω values. Therefore, we can state that Z-transform may exist anywhere on the Z-plane and DTFT exists only on the unit circle. One period of DTFT is equivalent to one loop around the unit circle on the Z-plane.

One of the important properties of Z-transform is with regards to time-shifting of discrete-time samples. A delay of K samples in the time domain is equivalent to multiplication by z-k in the Z-transform domain.

\[x[n-K] \rightleftharpoons z^{-k} X(z) \quad \quad (7)\]

Re-writing equation (1) in Z-domain,

Therefore, the transfer function of the generic digital filter is given by

From equations (5) and (6), the frequency response of the system can be computed by evaluating the transfer function H(z) on the unit circle (i.e, |z| = r =1 → z = ejω)

IIR filter

From implementation standpoint, there are two classes of digital filters: infinite impulse response (IIR) and finite impulse response (FIR) filters.

When ak ≠ 0 as in equation (2), the filter structure is characterized by the presence of feedback elements. Due to the presence of feedback elements, the impulse response of the filter may not become zero beyond certain point in time, but continues indefinitely and hence the name infinite impulse response (IIR) filter.

\[IIR: \quad \quad h[n] = \sum_{i=0}^{M} b_i \delta[n-i] – \sum_{k=1}^{N} a_k h[n-k], \quad a_k \neq 0 \quad \quad (11) \]

Therefore, equation (9) and (10) are essentially the transfer function and the frequency response of an IIR filtering system.

There exist different methods for implementing the filter structure. Examples include, direct form I structure, direct form II structure, lattice structure, transposition, state space representation etc.., Each method has its own advantage and disadvantage. For example, some method may be robust to precision issues, coefficient sensitivity, lesser memory and computation requirement. The methods are chosen depending on the application.

Figure 2 shows the direct form I signal flow graph for the generic IIR filter described by equation (1). The boxes represented by Z-1 are unit delays. This is because, the Z-transform of unit delay is Z-1

Figure 2: Direct form I signal flow graph for an IIR filter

FIR filter

When ak=0 for all ks, there is no feedback in the filter structure, no poles in the pole-zero plot (in fact all the poles sit at the origin for a causal FIR). The equation is no longer recursive. The impulse response of such filter dies out (becomes zero) beyond certain point in time and hence the name finite impulse response (FIR) filter.

Setting x[n] = δ[n] and ak=0 for all k, the impulse response of an FIR filter is given by,

\[FIR: \quad \quad h[n] = \sum_{i=0}^{M} b_i \delta[n-i] \quad \quad (13) \]

Evaluating the z-transform of impulse input x[n] = δ[n] , in Z-domain following mapping holds

From equation (13) and (14), the transfer function of the FIR filter in Z-domain is given by

The frequency response is given by evaluating the transfer function on the unit circle |z|=1.

The direct form I signal flow graph for the FIR filter is shown in Figure 3

Figure 3: Direct form I signal flow graph for FIR filter

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References

[1] Z- transform – MIT open course ware ↗

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Similar topics

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Linear regression using python – demystified

Key focus: Let’s demonstrate basics of univariate linear regression using Python SciPy functions. Train the model and use it for predictions.

Linear regression model

Regression is a framework for fitting models to data. At a fundamental level, a linear regression model assumes linear relationship between input variables () and the output variable (). The input variables are often referred as independent variables, features or predictors. The output is often referred as dependent variable, target, observed variable or response variable.

If there are only one input variable and one output variable in the given dataset, this is the simplest configuration for coming up with a regression model and the regression is termed as univariate regression. Multivariate regression extends the concept to include more than one independent variables and/or dependent variables.

Univariate regression example

Let us start by considering the following example of a fictitious dataset. To begin we construct the fictitious dataset by our selves and use it to understand the problem of linear regression which is a supervised machine learning technique. Let’s consider linear looking randomly generated data samples.

import numpy as np
import matplotlib.pyplot as plt #for plotting

np.random.seed(0) #to generate predictable random numbers

m = 100 #number of samples
x = np.random.rand(m,1) #uniformly distributed random numbers
theta_0 = 50 #intercept
theta_1 = 35 #coefficient
noise_sigma = 3

noise = noise_sigma*np.random.randn(m,1) #gaussian random noise

y = theta_0 + theta_1*x + noise #noise added target
 
plt.ion() #interactive plot on
fig,ax = plt.subplots(nrows=1,ncols=1)
plt.plot(x,y,'.',label='training data')
plt.xlabel(r'Feature $x_1$');plt.ylabel(r'Target $y$')
plt.title('Feature vs. Target')
Figure 1: Simulated data for linear regression problem

In this example, the data samples represent the feature and the corresponding targets . Given this dataset, how can we predict target as a function of ? This is a typical regression problem.

Linear regression

Let be the pair that forms one training example (one point on the plot above). Assuming there are such sample points as training examples, then the set contains all the pairs .

In the univariate linear regression problem, we seek to approximate the target as a linear function of the input , which implies the equation of a straight line (example in Figure 2) as given by

where, is the intercept, is the slope of the straight line that is sought and is always . The approximated target serves as a guideline for prediction. The approximated target is denoted by

Using all the samples from the training set , we wish to find the parameters that well approximates the relationship between the given target samples and the straight line function .

If we represent the variables s, the input samples for and the target samples as matrices, then, equation (1) can be expressed as a dot product between the two sequences

It may seem that the solution for finding is straight forward

However, matrix inversion is not defined for matrices that are not square. Moore-Penrose pseudo inverse generalizes the concept of matrix inversion to a matrix. Denoting the Moore-Penrose pseudo inverse for as , the solution for finding is

For coding in Python, we utilize the scipy.linalg.pinv function to compute Moore-Penrose pseudo inverse and estimate .

xMat = np.c_[ np.ones([len(x),1]), x ] #form x matrix
from scipy.linalg import pinv
theta_estimate = pinv(xMat).dot(y)
print(f'theta_0 estimate: {theta_estimate[0]}')
print(f'theta_1 estimate: {theta_estimate[1]}')

The code results in the following estimates for , which are very close to the values used to generate the random data points for this problem.

>> theta_0 estimate: [50.66645323]
>> theta_1 estimate: [34.81080506]

Now, we know the parameters of our example system, the target predictions for new values of feature can be done as follows

x_new = np.array([[-0.2],[0.5],[1.2] ]) #new unseen inputs
x_newmat = np.c_[ np.ones([len(x_new),1]), x_new ] #form xNew matrix
y_predict  = np.dot(x_newmat,theta_estimate)
>>> y_predict #predicted y values for new inputs for x_1
array([[43.70429222],
       [68.07185576],
       [92.43941931]])

The approximated target as a linear function of feature, is plotted as a straight line.

plt.plot(x_new,y_predict,'-',label='prediction')
plt.text(0.7, 55, r'Intercept $\theta_0$ = %0.2f'%theta_estimate[0])
plt.text(0.7, 50, r'Coefficient $\theta_1$ = %0.2f'%theta_estimate[1])
plt.text(0.5, 45, r'y= $\theta_0+ \theta_1 x_1$ = %0.2f + %0.2f $x_1$'%(theta_estimate[0],theta_estimate[1]))
plt.legend() #plot legend
Figure 2: Linear Regression – training samples and prediction

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References

[1] Boyd and Vandenberghe , “Convex Optimization”, ISBN: 978-0521833783, Cambridge University Press, 1 edition, March 2004.↗

Related topics

[1] Introduction to Signal Processing for Machine Learning
[2] Generating simulated dataset for regression problems - sklearn make_regression
[3] Hands-on: Basics of linear regression

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Generating simulated dataset for regression problems

Key focus: Generating simulated dataset for regression problems using sklearn make_regression function (Python 3) is discussed in this article.

Problem statement

Suppose, a survey is conducted among the employees of a company. In that survey, the salary and the years of experience of the employees are collected. The aim of this data collection is to build a regression model that could predict the salary from the given experience (especially for the values not seen by the model).

If you are developer, you often have no access to survey data. In this scenario, you wish you could simulate the data for building a regression model.

Generating the dataset

To construct a simulated dataset for this scenario, the sklearn.dataset.make_regression↗ function available in the scikit-learn library can be used. The function generates the samples for a random regression problem.

The make_regression↗ function generates samples for inputs (features) and output (target) by applying random linear regression model. The values for generated samples have to be scaled to appropriate range for the given problem.

import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt #for plotting

x, y, coef = datasets.make_regression(n_samples=100,#number of samples
                                      n_features=1,#number of features
                                      n_informative=1,#number of useful features 
                                      noise=10,#bias and standard deviation of the guassian noise
                                      coef=True,#true coefficient used to generated the data
                                      random_state=0) #set for same data points for each run

# Scale feature x (years of experience) to range 0..20
x = np.interp(x, (x.min(), x.max()), (0, 20))

# Scale target y (salary) to range 20000..150000 
y = np.interp(y, (y.min(), y.max()), (20000, 150000))

plt.ion() #interactive plot on
plt.plot(x,y,'.',label='training data')
plt.xlabel('Years of experience');plt.ylabel('Salary $')
plt.title('Experience Vs. Salary')
Figure 1: Simulated dataset for linear regression problem

If you want the data to be presented in pandas dataframe format:

import pandas as pd
df = pd.DataFrame(data={'experience':x.flatten(),'salary':y})
df.head(10)
Figure 2: Generated dataset presented as pandas dataframe

We have successfully completed generating simulated dataset for regression problems in Python3. Let’s move on to build and train a linear regression model using the generated dataset and use it for predictions.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Related topics

[1] Introduction to Signal Processing for Machine Learning
[2] Generating simulated dataset for regression problems - sklearn make_regression
[3] Hands-on: Basics of linear regression

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Plot audio file as time series using Scipy python

Often the most basic step in signal processing of audio files, one would like to visualize an audio sample file as time-series data.

Audio sounds can be thought of as an one-dimensional vector that stores numerical values corresponding to each sample. The time-series plot is a two dimensional plot of those sample values as a function of time.

Python’s SciPy library comes with a collection of modules for reading from and writing data to a variety of file formats. For example, the scipy.io.wavfile module can be used to read from and write to a .wav format file.

For the following demonstration, sample audio files given in this URL are used for the visualization task.

The read function in the scipy.io.wavefile module can be utilized to open the selected wav file. It returns the sample rate and the data samples.

>>> from scipy.io.wavfile import read #import the required function from the module

>>> samplerate, data = read('CantinaBand3.wav')

>>> samplerate #echo samplerate
22050

>>> data #echo data -> note that the data is a single dimensional array
array([   3,    7,    0, ...,  -12, -427, -227], dtype=int16)

Compute the duration and the time vector of the audio sample from the sample rate

>>> duration = len(data)/samplerate
>>> time = np.arange(0,duration,1/samplerate) #time vector

Plot the time-series data using matplotlib package

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> plt.plot(time,data)
>>> plt.xlabel('Time [s]')
>>> plt.ylabel('Amplitude')
>>> plt.title('CantinaBand3.wav')
>>> plt.show()
Figure 1: Time series plot of audio file using Python Scipy

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Plot FFT using Python – FFT of sine wave & cosine wave

Key focus: Learn how to plot FFT of sine wave and cosine wave using Python. Understand FFTshift. Plot one-sided, double-sided and normalized spectrum using FFT.

Introduction

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT).  Often we are confronted with the need to generate simple, standard signals (sine, cosine, Gaussian pulse, squarewave, isolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Python and how to represent them in frequency domain using FFT. If you are inclined towards Matlab programming, visit here.

This article is part of the book Digital Modulations using Python, ISBN: 978-1712321638 available in ebook (PDF) and Paperback (hardcopy) formats

Sine Wave

In order to generate a sine wave, the first step is to fix the frequency f of the sine wave. For example, we wish to generate a sine wave whose minimum and maximum amplitudes are -1V and +1V respectively. Given the frequency of the sinewave, the next step is to determine the sampling rate.

For baseband signals, the sampling is straight forward. By Nyquist Shannon sampling theorem, for faithful reproduction of a continuous signal in discrete domain, one has to sample the signal at a rate higher than at-least twice the maximum frequency contained in the signal (actually, it is twice the one-sided bandwidth occupied by a real signal. For a baseband signal bandwidth ( to ) and maximum frequency in a given band are equivalent).

For Python implementation, let us write a function to generate a sinusoidal signal using the Python’s Numpy library. Numpy is a fundamental library for scientific computations in Python. In order to use the numpy package, it needs to be imported. Here, we are importing the numpy package and renaming it as a shorter alias np.

import numpy as np

Next, we define a function for generating a sine wave signal with the required parameters.

def sine_wave(f,overSampRate,phase,nCyl):
	"""
	Generate sine wave signal with the following parameters
	Parameters:
		f : frequency of sine wave in Hertz
		overSampRate : oversampling rate (integer)
		phase : desired phase shift in radians
		nCyl : number of cycles of sine wave to generate
	Returns:
		(t,g) : time base (t) and the signal g(t) as tuple
	Example:
		f=10; overSampRate=30;
		phase = 1/3*np.pi;nCyl = 5;
		(t,g) = sine_wave(f,overSampRate,phase,nCyl)
	"""
	fs = overSampRate*f # sampling frequency
	t = np.arange(0,nCyl*1/f-1/fs,1/fs) # time base
	g = np.sin(2*np.pi*f*t+phase) # replace with cos if a cosine wave is desired
	return (t,g) # return time base and signal g(t) as tuple

We note that the function sine wave is defined inside a file named signalgen.py. We will add more such similar functions in the same file. The intent is to hold all the related signal generation functions, in a single file. This approach can be extended to object oriented programming. Now that we have defined the sine wave function in signalgen.py, all we need to do is call it with required parameters and plot the output.

"""
Simulate a sinusoidal signal with given sampling rate
"""
import numpy as np
import matplotlib.pyplot as plt # library for plotting
from signalgen import sine_wave # import the function

f = 10 #frequency = 10 Hz
overSampRate = 30 #oversammpling rate
fs = f*overSampRate #sampling frequency
phase = 1/3*np.pi #phase shift in radians
nCyl = 5 # desired number of cycles of the sine wave

(t,x) = sine_wave(f,overSampRate,phase,nCyl) #function call

plt.plot(t,x) # plot using pyplot library from matplotlib package
plt.title('Sine wave f='+str(f)+' Hz') # plot title
plt.xlabel('Time (s)') # x-axis label
plt.ylabel('Amplitude') # y-axis label
plt.show() # display the figure

Python is an interpreter based software language that processes everything in digital. In order to obtain a smooth sine wave, the sampling rate must be far higher than the prescribed minimum required sampling rate, that is at least twice the frequency – as per Nyquist-Shannon theorem. Hence, we need to sample the input signal at a rate significantly higher than what the Nyquist criterion dictates. Higher oversampling rate requires more memory for signal storage. It is advisable to keep the oversampling factor to an acceptable value.

An oversampling factor of is chosen in the previous function. This is to plot a smooth continuous like sine wave. Thus, the sampling rate becomes . If a phase shift is desired for the sine wave, specify it too.

Figure 1: A 10Hz sinusoidal wave with 5 cycles and phase shift 1/3π radians

Different representations of FFT:

Since FFT is just a numeric computation of -point DFT, there are many ways to plot the result. The FFT, implemented in Scipy.fftpack package, is an algorithm published in 1965 by J.W.Cooley and
J.W.Tuckey for efficiently calculating the DFT.

The SciPy functions that implement the FFT and IFFT can be invoked as follows

from scipy.fftpack import fft, ifft
X = fft(x,N) #compute X[k]
x = ifft(X,N) #compute x[n]

1. Plotting raw values of DFT:

The x-axis runs from to – representing sample values. Since the DFT values are complex, the magnitude of the DFT is plotted on the y-axis. From this plot we cannot identify the frequency of the sinusoid that was generated.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

NFFT=1024 #NFFT-point DFT      
X=fft(x,NFFT) #compute DFT using FFT    

fig1, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
nVals = np.arange(start = 0,stop = NFFT) # raw index for FFT plot
ax.plot(nVals,np.abs(X))      
ax.set_title('Double Sided FFT - without FFTShift')
ax.set_xlabel('Sample points (N-point DFT)')        
ax.set_ylabel('DFT Values')
fig1.show()
Figure 2: Double sided FFT – without FFTShift

2. FFT plot – plotting raw values against normalized frequency axis:

In the next version of plot, the frequency axis (x-axis) is normalized to unity. Just divide the sample index on the x-axis by the length of the FFT. This normalizes the x-axis with respect to the sampling rate . Still, we cannot figure out the frequency of the sinusoid from the plot.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

NFFT=1024 #NFFT-point DFT  
X=fft(x,NFFT) #compute DFT using FFT     

fig2, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
   
nVals=np.arange(start = 0,stop = NFFT)/NFFT #Normalized DFT Sample points         
ax.plot(nVals,np.abs(X))     
ax.set_title('Double Sided FFT - without FFTShift')        
ax.set_xlabel('Normalized Frequency')
ax.set_ylabel('DFT Values')
fig2.show()
Figure 3: Double sided FFT with normalized x-axis (0 to 1)

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

As you know, in the frequency domain, the values take up both positive and negative frequency axis. In order to plot the DFT values on a frequency axis with both positive and negative values, the DFT value at sample index has to be centered at the middle of the array. This is done by using FFTshift function in Scipy Python. The x-axis runs from to where the end points are the normalized ‘folding frequencies’ with respect to the sampling rate .

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,fftshift

NFFT=1024 #NFFT-point DFT      
X=fftshift(fft(x,NFFT)) #compute DFT using FFT  

fig3, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
    
fVals=np.arange(start = -NFFT/2,stop = NFFT/2)/NFFT #DFT Sample points        
ax.plot(fVals,np.abs(X))
ax.set_title('Double Sided FFT - with FFTShift')
ax.set_xlabel('Normalized Frequency')
ax.set_ylabel('DFT Values');
ax.autoscale(enable=True, axis='x', tight=True)
ax.set_xticks(np.arange(-0.5, 0.5+0.1,0.1))
fig.show()
Figure 4: Double sided FFT with normalized x-axis (-0.5 to 0.5)

4. FFT plot – Absolute frequency on the x-axis vs. magnitude on y-axis:

Here, the normalized frequency axis is just multiplied by the sampling rate. From the plot below we can ascertain that the absolute value of FFT peaks at and . Thus the frequency of the generated sinusoid is . The small side-lobes next to the peak values at and are due to spectral leakage.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,fftshift

NFFT=1024     
X=fftshift(fft(x,NFFT))

fig4, ax = plt.subplots(nrows=1, ncols=1) #create figure handle

fVals=np.arange(start = -NFFT/2,stop = NFFT/2)*fs/NFFT
ax.plot(fVals,np.abs(X),'b')
ax.set_title('Double Sided FFT - with FFTShift')
ax.set_xlabel('Frequency (Hz)')         
ax.set_ylabel('|DFT Values|')
ax.set_xlim(-50,50)
ax.set_xticks(np.arange(-50, 50+10,10))
fig4.show()
Figure 5: Double sided FFT – Absolute frequency on the x-axis vs. magnitude on y-axis

5. Power Spectrum – Absolute frequency on the x-axis vs. power on y-axis:

The following is the most important representation of FFT. It plots the power of each frequency component on the y-axis and the frequency on the x-axis. The power can be plotted in linear scale or in log scale. The power of each frequency component is calculated as

Where is the frequency domain representation of the signal . In Python, the power has to be calculated with proper scaling terms.

Figure 6: Power spectral density using FFT

Plotting the PSD plot with y-axis on log scale, produces the most encountered type of PSD plot in signal processing.

Figure 7: Power spectral density (y-axis on log scale) using FFT

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Phase demodulation via Hilbert transform: Hands-on

Key focus: Demodulation of phase modulated signal by extracting instantaneous phase can be done using Hilbert transform. Hands-on demo in Python & Matlab.

Phase modulated signal:

The concept of instantaneous amplitude/phase/frequency are fundamental to information communication and appears in many signal processing application. We know that a monochromatic signal of form x(t) = a cos(ω t + ɸ) cannot carry any information. To carry information, the signal need to be modulated. Different types of modulations can be performed – amplitude modulation, phase modulation / frequency modulation.

In amplitude modulation, the information is encoded as variations in the amplitude of a carrier signal. Demodulation of an amplitude modulated signal, involves extraction of envelope of the modulated signal. This was discussed and demonstrated here.

This article is part of the following books
● Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
● Digital Modulations using Python ISBN: 978-1712321638
● Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

In phase modulation, the information is encoded as variations in the phase of the carrier signal. In its generic form, a phase modulated signal is expressed as an information-bearing sinusoidal signal modulating another sinusoidal carrier signal

\[x(t) = A cos \left[ 2 \pi f_c t + \beta + \alpha sin \left( 2 \pi f_m t + \theta \right) \right]   \quad \quad \quad (1)\]

where, m(t) = α sin (2 π fm t + θ ) represents the information-bearing modulating signal, with the following parameters

α – amplitude of the modulating sinusoidal signal
fm – frequency of the modulating sinusoidal signal
θ – phase offset of the modulating sinusoidal signal

The carrier signal has the following parameters

A – amplitude of the carrier
fc – frequency of the carrier and fc>>fm
β – phase offset of the carrier

Demodulating a phase modulated signal:

The phase modulated signal shown in equation (1), can be simply expressed as

\[x(t) = A cos \left[ \phi(t) \right]    \quad\quad\quad (2) \]

Here,  ɸ(t) is the instantaneous phase  that varies according to the information signal m(t).

A phase modulated signal of form x(t) can be demodulated by forming an analytic signal by applying Hilbert transform and then extracting the instantaneous phase. This method is explained here.

We note that the instantaneous phase is ɸ(t) = 2 π fc t + β + α sin (2 π fm t + θ) is linear in time, that is proportional to 2 π fc t . This linear offset needs to be subtracted from the instantaneous phase to obtained the information bearing modulated signal. If the carrier frequency is known at the receiver, this can be done easily. If not, the carrier frequency term 2 π fc t needs to be estimated using a linear fit of the unwrapped instantaneous phase. The following Matlab and Python codes demonstrate all these methods.

Matlab code

%Demonstrate simple Phase Demodulation using Hilbert transform
clearvars; clc;
fc = 240; %carrier frequency
fm = 10; %frequency of modulating signal
alpha = 1; %amplitude of modulating signal
theta = pi/4; %phase offset of modulating signal
beta = pi/5; %constant carrier phase offset 
receiverKnowsCarrier= 'False'; %If receiver knows the carrier frequency & phase offset

fs = 8*fc; %sampling frequency
duration = 0.5; %duration of the signal
t = 0:1/fs:1-1/fs; %time base

%Phase Modulation
m_t = alpha*sin(2*pi*fm*t + theta); %modulating signal
x = cos(2*pi*fc*t + beta + m_t ); %modulated signal

figure(); subplot(2,1,1)
plot(t,m_t) %plot modulating signal
title('Modulating signal'); xlabel('t'); ylabel('m(t)')

subplot(2,1,2)
plot(t,x) %plot modulated signal
title('Modulated signal'); xlabel('t');ylabel('x(t)')

%Add AWGN noise to the transmitted signal
nMean = 0; %noise mean
nSigma = 0.1; %noise sigma
n = nMean + nSigma*randn(size(t)); %awgn noise
r = x + n;  %noisy received signal

%Demodulation of the noisy Phase Modulated signal
z= hilbert(r); %form the analytical signal from the received vector
inst_phase = unwrap(angle(z)); %instaneous phase

%If receiver don't know the carrier, estimate the subtraction term
if strcmpi(receiverKnowsCarrier,'True')
    offsetTerm = 2*pi*fc*t+beta; %if carrier frequency & phase offset is known
else
    p = polyfit(t,inst_phase,1); %linearly fit the instaneous phase
    estimated = polyval(p,t); %re-evaluate the offset term using the fitted values
    offsetTerm = estimated;
end
    
demodulated = inst_phase - offsetTerm;

figure()
plot(t,demodulated); %demodulated signal
title('Demodulated signal'); xlabel('n'); ylabel('\hat{m(t)}');

Python code

import numpy as np
from scipy.signal import hilbert
import matplotlib.pyplot as plt
PI = np.pi

fc = 240 #carrier frequency
fm = 10 #frequency of modulating signal
alpha = 1 #amplitude of modulating signal
theta = PI/4 #phase offset of modulating signal
beta = PI/5 #constant carrier phase offset 
receiverKnowsCarrier= False; #If receiver knows the carrier frequency & phase offset

fs = 8*fc #sampling frequency
duration = 0.5 #duration of the signal
t = np.arange(int(fs*duration)) / fs #time base

#Phase Modulation
m_t = alpha*np.sin(2*PI*fm*t + theta) #modulating signal
x = np.cos(2*PI*fc*t + beta + m_t ) #modulated signal

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,m_t) #plot modulating signal
plt.title('Modulating signal')
plt.xlabel('t')
plt.ylabel('m(t)')
plt.subplot(2,1,2)
plt.plot(t,x) #plot modulated signal
plt.title('Modulated signal')
plt.xlabel('t')
plt.ylabel('x(t)')

#Add AWGN noise to the transmitted signal
nMean = 0 #noise mean
nSigma = 0.1 #noise sigma
n = np.random.normal(nMean, nSigma, len(t))
r = x + n  #noisy received signal

#Demodulation of the noisy Phase Modulated signal
z= hilbert(r) #form the analytical signal from the received vector
inst_phase = np.unwrap(np.angle(z))#instaneous phase

#If receiver don't know the carrier, estimate the subtraction term
if receiverKnowsCarrier:
    offsetTerm = 2*PI*fc*t+beta; #if carrier frequency & phase offset is known
else:
    p = np.poly1d(np.polyfit(t,inst_phase,1)) #linearly fit the instaneous phase
    estimated = p(t) #re-evaluate the offset term using the fitted values
    offsetTerm = estimated
                           
demodulated = inst_phase - offsetTerm 

plt.figure()
plt.plot(t,demodulated) #demodulated signal
plt.title('Demodulated signal')
plt.xlabel('n')
plt.ylabel('\hat{m(t)}')

Results

Figure 1: Phase modulation – modulating signal and modulated (transmitted) signal
Figure 2: Demodulated signal from the noisy received signal

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

For further reading

[1] V. Cizek, “Discrete Hilbert transform”, IEEE Transactions on Audio and Electroacoustics, Volume: 18 , Issue: 4 , December 1970.↗

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Extract envelope, phase using Hilbert transform: Demo

Key focus: Learn how to use Hilbert transform to extract envelope, instantaneous phase and frequency from a modulated signal. Hands-on demo using Python & Matlab.

If you would like to brush-up the basics on analytic signal and how it related to Hilbert transform, you may visit article: Understanding Analytic Signal and Hilbert Transform

Introduction

The concept of instantaneous amplitude/phase/frequency are fundamental to information communication and appears in many signal processing application. We know that a monochromatic signal of form cannot carry any information. To carry information, the signal need to be modulated. Take for example the case of amplitude modulation, in which a positive real-valued signal modulates a carrier . That is, the amplitude modulation is effected by multiplying the information bearing signal with the carrier signal .

Here, is the angular frequency of the signal measured in radians/sec and is related to the temporal frequency as . The term is also called instantaneous amplitude.

Similarly, in the case of phase or frequency modulations, the concept of instantaneous phase or instantaneous frequency is required for describing the modulated signal.

Here, is the constant amplitude factor (no change in the envelope of the signal) and is the instantaneous phase which varies according to the information. The instantaneous angular frequency is expressed as the derivative of instantaneous phase.

This article is part of the following books
● Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
● Digital Modulations using Python ISBN: 978-1712321638
● Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Definition

Generalizing these concepts, if a signal is expressed as

  • The instantaneous amplitude or the envelope of the signal is given by
  • The instantaneous phase is given by 
  • The instantaneous angular frequency is derived as
  • The instantaneous temporal frequency is derived as

Problem statement

An amplitude modulated signal is formed by multiplying a sinusoidal information and a linear frequency chirp. The information content is expressed as and the linear frequency chirp is made to vary from to . Given the modulated signal, extract the instantaneous amplitude (envelope), instantaneous phase and the instantaneous frequency.

Solution

We note that the modulated signal is a real-valued signal. We also take note of the fact that amplitude/phase and frequency can be easily computed if the signal is expressed in complex form. Which transform should we use such that the we can convert a real signal to the complex plane without altering the required properties ?? Answer: Apply Hilbert transform and form the analytic signal on the complex plane.  Figure 1 illustrates this concept.

Figure 1: Converting a real-valued signal to complex plane using Hilbert Transform

If we express the real-valued modulated signal as an analytic signal, it is expressed in complex plane as

where, (HT[.]) represents the Hilbert Transform operation. Now, the required parameters are very easy to obtain.

The instantaneous amplitude (envelope extraction) is computed in the complex plane as

The instantaneous phase is computed in the complex plane as

The instantaneous temporal frequency is computed in the complex plane as

Once we know the instantaneous phase, the carrier can be regenerated as . The regenerated carrier is often referred as Temporal Fine Structure (TFS) in Acoustic signal processing [1].

Matlab

fs = 600; %sampling frequency in Hz
t = 0:1/fs:1-1/fs; %time base
a_t = 1.0 + 0.7 * sin(2.0*pi*3.0*t) ; %information signal
c_t = chirp(t,20,t(end),80); %chirp carrier
x = a_t .* c_t; %modulated signal

subplot(2,1,1); plot(x);hold on; %plot the modulated signal

z = hilbert(x); %form the analytical signal
inst_amplitude = abs(z); %envelope extraction
inst_phase = unwrap(angle(z));%inst phase
inst_freq = diff(inst_phase)/(2*pi)*fs;%inst frequency

%Regenerate the carrier from the instantaneous phase
regenerated_carrier = cos(inst_phase);

plot(inst_amplitude,'r'); %overlay the extracted envelope
title('Modulated signal and extracted envelope'); xlabel('n'); ylabel('x(t) and |z(t)|');
subplot(2,1,2); plot(cos(inst_phase));
title('Extracted carrier or TFS'); xlabel('n'); ylabel('cos[\omega(t)]');

Python

import numpy as np
from scipy.signal import hilbert, chirp
import matplotlib.pyplot as plt

fs = 600.0 #sampling frequency
duration = 1.0 #duration of the signal
t = np.arange(int(fs*duration)) / fs #time base

a_t =  1.0 + 0.7 * np.sin(2.0*np.pi*3.0*t)#information signal
c_t = chirp(t, 20.0, t[-1], 80) #chirp carrier
x = a_t * c_t #modulated signal

plt.subplot(2,1,1)
plt.plot(x) #plot the modulated signal

z= hilbert(x) #form the analytical signal
inst_amplitude = np.abs(z) #envelope extraction
inst_phase = np.unwrap(np.angle(z))#inst phase
inst_freq = np.diff(inst_phase)/(2*np.pi)*fs #inst frequency

#Regenerate the carrier from the instantaneous phase
regenerated_carrier = np.cos(inst_phase)

plt.plot(inst_amplitude,'r'); #overlay the extracted envelope
plt.title('Modulated signal and extracted envelope')
plt.xlabel('n')
plt.ylabel('x(t) and |z(t)|')
plt.subplot(2,1,2)
plt.plot(regenerated_carrier)
plt.title('Extracted carrier or TFS')
plt.xlabel('n')
plt.ylabel('cos[\omega(t)]')

Results

Figure 2: Amplitude modulation using a chirp signal and extraction of envelope and TFS

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Reference

[1] Moon, Il Joon, and Sung Hwa Hong. “What Is Temporal Fine Structure and Why Is It Important?” Korean Journal of Audiology 18.1 (2014): 1–7. PMC. Web. 24 Apr. 2017.↗

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Understanding Analytic Signal and Hilbert Transform

Key focus of this article: Understand the relationship between analytic signal, Hilbert transform and FFT. Hands-on demonstration using Python and Matlab.

Introduction

Fourier Transform of a real-valued signal is complex-symmetric. It implies that the content at negative frequencies are redundant with respect to the positive frequencies. In their works, Gabor [1] and Ville [2], aimed to create an analytic signal by removing redundant negative frequency content resulting from the Fourier transform. The analytic signal is complex-valued but its spectrum will be one-sided (only positive frequencies) that preserved the spectral content of the original real-valued signal. Using an analytic signal instead of the original real-valued signal, has proven to be useful in many signal processing applications. For example, in spectral analysis, use of analytic signal in-lieu of the original real-valued signal mitigates estimation biases and eliminates cross-term artifacts due to negative and positive frequency components [3].

This article is part of the following books
● Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
● Digital Modulations using Python ISBN: 978-1712321638
● Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Continuous-time analytic signal

Let be a real-valued non-bandlimited finite energy signal, for which we wish to construct a corresponding analytic signal .  The Continuous Time Fourier Transform of is given by

Lets say the magnitude spectrum of is as shown in Figure 1(a). We note that the signal is a real-valued and its magnitude spectrum is symmetric and extends infinitely in the frequency domain.

Figure 1: (a) Spectrum of continuous signal x(t) and (b) spectrum of analytic signal z(t)

As mentioned in the introduction, an analytic signal can be formed by suppressing the negative frequency contents of the Fourier Transform of the real-valued signal. That is, in frequency domain, the spectral content of the analytic signal is given by

The corresponding spectrum of the resulting analytic signal is shown in Figure 1(b).

Since the spectrum of the analytic signal is one-sided, the analytic signal will be complex valued in the time domain, hence the analytic signal can be represented in terms of real and imaginary components as . Since the spectral content is preserved in an analytic signal, it turns out that the real part of the analytic signal in time domain is essentially the original real-valued signal itself . Then, what takes place of the imaginary part ? Who is the companion to x(t) that occupies the imaginary part in the resulting analytic signal ? Summarizing as equation,

It is interesting to note that Hilbert transform [4] can be used to find a companion function (imaginary part in the equation above) to a real-valued signal such that the real signal can be analytically extended from the real axis to the upper half of the complex plane . Denoting Hilbert transform as , the analytic signal is given by

From these discussion, we can see that an analytic signal for a real-valued signal , can be constructed using two approaches.

●  Frequency domain approach: The one-sided spectrum of is formed from the two-sided spectrum of the real-valued signal by applying equation (2)
● Time domain approach: Using Hilbert transform approach given in equation (4)

One of the important property of an analytic signal is that its real and imaginary components are orthogonal

Discrete-time analytic signal

Since we are in digital era, we are more interested in discrete-time signal processing. Consider a continuous real-valued signal gets sampled at interval seconds and results in real-valued discrete samples , i.e, . The spectrum of the continuous signal is shown in Figure 2(a). The spectrum of that results from the process of periodic sampling is given in Figure 2(b) (Refer here more details on the process of sampling).  The spectrum of discrete-time signal can be obtained by Discrete-Time Fourier Transform (DTFT).

Figure 2: (a) CTFT of continuous signal x(t), (b) Spectrum of x[n] resulted due to periodic sampling and (c) Periodic one-sided spectrum of analytical signal z[n]

At this point, we would like to construct a discrete-time analytic signal from the real-valued sampled signal . We wish the analytic signal is complex valued and should satisfy the following two desired properties

● The real part of the analytic signal should be same as the original real-valued signal.
● The real and imaginary part of the analytic signal should satisfy the following property of orthogonality

In Frequency domain approach for the continuous-time case, we saw that an analytic signal is constructed  by suppressing the negative frequency components from the spectrum of the real signal. We cannot do this for our periodically sampled signal . Periodic mirroring nature of the spectrum prevents one from suppressing the negative components. If we do so, it will vanish the entire spectrum. One solution to this problem is to set the negative half of each spectral period to zero. The resulting spectrum of the analytic signal is shown in Figure 2(c).

Given a record of samples of even length , the procedure to construct the analytic signal is as follows. This method satisfies both the desired properties listed above.

● Compute the -point DTFT of using FFT
● N-point periodic one-sided analytic signal is computed by the following transform

● Finally, the analytic signal (z[n]) is obtained by taking the inverse DTFT of

Matlab

The given procedure can be coded in Matlab using the FFT function. Given a record of real-valued samples , the corresponding analytic signal can be constructed as given next. Note that the Matlab has an inbuilt function to compute the analytic signal. The in-built function is called hilbert.

function z = analytic_signal(x)
%x is a real-valued record of length N, where N is even %returns the analytic signal z[n]
x = x(:); %serialize
N = length(x);
X = fft(x,N);
z = ifft([X(1); 2*X(2:N/2); X(N/2+1); zeros(N/2-1,1)],N);
end

To test this function, we create a 5 seconds record of a real-valued sine signal. The analytic signal is constructed and the orthogonal components are plotted in Figure 3. From the plot, we can see that the real part of the analytic signal is exactly same as the original signal (which is the cosine signal) and the imaginary part of the analytic signal is phase shifted version of the original signal. We note that the imaginary part of the analytic signal is a cosine function with amplitude scaled by which is none other than the Hilbert transform of sine function.

t=0:0.001:0.5-0.001;
x = sin(2*pi*10*t); %real-valued f = 10 Hz
subplot(2,1,1); plot(t,x);%plot the original signal
title('x[n] - original signal'); xlabel('n'); ylabel('x[n]');

z = analytic_signal(x); %construct analytic signal
subplot(2,1,2); plot(t, real(z), 'k'); hold on;
plot(t, imag(z), 'r');
title('Components of Analytic signal'); 
xlabel('n'); ylabel('z_r[n] and z_i[n]');
legend('Real(z[n])','Imag(z[n])');

Python

Equivalent code in Python is given below (tested with Python 3.6.0)

import numpy as np
def main():
    t = np.arange(start=0,stop=0.5,step=0.001)
    x = np.sin(2*np.pi*10*t)
    
    import matplotlib.pyplot as plt
    plt.subplot(2,1,1)
    plt.plot(t,x)
    plt.title('x[n] - original signal')
    plt.xlabel('n')
    plt.ylabel('x[n]')
    
    z = analytic_signal(x)
    
    plt.subplot(2,1,2)
    plt.plot(t,z.real,'k',label='Real(z[n])')
    plt.plot(t,z.imag,'r',label='Imag(z[n])')
    plt.title('Components of Analytic signal')
    plt.xlabel('n')
    plt.ylabel('z_r[n] and z_i[n]')
    plt.legend()

def analytic_signal(x):
    from scipy.fftpack import fft,ifft
    N = len(x)
    X = fft(x,N)
    h = np.zeros(N)
    h[0] = 1
    h[1:N//2] = 2*np.ones(N//2-1)
    h[N//2] = 1
    Z = X*h
    z = ifft(Z,N)
    return z

if __name__ == '__main__':
    main()
Figure 3: Components of analytic signal for a real-valued sine function

Hilbert Transform using FFT

We should note that the hilbert function in Matlab returns the analytic signal $latex z[n]$ not the hilbert transform of the signal . To get the hilbert transform, we should simply get the imaginary part of the analytic signal. Since we have written our own function to compute the analytic signal, getting the hilbert transform of a real-valued signal goes like this.

x_hilbert = imag(analytic_signal(x))

In the coming posts, we will some of the applications of constructing an analytic signal. For example: Find the instantaneous amplitude and phase of a signal, envelope detector for an amplitude modulated signal, detecting phase changes in a sine wave.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References:

[1] D. Gabor, “Theory of communications”, Journal of the Inst. Electr. Eng., vol. 93, pt. 111, pp. 42-57, 1946. See definition of complex signal on p. 432.↗
[2] J. A. Ville, “Theorie et application de la notion du signal analytique”, Cables el Transmission, vol. 2, pp. 61-74, 1948.↗
[3] S. M. Kay, “Maximum entropy spectral estimation using the analytical signal”, IEEE transactions on Acoustics, Speech, and Signal Processing, vol. 26, pp. 467-469, October 1978.↗
[4] Frank R. Kschischang, “The Hilbert Transform”, University of Toronto, October 22, 2006.↗
[5] S. L. Marple, “Computing the discrete-time ‘analytic’ signal via FFT,” Conference Record of the Thirty-First Asilomar Conference on Signals, Systems and Computers , Pacific Grove, CA, USA, 1997, pp. 1322-1325 vol.2.↗

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Choosing FIR or IIR ? Understand design perspective

“What is the best filter that I should use? FIR or IIR ?” is often the question asked by many. There exists two different types of Linear Time Invariant (LTI) filters from transfer function standpoint : FIR (Finite Impulse Response) and IIR (Infinite Impulse Response) filters and myriad design techniques for designing them.  The mere fact that there exists so many techniques for designing a filter, suggests that there is no single optimal filter design.  One has to weigh-in the pros and cons of choosing a filter design by considering the factors discussed here.

Before proceeding to know the secrets of choosing a filter, I urge you to brush up the fundamentals of digital filter design.

Design specification:

A filter design starts with a specification. We may have a simple specification that just calls for removing an unwanted frequency component or it can be a complete design specification that calls for various parameters like – amount of ripples allowed in passband, stop band attenuation, transition width etc.., The design specification usually calls for satisfying one or more of the following:

● Desired Magnitude response –
● Desired Phase response –
● Tolerance specifications – that specifies how much the filter response is allowed to vary when compared with ideal response. Examples include how much ripples allowed in passband, stop band etc..,

Figure 1: FIR and IIR filters can be realized as direct-form structures. (a) Structure with feed-forward elements only – typical for FIR designs. (b) Structure with feed-back paths – generally results in IIR filters. Other structures like lattice implementations are also available.

Given the specification above, the goal of a filter design process to choose the parameters , , and , such that the transfer function of the filter

yields the desired response: . In other words, the design process also involves choosing the number and location of the zeros zeros () and poles () in the pole-zero plot.

Two types of filter can manifest from the given transfer function above.

● When , there is no feedback in the filter structure, no poles in the pole-zero plot (in fact all the poles sit at the origin). The impulse response of such filter dies out (becomes zero) beyond certain point of time and it is classified as Finite Impulse Response (FIR) filter. It provides linear phase characteristic in the passband.
● When , the filter structure is characterized by the presence of feedback elements. Due to the presence of feedback elements, the impulse response of the filter may not become zero beyond certain point, but continues indefinitely and hence the name Infinite Impulse Response (IIR) filter.
● Caution: In most cases, the presence of feedback elements provide infinite impulse response. It is not always true. There are some exceptional cases where the presence of feedback structure may result in finite impulse response. For example, a moving average filter will have a finite impulse response. The output of a moving average filter can be described using a recursive formula, which will result in a structure with feedback elements.

General considerations in design:

As specified earlier, the choice of filter and the design process depends on design specification, application and the performance issues associates with them. However, the following general considerations are applied in practical design.

Minimizing number of computations

In order to minimize memory requirements for storing the filter co-efficients , and to minimize the number of computations, ideally we would like to be as small as possible. For the same specification, IIR filters result in much lower order when compared to its FIR counter part. Therefore, IIR filters are efficient when viewed from this standpoint.

Need for real-time processing

The given application may require processing of input samples in real-time or the input samples may exist in a recorded state (example: video/audio playback, image processing applications, audio compression). From this perspective, we have two types of filter systems

● Causal Filter
— Filter output depends on present and past input samples, not on the future samples. The output may also depend on the past output samples, as in IIR filters. Strictly no future samples.
— Such filters are very much suited for real-time applications.

● Non-Causal Filter
— There are many practical cases where a non-causal filter is required. Typically, such application warrants some form of post-processing, where the entire data stream is already stored in memory.
— In such cases, a filter can be designed that can take in all type of input samples : present, past and future, for processing. These filters are classified as non-causal filters.
— Non-causal filters have much simpler design methods.

It can be often seen in many signal processing texts, that the causal filters are practically realizable. That does not mean non-causal filters are not practically implementable. In fact both types of filters are implementable and you can see them in many systems today. The question you must ask is : whether your application requires real-time processing or processing of pre-recorded samples. If the application requires real-time processing, causal filters must be used. Otherwise, non-causal filters can be used.

Consequences of causal filter:

If the application requires real-time processing, causal filters are the only choice for implementation. Following consequences must be considered if causality is desired.

Ideal filters with finite bands of zero response (example: brick-wall filters), cannot be implemented using causal filter structure. A direct consequence of causal filter is that the response cannot be ideal. Hence, we must design the filter that provides a close approximation to the desired response . If tolerance specification is given, it has to be met.

For example, in the case of designing a low pass filter with given passband frequency () and stopband frequencies (), additional tolerance specifications like allowable passband ripple factor (), stopband ripple factor () need to be considered for the design,. Therefore, the practical filter design involves choosing and and then designing the filter with and that satisfies all the given requirements/responses. Often, iterative procedures may be required to satisfy all the above (example: Parks and McClellan algorithm used for designing optimal causal FIR filters [1]).

Figure 2: A sample filter design specification

For a causal filter, frequency response’s real part and the imaginary part become Hilbert transform pair [2]. Hence the magnitude and phase responses become interdependent.

Stability

For a causal LTI digital filter will be BIBO (Bounded Input Bounded Output) stable, if and only if the impulse response is absolutely summable.

Impulse response of FIR filters are always bounded and hence they are inherently stable. On the other hand, an IIR filter may become unstable if not designed properly.

Consider an IIR filter implemented using a floating point processor that has enough accuracy to represent all the coefficients in the transfer function below

The corresponding impulse response is plotted in Figure 3(a). The plot shows that the impulse response decays rapidly to zero as increases. For this case, the sum in equation (2) will be finite. Hence this IIR filter is stable.

Suppose, if we were to implement the same filter in a fixed point processor and we are forced to round-off the co-efficients to 2 digits after the decimal point, the  same transfer function looks like this

The corresponding impulse response  plotted in Figure 3(b)  shows that the impulse response increases rapidly towards a constant value as  increases. For this case, the sum in equation (2) will tend to infinity. Hence the implemented IIR filter is unstable.

Figure 3: Impact of poorly implemented IIR filter on stability. (a) Stable IIR filter, (b) The same IIR filter becomes unstable due to rounding effects.

Therefore, it is imperative that an IIR filter implementation need to be tested for stability.  To analyze the stability of the filter, the infinite sum in equation (2) need to be computed and it is often difficult to compute this sum. Analysis of pole-zero plot is an alternate solution for this problem. To have a stable causal filter, the poles of the transfer function should lie completely strictly inside the unit circle on the pole-zero plot. The pole-zero plot for the above given transfer functions , are plotted in Figure 4. It shows that for the transfer function , all the poles lie within the unit circle (the region of stability) and hence it is a stable IIR filter. On the other hand,  for the transfer function ,  one poles lie exactly on the unit circle (ie, it is just out of the region of stability) and hence it is an unstable IIR filter.

Figure 4: Impact of poorly implemented IIR filter on stability. (a) Stable IIR filter since all poles lie inside the unit circle, (b) Due to rounding effects, the one pole lies on the unit circle (just outside the region of stability), hence it becomes unstable.

Linear phase requirement

In many signal processing applications, it is needed that a digital filter should not alter the angular relationship between the real and imaginary components of a signal, especially in the passband. In otherwords, the phase relationship between the signal components should be preserved in the filter’s passband.  If not, we have phase distortion.

Phase distortion is a concern in many signal processing applications. For example, in phase modulations like GMSK [3],the entire demodulation process hinges on the phase relationship between the inphase and quadrature components of the incoming signal. If we have a phase distorting filter in the demodulation chain, the entire detection process goes for a toss. Hence, we have to pay attention to the phase characteristics of such filters. To have no phase distortion when processing a signal through a filter, every spectral component of the signal inside the passband should be delayed by the same amount time delay measured in samples. In other words, the phase response in the passband should be a linear function (straight line) of frequency (except for the phase wraps at the band edges).  A filter that satisfies this property is called a linear phase filter. FIR filters provide perfect linear phase characteristic in the passband region (Figure 5) and hence avoids phase distortion. All IIR filters provide non-linear phase characteristic. If a real-time application warrants for zero phase distortion, FIR filters are the immediate choice for design.

It is intuitive to see the phase response of a generalized linear phase filter should follow the relationship , where is the slope and is the intercept when viewing the linear relationship between the frequency and the phase response in the passband (Figure 5). The phase delay and group delay are the important filter characteristics considered for ascertaining the phase distortion and they relate to the intercept and the slope of the phase response in the passband. Linear phase filters are characterized by constant group delay. Any deviation in the group delay from the constant value inside the passband, indicates presence of certain degree of non-linearity in the phase and hence causes phase distortion.

Figure 5: FIR filter showing linear phase characteristics in the passband

Phase delay is the time delay experienced by each spectral component of the input signal. For a filter with the frequency response , the phase delay response defined in terms of phase response as

Group delay is the delay experienced by a group of spectral components within a narrow frequency interval about [4]. The group delay response is defined as the negative derivative of the phase response .

For the generalized linear phase filter, the group delay and phase delay are given by

Summary of design choices

● IIR filters are efficient, they can provide similar magnitude response for fewer coefficients or lower sidelobes for same number of coefficients
● For linear phase requirement, FIR filters are the immediate choice for the design
● FIR filters are inherently stable. IIR filters are susceptible to finite length words effects of fixed point arithmetic and hence the design has to be rigorously tested for stability.
● IIR filters provide less average delay compared to its equivalent FIR counterpart. If the filter has to be used in a feedback path in a system, the amount of filter delay is very critical as it affects the stability of the overall system.
● Given a specification, an IIR design can be easily deduced based on closed-form expressions. However, satisfying the design requirements using an FIR design generally requires iterative procedures.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References:

[1] J. H. McClellan, T. W. Parks, and L. R. Rabiner, “A Computer Program for Designing Optimum FIR Linear Phase Digital Filters,” IEEE Trans, on Audio and Electroacoustics, Vol. AU-21, No. 6, pp. 506-526, December 1973.↗

[2] Frank R. Kschischang, ‘The Hilbert Transform’, Department of Electrical and Computer Engineering, University of Toronto, October 22, 2006.↗

[3]  Thierry Turletti, ‘GMSK in a nutshell’, Telemedia Networks and Systems Group Laboratory for Computer Science, Massachussets Institute of Technology April, 96.↗

[4] Julius O. Smith III, ‘Introduction to digital filters – with audio applications’, Center for Computer Research in Music and Acoustics (CCRMA), Department of Music, Stanford University.↗

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)
Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing