# Building a FIR filter from an arbitrary frequency-magnitude response curve (eg. Least Squares fitting) [duplicate]

I'm trying to create a FIR filter from a magnitude equation, where my starting equation provides the magnitude (amplitude) between 0 and 1 for any given frequency in Hz.

I posted the magnitude equations I would ideally like to convert into an FIR filter here:

https://math.stackexchange.com/questions/3666215/how-to-do-an-inverse-fourier-transform-of-this-equation

It seems there are more than one ways to get the result, but as I have never made a FIR filter before I'm not sure of the best way.

Option #1: Integral-Based Solution

One suggestion I was told was that I need to take the magnitude equation, do an Inverse Fourier Transform on it (which depends on the integrals), and then from that equation I can somehow construct a FIR filter.

But the integrals involve crazy operations like "Ei" and "erf" that probably aren't practical. So to do that, I would need to make a polynomial line of best fit for the magnitude equation and then use that for the reverse Fourier and FIR build.

Option #2: Least Squares Fitting

Alternatively, I read here a different suggestion for what I think is essentially the same problem:

Design a linear-phase FIR filter approximating the magnitude of a given IIR filter

That shares instructions on how to build a FIR filter to copy a IIR filter. If I understand, they used the "freqz" operation in Matlab to create a vector of frequency-vs-magnitude data points from the IIR filter. Those points were then fed into a least squares FIR filter design called "lslevin.m" which would then compute the FIR coefficients (or provide some output needed to get the coefficients).

If I'm reading that correctly, that approach would entirely solve my problem. I wouldn't even need the "freqz" operation, since I already have the amplitude/magnitudes at every frequency. I can just input those directly to the least squares FIR filter algorithm. It would also be linear phase which would be ideal.

Am I understanding this correctly, and if so, will that give me what I need?

Option #3: Directly Mapping Frequency-Amplitudes to FIR Bins

As a third option, I was recently told another approach:

If you can get a value for the desired magnitude at any given frequency, the easiest approach is to sample the spectrum: simply compute an array of desired magnitudes with frequencies matching the FFT bin-frequencies for an FFT of your array size. No integration required. Essentially you just build the array that you would like to get out of FFT if you had the perfect filter.

Then you can take IFFT to get a (periodic) signal that matches the desired spectrum at the compute bin-frequencies exactly. At this point, you want to apply some window function (with basically the same considerations as for any other windowed filter design), so you get a single IR (rather than a cycle of a periodic signal). This will smooth out the response slightly, so it won't be exact anymore, but assuming your original magnitude is reasonably smooth with respect to the length of the filter (ie. the frequency samples we took don't vary widely from one bin to the next) the windowed approximation is usually decent and you can make it better by computing a longer filter (by sampling more points and computing a longer IFFT).

This type of design isn't "optimal" in any sense, but it's nice and simple and it works as long as you can tell the desired magnitude at any given frequency.

Is that going to be the same or different from the least squares method and in what ways? Which approach will be cheaper or more accurate?

• Yes I believe you are reading it correctly. May 9 '20 at 11:28
• @lennon310 that's the exact page I linked in my question, so that's what I was asking about. Thanks though. Sounds like it is what I need.
– mike
May 9 '20 at 12:50
• So why haven't you yet tried it? May 9 '20 at 18:39
• @Matt L. I'm a beginner with FIR filters and I have zero education in this subject so these things take me a lot of time to figure out. Also I was just told that alternatively I should just take my magnitude equation to assign the values directly to the bin frequencies of my FIR based on the number of taps I have, and I can somehow use the FFT that to get my coefficients that way more directly. Is that a better method than the least squares one or the same?
– mike
May 10 '20 at 6:00
• I added the third option to the question. Any guidance is appreciated.
– mike
May 10 '20 at 6:06

Am I understanding this correctly, and if so, will that give me what I need?

