FIR Filter to Match Any Magnitude and Phase Response

• Python
• DFT_FFT
• Tutorials

Note: At the upcoming DSP Online Conference, I’ll be giving a workshop where I’ll dive into practical DSP design techniques like the one described below. This post is a small taste of the kind of insights and hands-on learning we’ll explore together. 

Introduction

Can FIR filters be made to match arbitrary magnitude and phase responses? One of my students asked me this very question, so I created this example "Pi-Filter" as a demonstration of the technique. (Well, my less useful version of a Pi-Filter, not this: Wikipedia: Pi FIlter)

On the left is the magnitude and phase of the frequency response for the filter plotted on a complex plane, and on the right is the same response with the magnitude and phase plotted separately, consistent with how we would typically review frequency responses.  This plot shows both the original target as an arbitrary magnitude and phase response along with the achieved filter as a 165 tap filter. On the complex plot a star marks the start of the frequency response corresponding to DC, and an arrow shows the direction of the response as we sweep frequency along the positive frequency axis.  The frequency response sweeps from DC to the sampling rate, and thus the coefficients of the filter are complex as well to allow for a unique frequency response over this range. (We can also create a filter with real coefficients which would then have a unique response from DC to half the sampling rate, corresponding to 0.5 cycles/sample as a normalized frequency.)

What I will show in this article is the further details on how I created this response, which can be applied to making high quality FIR filters meeting any arbitrary (single-valued) frequency response. The trade space in terms of the accuracy of the match is the total number of coefficients required, and for causal filter implementations, the total delay through the filter. Allowing for that, we can meet any requirement needed to an acceptable error limit (errors which can be 100s of dB lower if we desired!).  

Background for Design Approach

The coefficients of an FIR filter are the samples of the filter's time-domain impulse response. The Fourier Transform of the impulse response is the frequency response. When we first learn about these fundamental filter design and Fourier properties, we may naively conclude that we can simply use the inverse FFT (IFFT) to create the filter coefficients from the target frequency response. This is called the "Frequency Sampling" method of filter design where N samples in frequency results in N filter coefficients in time, and this approach is inferior compared to other filter design techniques. 

A superior and simple filter design method is the "Windowing" method, where samples of the true time domain impulse response for a target frequency response are multiplied with a time-domain window resulting in a finite duration filter length with reduced truncation errors due to tapering the truncated impulse response at its edges. The true time domain impulse response for a target frequency response is given by the inverse Fourier Transform, not the inverse Discrete Fourier Transform (IDFT which the IFFT provides). The IFFT suffers from time domain aliasing which occurs when the true time domain impulse response exceeds the length of the desired filter.  And any sharp transitions in the frequency domain (a response that changes relatively rapidly over a short frequency range), necessitate long impulse response durations. 

For example, a lowpass filter with a "brick-wall" (rectangular pulse shaped) frequency response), requires a Sinc function as the time domain impulse response (a Sinc function in time is the inverse Fourier Transform of a rectangular function in the frequency domain). A Sinc function impulse response extends to plus or minus infinity, meaning we would need an infinitely long filter to realize a true brick-wall frequency response (hence not feasible). The inverse FFT of a sampled rectangular frequency response is the Dirichlet Kernel: the Dirichlet Kernel is an aliased Sinc function and this is a simple demonstration of time domain aliasing that results with the DFT or FFT. The plot below compares the Sinc function to the Dirichlet Kernel: truncating the Sinc as the coefficients of the filter will result in a better approximation to the brick-wall filter than using the Dirichlet Kernel (and likewise with windowing: windowing the Sinc rather than the Dirichlet Kernel would be better). 

To create a low pass filter, we could either use the samples of the Dirichlet Kernel as the impulse response or the samples of the Sinc function. In either case we could multiply the impulse responses with a window to reduce truncation error, but if we used the Dirichlet Kernel (from the IFFT), the filter would have the errors due to time domain aliasing as well as the truncation error (which itself is minimized with the window). It is a much better approach (less error for the same number of coefficients) to use the true time domain impulse response (in this case a Sinc) and then multiply that with a window. 

So that's straightforward for a low pass filter and other shapes where we can easily mathematically compute the inverse Fourier Transform or are otherwise given the target true impulse response needed.  So what do we do when have a target frequency response where we can't easily determine the true impulse response?  The trick is to more densely sample the target frequency response prior to computing the IFFT.  As the number of samples in the target response as an FFT increases, the IFFT approaches the true inverse Fourier Transform. So in this regard we are using the IFFT as a tool to estimate the actual inverse Fourier Transform free of time domain aliasing. More samples in the frequency domain (without changing the sampling rate, the frequency range is still DC to the desired sampling rate for the filter), means more samples in the time domain for the resulting IFFT.  As long as the impulse response is decaying, we can choose a length based on our acceptable distortion limits such that the resulting time domain aliasing is negligible.  With this we get the "true time domain impulse response" needed to proceed with the window design approach for FIR filters. 

Pi-Filter Example

For this example I used a graphic of Pi downloaded from Pinterest.com as the target frequency response (a Pi- shape for fun, after doing this I found out later that Rene Christensen of Acculution used a very similar graphic years earlier to demonstrate epicycles, also very interesting and related to the Fourier Transform.)

I interpolated the graphic to 1501 samples (odd samples are just convenient for the subsequent processing I used, but not necessary), resulting in the starting graphic below as the filter frequency response as viewed on a complex plane (magnitude and phase as frequency is swept from DC to the sampling rate). 

This same response as separate plots of magnitude (in dB) and phase (in radians) is shown below. This represents an arbitrary frequency response for which we would like to design an FIR filter.

As the next step, and done iteratively with the interpolation of the graphic above, I reviewed the resulting impulse response by computing the IFFT (python code shown below along with all imports used, where z_new corresponds to the 1501 complex samples of the target frequency response)

import numpy as np
import scipy.signal as sig
import scipy.fft as fft
impulse_response = fft.fftshift(fft.ifft(z_new))

In particular I am confirming that the impulse response is long enough to have decayed below a certain level of distortion (that level depends on the filter design requirements, but for this example I chose -120 dB to demonstrate how the impulse response behaves, and suggests any time domain aliasing that has folded back into the response area we will ultimately select will be more than insignificant).

Notice the use of 'fftshift' when computing the IFFT. This converts a non-causal time domain result to be causal starting at time index = 0. 

With an acceptable reduction in time domain aliasing (which goes down as we increase the original interpolation in the frequency domain), we can then do the final step as the Window method to truncate and window the impulse response, centered on the dominant area of the response. 

ntaps = 165
center = int(impulse_response.size //2)
coeff = impulse_response[center - ntaps//2: center + ntaps//2 + 1]
coeff *= sig.windows.kaiser(coeff.size, beta=8)

The view the frequency response as first introduced, we need to remove a linear phase in the frequency domain corresponding to the delay of the filter. If we were to implement this filter, we would have a matching delay in a reference path, which is consistent with how 90 degree splitters (Hilbert transforms) are also implemented in practice, and sometimes referred to as phase tracking networks. The frequency response is determined in Python as follows:

w, h = sig.freqz(coeff, whole=True) 
# remove a linear phase corresponding to filter delay
h *= np.exp(1j *2*np.pi* np.linspace(0, int((coeff.size-1)/2), h.size, endpoint=False))

This results in the matching frequency response as first introduced.

If this example caught your interest, you’ll be glad to know I’ll be diving into techniques like this in much more detail at my workshop at the 2025 DSP Online Conference.

Sign in