1
$\begingroup$

I have been attempting to design a digital narrow band-pass filter with following parameters. Sampling frequency $f_s = 10^4\,\mathrm{Hz}$, low cutoff frequency $f_{cl} = 290\,\mathrm{Hz}$ and high cutoff frequency $f_{ch} = 310\,\mathrm{Hz}$.

After problems with stability of the iir type I have decided to use the fir type filter. Due to the fact that I am very beginner in the dsp (and I don't have Matlab) I have attempted to exploit the DSP Guide and the Scilab wfir command.

The problem which I have is the filter order. Let's say I will use the design approximation from the DSP Guide (Chapter 16) $M\approx\frac{4}{bw}$. For my filter parameters following holds $bw = \frac{310 - 290}{10^4} = 0.002$ i.e. $M \approx \frac{4}{bw} = 2000$. It seems to me to be a huge filter order. Neverthless I have finished the design

Ts = 100*10^(-6); // sampling period (s) fs = 1/Ts; // sampling frequency (Hz) fbp_low = 290; // low cut-off frequency (Hz) fbp_high = 310; // high cut-off frequency (Hz) // difference equation calculation [filter_coeffs] = wfir('bp', 2000, [fbp_low/fs, fbp_high/fs], 'hm', [0, 0]); // coefficients of the polynomial in the numerator of the transfer funtion B = filter_coeffs; // coefficients of the polynomial in the denominator of the transfer funtion A = 1.0; // Frequency response // polynomial in the numerator of the transfer function in z^-1 num = poly(B, 'invz', 'c'); // polynomial in the denominator of the transfer function in z^-1 den = poly(A, 'invz', 'c'); // relative frequency - w*k*Ts = 2*pi*f/fs*k, maximum of the frequency is f = fs/2 // fr_max = f_max/fs = 0.5 fr = (0:0.0001:0.5); // complex frequency response (transfer function in z^(-1)) // z = exp(s*T) = exp([sigma + j*omega]*T) hf = freq(num, den, exp(-%i*2*%pi*fr)); // magnitude magnitude = abs(hf); // phase hf_imag = imag(hf); hf_real = real(hf); phase = atand(hf_imag, hf_real); scf(); plot(fr, magnitude, 'Linewidth', 2); title('Magnitude'); xlabel('$f_r = \frac{f}{f_s}$'); ylabel('$Mag(H(z))\,(\mathrm{-})$'); xgrid; scf(); plot(fr, phase, 'Linewidth', 2); title('Phase'); xlabel('$f_r = \frac{f}{f_s}$'); ylabel('$Arg(H(z))\,(\mathrm{^\circ})$'); xgrid;

and evaluated the frequency response of this filter

magnitude

phase

As far as the magnitude part I would say that the design specifications have been met. But the filter order seems to me to be so huge that it rules out this filter for practical usage (real time control software). May I ask you for an advice how to set the filter order to get a practically realizable narrow band-pass filter which mets the specification?

Improve this question
$\endgroup$
2
  • $\begingroup$ For real-time control, an IIR filter would be a better option. $\endgroup$ Commented Nov 12, 2024 at 13:54
  • 1
    $\begingroup$ It's gonna be a really long FIR. You might wanna go IIR on this. $\endgroup$ Commented Nov 12, 2024 at 16:13

3 Answers 3

Reset to default
1
$\begingroup$

I see this question, or a variant of it, asked quite a bit, so I figure I'd write out an answer. The most narrowband possible filter you can have is a complex exponential (cosine for real-valued signals) centered at frequency $\omega_{0}$. This is really a pretty straightforward proof. A finite length, periodic signal is equivalent to multiplying the infinite length signal by a rectangular window. This means you are convolving the spectrum of the signal with a sinc function which has the effect of broadening the spectrum. For a fixed size rectangular window, the most narrowband signal after windowing is the most narrowband filter prior to windowing, which is a complex exponential at that frequency since its DTFT is a delta function. Thus, an impulse response that is a complex exponential will have the most narrowband frequency response.

The DFT of a finite-length complex exponential is \begin{align} X[k] &= \sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn} \\ &= \sum_{n=0}^{N-1}e^{j\left(\omega_{0}-\frac{2\pi k}{N}\right)n} \\ &= \frac{1-e^{j\left(\omega_{0}-\frac{2\pi k}{N}\right)N}}{1-e^{j\left(\omega_{0}-\frac{2\pi k}{N}\right)}} \\ &= e^{j\left(\omega_{0}-\frac{2\pi k}{N}\right)\frac{N}{2}}\frac{\sin\left(\frac{N}{2}\left(\omega_{0}-\frac{2\pi}{N}k\right)\right)}{\sin\left(\frac{1}{2}\left(\omega_{0}-\frac{2\pi}{N}k\right)\right)} \end{align} The magnitude of this is obviously \begin{equation} \left\lvert X[k]\right\rvert=\left\lvert\frac{\sin\left(\frac{N}{2}\left(\omega_{0}-\frac{2\pi}{N}k\right)\right)}{\sin\left(\frac{1}{2}\left(\omega_{0}-\frac{2\pi}{N}k\right)\right)}\right\rvert \end{equation} Letting $\theta = \frac{1}{2}\left(\omega_{0}-\frac{2\pi}{N}k\right)$, we get \begin{equation} \left\lvert X[k]\right\rvert = \left\lvert\frac{\sin\left(N\theta\right)}{\sin\left(\theta\right)}\right\rvert \end{equation} The peak of this function occurs at $\theta=0$ \begin{equation} \left\lvert X_{max}\right\rvert = \lim_{\theta\to 0}\left\lvert\frac{\sin\left(N\theta\right)}{\sin\left(\theta\right)}\right\rvert=N \end{equation} The half-power beamwidth occurs when \begin{equation} \left\lvert\frac{\sin\left(N\theta\right)}{\sin\left(\theta\right)}\right\rvert=\frac{N}{\sqrt{2}} \end{equation} ie, \begin{equation} \left\lvert\frac{\sin\left(N\theta\right)}{N\sin\left(\theta\right)}\right\rvert=\frac{1}{\sqrt{2}} \end{equation} We can use a Taylor expansion to approximate this \begin{equation} \sin(x) \approx x - \frac{x^{3}}{6} \end{equation} \begin{align} \left\lvert\frac{\sin\left(N\theta\right)}{N\sin\left(\theta\right)}\right\rvert &\approx \left\lvert\frac{N\theta - \frac{\left(N\theta\right)^{3}}{6}}{N\left(\theta-\frac{\theta^{3}}{6}\right)}\right\rvert \\ &= \left\lvert\frac{1-\frac{N^{2}\theta^{2}}{6}}{1-\frac{\theta^{2}}{6}}\right\rvert \end{align} Using a small angle approximation, we can say \begin{equation} \left\lvert X[k]\right\rvert = \left\lvert 1-\frac{N^{2}\theta^{2}}{6}\right\rvert \end{equation} Assume $\theta < \frac{\sqrt{6}}{N}$, we then have \begin{equation} 1-\frac{N^{2}\theta^{2}}{6} = \frac{1}{\sqrt{2}} \end{equation} which means \begin{equation} \theta \approx \frac{1.325}{N} \end{equation} So, if we want a half-power beamwidth of 0.002, we then get \begin{equation} N = \frac{1.325}{0.002} \approx 663 \end{equation} Here is a plot showing, for a 662 length sinusoid, the half-power beamwidth occurs at approximately the right point.

