forked from TcheL/Road2Filter
-
Notifications
You must be signed in to change notification settings - Fork 0
/
lpFIR.py
executable file
·48 lines (33 loc) · 1.17 KB
/
lpFIR.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#!/usr/bin/env python3
import numpy as np
def lowPassFIR(fc, b = 0.08):
#fc: Cutoff frequency as a fraction of the sampling rate (in (0, 0.5)).
#b: Transition band, as a fraction of the sampling rate (in (0, 0.5)).
N = int(np.ceil((4 / b)))
if not N % 2: N += 1 # Make sure that N is odd.
n = np.arange(N)
# Compute sinc filter.
h = np.sinc(2 * fc * (n - (N - 1) / 2))
# Compute Blackman window.
w = 0.42 - 0.5 * np.cos(2 * np.pi * n / (N - 1)) + \
0.08 * np.cos(4 * np.pi * n / (N - 1))
# Multiply sinc filter by window.
h = h * w
# Normalize to get unity gain.
h = h / np.sum(h)
return N, h
#-------------------------------------------------------------------------------
import matplotlib.pyplot as plt
fs = 1000 # sampling frequency
# generate properly the time vector
t = np.arange(1000)/fs
sga = np.sin(2*np.pi*2*t) # signal with f = 2
sgb = np.sin(2*np.pi*10*t) # signal with f = 10
sgc = sga + sgb
N, h = lowPassFIR(fc = 6/fs, b = 2/fs)
sgf = np.convolve(sgc, h)
plt.plot(t, sgc, label = 'original')
plt.plot(t, sgf[int(N/2):1000 + int(N/2)], label = 'low-pass')
plt.plot(t, sga, label = 'f = 2')
plt.legend()
plt.show()