Yes you are understanding everything correctly. If your "need" is a filter that will match the magnitude response then all choices should work. If you are interested in a very quick and move-on solution, then proceed with Frequency Sampling with a very large number of taps and move on. If you actually have a concern with efficiency (meaning you will be implementing the filter in hardware), then I recommend Least-Squares (although you could run into convergence issues). If your interest is education, then try all and compare the results yourself as nothing will beat the first-hand experience. All approaches will result in a solution.

The descriptions of the different approaches for designing FIR filters as stated by the OP are reasonably accurate, and given no constraint on total number of taps (and total delay) all could achieve the desired objective, so often the decision is just a trade of how much time you want to spend on designing the filter to achieve the best efficiency (least number of taps).

A fourth approach not mentioned is the minimax solution (equiripple) using the Remez exchange algorithm is also a common design approach.

The graphic below lays out my flow-chart for the typical steps in the FIR design process, with further details specific to estimating the order of the filter and the 4 major design approaches. Estimating Order: Estimating the number of taps required in the FIR filter is determined before any further design efforts are made in any of the approaches. The number of taps required is driven from the derivative of the filter response versus frequency, both in magnitude and phase. The derivative of the magnitude response for typical low pass or high pass filter structures translates to very tight transition bands require a large number of taps, and fred harris' rule which is my favorite due to it's simplicity:

$$N \approx \frac{A}{22}\frac{2\pi}{\Delta\omega}$$

Where:
$$A$$: filter attenuation in dB
$$\Delta\omega$$: Transition band of filter in radians/sample ($$2\pi =$$ sampling rate).

The OP's case isn't a simple low pass or high pass filter, but this illuminates that the tighter the features are in frequency, the more taps will be needed.) Similarly the negative derivative of phase with respect to frequency is the group delay of the filter by definition. The OP does not care about phase but for general purposes we see how the slope of the phase vs frequency also drives the number of taps since the filter delay must be within the span of the filter.

The approaches are further detailed below by their common name, but as far as tap estimation, the approach above is reasonably close for the windowing, least squares, and equi-ripple design approaches. For frequency sampling a significantly higher number of taps may be needed (could be 2x - 3x).

It would be common to start with a much higher number $$N$$ and then reduced the number of taps once the relative magnitude of the tails at the edge of the filter is reviewed. This can be reduced by pursuing the same algorithm with a smaller $$N$$, or truncating and windowing the solution achieved with a higher $$N$$.

Additional details on estimating the number of taps needed:

How many taps does an FIR filter need?

Design Coefficients:

The four common approaches for designing the FIR filter coefficients from target specifications are as follows:

