How to invert FFT these magnitude equations to get FIR filters?

Question

I need filter designs that will reproduce or approximate each of the following:

1) ${\rm magnitude} = e^{-cf^{2}}$

2) ${\rm magnitude} = e^{-\frac{a+b\sqrt{gf}}{df}}$

Where: $a$, $b$, $c$, $d$ and $g$ are constants $>0$, and $f$ is frequency in Hz.

Or if it's easier or an equivalent problem, the multiple of the two filters as one combined filter would even be better:

3) ${\rm magnitude} = e^{-cf^{2} -\frac{a+b\sqrt{gf}}{df}}$

I have never done an inverse FFT or designed a FIR filter though so I'm looking for any basic guidance on how this might be done.

This is for physical modeling so the closer I can get the better. Accuracy is only important with audio frequencies (20-20,000 Hz).

I really need something that will work one way or another. Any guidance is appreciated.

Answer 1

As @RBJ commented the first one is a Gaussian filter. The Fourier Transform of a Gaussian is a Gaussian so this one can be easily created by sampling the impulse response given by the inverse Fourier Transform.

For case (1), an FIR filter with a magnitude response given by:

$$|H(f)| = e^{-cf^2}\tag{1}\label{1}$$

Given the Fourier Transform of a Gaussian:

$$ e^{-\pi t^2} \leftrightarrow e^{-\pi f^2}\tag{2}\label{2}$$

And using the Fourier Transform scaling property:

$$\frac{1}{|a|}x(t/a) \leftrightarrow X(af) \tag{3}\label{3}$$

We can rewrite $\ref{1}$ as:

$$X(af) = |H(af)| = e^{-\pi (af)^2}\tag{4}\label{4}$$

with $a = \sqrt{c/\pi}$

If we restrict $c$ to be a positive real number, then from $\ref{3}$ we get:

$$\mathscr{F}^{-1}\{ e^{\pi(af)^2}\}=\frac{1}{|a|}e^{-\pi(t/a)^2} = \sqrt{\frac{\pi}{c}}\left( e^{\frac{(\pi t)^2}{c}}\right)\tag{5}\label{5}$$

This is demonstrated with an increasing number of coefficients showing the excellent match that is achieved.

The MATLAB/Octave code used to create the filter coefficients was:

fs = 1;   # sampling rate
c = 100;
N = 27;   # number of taps, must be odd
t = [-(N-1)/2: (N-1)/2]/fs;   # time axis
coeff = sqrt(pi/c)*exp(-pi^2*t.^2/c); # sampled imp resp

Case (2) appears to be a high-pass filter which could be solved using the approach above by sampling the inverse Fourier Transform for $e^{-a/f}$ and $e^{-a/\sqrt{f}}$.

I couldn't easily proceed from this point and posted on the mathematics site for a possible solution:

https://math.stackexchange.com/questions/3666703/inverse-fourier-transform-for-e-a-f-and-e-a-sqrtf

A very simple but sub-optimum approach is to derive the filter coefficients from the inverse DFT. This is the frequency sampling approach to FIR design and provides an exact match at the frequency sampling locations, but results in significantly more ripple than optimized algorithms such as least-square and Parks-McLellan. The ripple can be significantly reduced by windowing the resulting coefficients, but for this application any windowing will significantly reduce the null at DC. However, in the spirit of "needing something that will work one way or another" the Frequency Sampling approach can be considered if a larger number of taps is feasible. Regardless of approach, the number of taps for an FIR filter is proportional to the inverse of the width of the transition band of the filter, so if the parameters are such that the filter is a tight DC-notch filter, a large number of taps with any traditional FIR approach will be necessary (an IIR filter may make a lot more sense, or a modification of Richard Lyons clever implementation of a linear=phase DC notch filter with CIC structure: https://www.dsprelated.com/showarticle/58.php ).

The example plots below show the result for a 275 and 501 tap FIR filter using the Frequency Sampling approach with the following example parameters as given in the code below to implement an FIR for Case (2):

fs = 1;
a = 500; b = 50; d = 50000; g = 100;
N = 275;
f = linspace(-fs/2, fs/2, N);
magC2 = e.^(-(a+b.*(sqrt(g*abs(f))))./(d*abs(f)));
coeffC2 = fftshift(ifft(ifftshift(magC2)));

The above plot appears to show an excellent match in both cases, however zooming in on the stop band shows that the frequency sampling approach only maintains a match (with the tighter parameters used) down to around -25 dB with 275 taps, while it maintains a match down to around -65 dB with 501 taps. Pursuing either an approach of sampling the computed impulse response from the inverse Fourier Transform, or a least squares algorithm to match the magnitude response will allow for a solution with significantly fewer taps. Also if the notch needed is actually not as tight as I have demonstrated here, the number of taps required will be significantly reduced as well.

Also observe how we can estimate the number of taps needed for any approach from the ideal in the plot above: if we want/need to maintain high accuracy out to -80 dB rejection, we observe for this case with the particular parameters chosen that the steepest transition the slope is approximately 65 dB over a frequency of 0.005 samples/cycle. From fred harris' rule of thumb this would be:

$$N \approx \frac{A}{22}\frac{1}{\Delta f} = (65/22)(1/0.005) = 590 $$

We can also see how increasingly challenging it will be to maintain a tight match at deeper rejections given the slope is continuously increasing.