0
$\begingroup$

I want to design a sinc filter in the time domain (a bandpass filter in the frequency domain).

I wrote test code in python to investigate, and would appreciate help understanding it and answer my question. Here is my script in python:

import numpy as np
from matplotlib import pyplot as plt

fs=1200
ww = np.linspace(-4, 4, fs)
a = 2*25 *np.sinc(2*25*ww)
b= 2*3*np.sinc(2*3*ww)
#a = a/sum(a)
#b = b/sum(b)
ww2 = a-b
#ww2 = 2*25 *np.sinc(2*25*ww) - 2*3*np.sinc(2*3*ww)
print(np.sum(ww2))
#ww2= ww2/np.sum(ww2)
yf = np.fft.fft(ww2)
N = len(yf)
n = np.arange(N)
T = N * fs
freq = N* n / T
n_oneside = N // 2
f_oneside = freq[:n_oneside]

plt.plot(f_oneside, np.abs(yf[:n_oneside]), 'b')
plt.show()

This is written based on this wikipedia link about Sinc filters.

But the point is when I plot the output, I get a gained filter (around 150), not a filter with gain (1). I appreciate if someone can explain this to me.

enter image description here

$\endgroup$

2 Answers 2

1
$\begingroup$

One problem I see is the very low bandwidth edges compared to the sampling frequency. As long as you keep to this being an example here, you shouldn't have too many problems but, you should know that it will cause numeric errors due to cosines being too close to 1.

At any rate, just because you're designing in the time domain it doesn't mean you need to think in terms of the impulse length in seconds. Also, you are calculating the sinc() with frequencies that exceed 0.5, you should divide those to fs. The result you get is not only an aliased version but, also a large gain response (you're multiplying each sinc() with 2*3 and 2*25, so a gain of 600).

Try thinking in terms of samples, the length will come out implicitly. Judging by the way you calculate ww, you want a 1200 taps filter:

fs = 1200
# normalize the frequencies to fs
w1, w2 = 2*3/fs, 2*25/fs
# if you need 1201 taps use (-600,600)
t = np.linspace(-595.5, 595.5, 1200)
a, b = w1*np.sinc(w1*t), w2*np.sinc(w2*t)
c = a - b
# large FFT to account for the large N
Nfft = 32768
C = np.fft.fft(c, Nfft)
w = np.linspace(0, fs, Nfft)
pp.plot(w, abs(C))
pp.show()

bandpass FFT

$\endgroup$
2
  • $\begingroup$ Thank you. So, I appreciate if you tell me I understand the following statement correctly or not. For example, when we want to design a filter in time domain. we need to consider 'fs' taps in time domain so what will happen in frequency domain. because the descrete furier transform is periodic. Thus, we only consider it in one period (length = $\frac{f_s}{2}$). Is this two compatible ? I mean having 'fs' taps in time domain and then our maximum frequency based on Nyquist rate in frequency domain will be $\frac{f_s}{2}$? $\endgroup$
    – Captain
    Commented Nov 4, 2022 at 15:13
  • $\begingroup$ @Captain If you impose the impulse length T then the order will be N=T*fs. If you impose the order then the impulse length will be T=N/fs. If you want to impose both T and N then fs will have to be fs=N/T. In the OP you imposed a (non-causal) imp. length of 8s (your ww, here T; strange notation, BTW) and fs=1200, but the spacing between samples is not 1/fs, but 8/fs. I chose T to be 1 s but, if you want either N or T changed, you have to do it acording to the relations that I said in the beginning. The fs/2 limitation is inherent to sampling, in general. $\endgroup$ Commented Nov 4, 2022 at 17:18
1
$\begingroup$

First of all, you need to:

  1. Clarify your requirements (What is your sampling frequency? What is the center frequency and bandwidth of the BP filter?)
  2. Use a more convenient notation (for instance, you are using fs to denote a number of samples, not a frequency, and ww for a time axis, which is confusing).

I assume in the following that your sampling frequency is $f_S=50\, \mathrm{kHz}$ and you need to design a BP filter centered at $f_C=7\,\mathrm{kHz}$ with (one-sided) bandwidth $B=5.5\,\mathrm{kHz}$.

This results in (one-sided) bandwidth $B_p=12.5\,\mathrm{kHz}$ and $B_n=1.5\,\mathrm{kHz}$ for the positive and the negative sinc components repectively, which should be consistent with your code.

fs = 50e3
fc = 7e3
B  = 5.5e3

Bp = (fc+B)
Bn = (fc-B)

Now, assume that your filter shall have 511 taps (chosen odd for symmetry) and that you will be using a FFT of size 1024 (this implies zero padding from 511 to 1024). With this information, you can define a time axis (in seconds) and a frequency axis (in Hertz).

N = 511
t_axis = (np.arange(N)-(N-1)/2)/fs

Nfft = 1024
f_axis = np.arange(Nfft)/Nfft*fs

Now we can finally answer your question. The formula you found in Wikipedia was for a continuous time filter. When working with digital signal processing, you have to take the discretization into account, which changes the scaling factors.

For a continuous time sinc filter, the value at $t=0$ equals the integral of the frequency response:

$$ h(0)=\int_{-\infty}^{\infty}H(f)\,df=2B $$

since the rectangle has basis $2B$ and height $1$.

In the discrete world, the analogous property is

$$ h[0]=\frac{1}{N_{fft}}\sum_{k=0}^{N_{fft}-1}H[k] $$

Instead of a rectangle, you can now visualize $N_{fft}$ discrete bins, which take either the value 0 or 1. The result of the sum in the above equation is thus just the number of non-zero frequency bins, which depends on the ratio between $2B$ and $f_S$.

$$ h[0] = \frac{1}{N_{fft}} \left\lfloor \frac{2B}{f_S} \cdot N_{fft} \right\rfloor \approx \frac{2B}{f_S} $$

Note that I have used the symbol $N_{fft}$ to keep the reference to the Python script. Recall that FFT is just an efficient algorithm to compute the DFT, but the result is obviously independent of the implementation.

a = 2*Bp/fs * np.sinc(2*Bp*t_axis)
b = 2*Bn/fs * np.sinc(2*Bn*t_axis)
h = a-b

yf = np.fft.fft(h,Nfft)
n_oneside = Nfft // 2
f_oneside = f_axis[:n_oneside]

plt.plot(f_oneside/1e3, np.abs(yf[:n_oneside]), 'b')
plt.xlabel("Frequency [kHz]")
plt.grid()
plt.xlim([0,fs/2e3])
plt.show()

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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