• Windowing (OP Option #1)
• Frequency Sampling (OP Option #3)
• Equiripple (Parks-McClellan, Remez Exchange Algorithm)
• Least Squares (OP Option #2)

Windowing: (OP Option #1)

This method is to window the sampled continuous-time impulse response of the filter. The impulse response is determined from the inverse Fourier Transform of the frequency response, which is sampled and delayed to be causal and windowed to select the dominant response features based on the number of taps used. Just selecting samples of the impulse response with no further modification is multiplying the impulse response by a rectangular window in time. Multiplication in time is convolution in frequency, so the desired frequency response will be convolved with a Sinc function in frequency. The Sinc function has the tightest main lobe for a given time span ($$N$$ samples), but has relatively high sidelobes. This results in the closest match of transitions but would have more ripple in the passband and stopband compared to specialized windows.

This is demonstrated below with an example of designing a perfect (target design, not achieved) low-pass filter. The impulse response of the low pass filter is a Sinc function, which is delayed in time in order to be causal (resulting in a linear negative phase slope versus frequency) and then truncated with a rectangular window to given the number of taps desired. Once truncated, the other Sinc function in frequency due to the window (not the Sinc in time due to the desired filter response; in this case slightly confusing since two Sinc functions are involved) will convolve with the desired response, widening the bandwidth slightly and causing ripple in the passband and stopband. This is often improved by simply multiplying the desired impulse response with a different window shape that by design has very low side-lobes. The trade of doing this is a wider main-lobe which will serve to shift the exact locations of frequency transitions but will provide a much better match in areas where the frequency response is not changing rapidly. Increasing the number of taps resolves both cases.

Summary: Windowing is quick and easy if the inverse Fourier Transform of the target response is easily solved, and provides a reasonable result. It is sub-optimal in that other design approaches can achieve closer accuracy to target specifications with less number of taps. Frequency Sampling: (OP Option #3)

This is the method I mostly discourage anyone from using (although I ultimately used that to demonstrate a quick approach for the OP's more challenging filter at his other post here: How to invert FFT these magnitude equations to get FIR filters?. showing that it does have utility if a quick and useful result is needed).

The benefit of this approach is that it is the simplest: You use the inverse DFT to create the coefficients of the filter. This will result in an exact match at the DFT sample locations, but the resulting response will have significantly higher ripple (deviation) for all frequencies in between.

A hybrid approach that would be very viable in this specific application of the OP's is to zero pad the high frequency portion of the target frequency response (which is the middle of the DFT array) allowing for a much longer duration for the impulse response (eliminating time domain aliasing). The resulting impulse response can then be windowed as in the window approach given above. Optionally the achieved impulse response can be resampled and then truncated/windowed. Basically we would be using a longer Inverse DFT to approximate the Inverse Fourier Transform when the result would be too difficult to compute directly (as in the OP's case).

When deriving the coefficients from the ifft, the frequency response needs to include the "negative frequency" components: for a filter with real coefficients this would be a mirror of the DFT samples from $$0$$ to $$N/2-1$$ at $$N/2$$ to $$N-1$$. Further, when the magnitude response is of only concern, the following will properly center the impulse response in the span of the filter as depicted in the diagrams above:

coeff = fftshift(ifft(ifftshift(mag_response)));


Here mag_response was generate with the frequency axis going from $$-f_s/2$$ to $$+f_s/2$$ so the inner ifftshift translates it to the format expected in the DFT going from $$0$$ to $$f_s$$, where $$f_s$$ is the sampling rate. (corresponding to $$0$$ to $$N-1$$). The final fftshift converts the noncausal result in the time domain (the coefficients are centered on $$t=0$$) to a causal filter by centering the coefficients in the center of the filter span.

Filtering by zeroing out FFT bins is the frequency sampling approach, which is further described here including why it is generally a bad idea (meaning will require a lot more taps than the optimal approaches equiripple and least-squares):

Why is it a bad idea to filter by zeroing out FFT bins?

Equiripple:

This is the go-to design algorithm when the desired accuracy result is constrained by peak limits. The equiripple design approach developed by Thomas Parks and James McClellan in 1972 https://en.wikipedia.org/wiki/Parks%E2%80%93McClellan_filter_design_algorithm uses the Remez-exchange algorithm resulting in an optimum filter design given a minimax constraint (minimize the maximum error). Optimum meaning the solution will be the mimimum number of taps given the peak error design constraint. It can be used to match arbitrary shapes and multi-band filters and is supported by functions in MATLAB (firpm), Octave (remez) and Python (scipy.signal.remez). For a more complicated filter it may have challenges converging which would be my only reason to pursue one of the earlier approaches above.

Least-Squares: (OP Option #2)

This is the go-to design algorithm when the desired accuracy result is constrained by root-mean-square (rms) limits. Given the performance of most applications where I need to use filters (wireless communications) is optimized by rms performance and not peak error, this is my filter design of choice. Like the equiripple design, this is an optimal filter design in that the solution will be the mimimum number of taps given the rms design constraint. Unless I was effectively given a "though shall not cross this line" constraint in the filter specifications, I would use least-squares. It can be used to match arbitrary shapes and multi-band filters and is supported by functions in MATLAB/Octave (firls), and Python (scipy.signal.firls). For a more complicated filter it may have challenges converging which would be my only reason to pursue one of the non-optimal design approaches above.

Further posts of interest:

FIR Filter Design: Window vs Parks McClellan and Least Squares