0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ can you make your graphic into a .png or .jpg and then upload it into your question? you don't want to inconvenience other people that might help you. $\endgroup$ – robert bristow-johnson May 8 at 19:05
  • $\begingroup$ these are called "gaussian filters" and there is a way to approximate them, sorta, with an IIR filter but i don't remember the details. if you're doing this with an FIR filter, it's straight forward, if you make the FIR long enough. just iFFT your frequency response and get your impulse response. the Fourier Transform of a gaussian is also a gaussian. $\endgroup$ – robert bristow-johnson May 8 at 19:07
  • 1
  • $\begingroup$ Thanks RBJ. I think I need to ask something different then, so I will reword it. It sounds like my question is then just how do I make a FIR filter from a magnitude plot, which I need to learn how to do. Seems the basics are covered here: dspcan.homestead.com/files/Fir/fourier_series_intro.htm $\endgroup$ – mike May 8 at 20:35
  • $\begingroup$ What would be your sampling rate and typical values for a, b, c, d and g? If the range for a,b,c,d and g are constrained there may be a simple solution $\endgroup$ – Dan Boschen May 9 at 16:44
1
$\begingroup$

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.

Frequency Response

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)));

Notch

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.

Zoom in

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.