Plot showing frequency response of sinusoid

Now, this is better than 2000, but doesn't take into account windowing, and is still quite long and maybe not suitable for real time, so IIRs may be your only option.

Improve this answer
$\endgroup$
4
$\begingroup$

You've run into an incarnation of the "Fourier Uncertainty Principle" : A signal that is short in one domain will be very wide on the other. The corollary for filter design is: wherever you have a steep frequency transition (relative to the sample rate) you will end up with a very large group delay in that frequency range.

There is really no way around that. Best you can do is to find the best trade off between frequency domain specifications and group delay requirements for your specific application.

Please note, that wfir() implements a linear phase filters. That has constant group delay over all frequencies but the group delay is on the larger side. You get less group delay with a minimum phase filter, but it will be frequency dependent.

Improve this answer
$\endgroup$
4
$\begingroup$

The dependency of filter order on bandwidth is just a fact of life, and a consequence of living in a universe where the Fourier transform represents reality well. You're probably going to have to live with it.

  1. Move to a different universe*.
  2. Make that 2000-point FIR filter.
  3. Fix your IIR filter**.
  4. This will result in the most code, but could result in the least number of execution cycles: running the signal through a number of CIC filter-decimation stages until you're sampling at 1250Hz or even 625Hz, then bandpass, then either use the signal as-is at the lower sampling rate, or resample it back up to 10kHz

* I don't recommend this -- you may miss, and end up in a more mathematically difficult one than this.

** Ask how to fix your IIR filter here -- you're probably underestimating the necessary precision either in the coefficients or the data paths, but it's a separate question.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ There's an extension of #3. You can frequency shift and decimate in steps using fairly simple low-pass filters. Then you can upsample and interpolate to get back to the original sampling rate. It is more computationally efficient than the 2000 point filter, but not as efficient as the CIC. It is a little more general than the CIC approach which needs to be done if fixed point. I'd say it's really just a more generalized approach. It is discussed in R. Crochiere, and L. Rabiner's book on Multirate DSP. $\endgroup$ Commented Nov 15, 2024 at 14:23
  • $\begingroup$ @David: true. Given that your comment exists, I think I'm going to let my answer stand -- but you're still correct. $\endgroup$ Commented Nov 16, 2024 at 4:50
  • $\begingroup$ No need to modify your answer and I didn't want to write a full answer where yours already seems to address the issues. $\endgroup$ Commented Nov 17, 2024 at 14:04

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.