-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathFilters.py
154 lines (124 loc) · 5.31 KB
/
Filters.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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
import numpy as np
from scipy.signal import kaiserord, firwin
from scipy.signal import savgol_filter
def BandStopFilter():
# minima frequenza, detta frequenza di Nyquist (o anche cadenza di Nyquist), necessaria per campionare un segnale
# analogico senza perdere informazioni
nyq_rate = 1000 / 2
# regione di transizione del filtro dove il filtro passa da attenuare a non attenuare
width = 0.2 / nyq_rate
# oscillazione del filtro
ripple_db = 12
# ordine del filtro, parametro del filtro ottenuto con pigreco * shape of the window
num_of_taps, beta = kaiserord(ripple_db, width)
if num_of_taps % 2 == 0:
num_of_taps = num_of_taps + 1
# intervallo all'interno del quale attenuare
cutoff_hz = np.array([59.5, 60.5])
# filtro
filter_bs = firwin(num_of_taps, cutoff_hz / nyq_rate, window=('kaiser', beta), pass_zero='bandstop')
# valutare il filtro
# w, h = freqz(filter_bs, worN=40000)
#
# plt.plot((w / np.pi) * nyq_rate, 20 * np.log10(np.abs(h)), linewidth=2)
#
# plt.axhline(-ripple_db, linestyle='--', linewidth=1, color='c')
# delta = 10 ** (-ripple_db / 20)
# plt.axhline(20 * np.log10(1 + delta), linestyle='--', linewidth=1, color='r')
# plt.axhline(20 * np.log10(1 - delta), linestyle='--', linewidth=1, color='r')
#
# plt.xlabel('Frequency (Hz)')
# plt.ylabel('Gain (dB)')
# plt.title('Frequency Response')
# plt.xlim(10, 100)
# plt.ylim(-15, 5)
# plt.xscale('log')
# plt.grid(True)
# plt.tight_layout()
# plt.savefig("plot/BandStopFilter.svg", dpi=1200)
return filter_bs
def HighPassFilter():
# minima frequenza, detta frequenza di Nyquist (o anche cadenza di Nyquist), necessaria per campionare un segnale
# analogico senza perdere informazioni
nyq_rate = 1000 / 2
# regione di transizione del filtro dove il filtro passa da attenuare a non attenuare
width = 0.2 / nyq_rate
# oscillazione del filtro
ripple_db = 12
# ordine del filtro, parametro del filtro ottenuto con pigreco * shape of the window
num_of_taps, beta = kaiserord(ripple_db, width)
if num_of_taps % 2 == 0:
num_of_taps = num_of_taps + 1
# frequenza al di sotto della quale attenuare
cutoff_hz = 0.5
# filtro
filter_hf = firwin(num_of_taps, cutoff_hz / nyq_rate, window=('kaiser', beta), pass_zero='highpass')
# valutare il filtro
# w, h = freqz(filter_hf, worN=40000)
#
# plt.plot((w / np.pi) * nyq_rate, 20 * np.log10(np.abs(h)), linewidth=2)
#
# plt.axvline(cutoff_hz + width * nyq_rate, linestyle='--', linewidth=1, color='g')
# plt.axvline(cutoff_hz - width * nyq_rate, linestyle='--', linewidth=1, color='g')
# plt.axhline(-ripple_db, linestyle='--', linewidth=1, color='c')
# delta = 10 ** (-ripple_db / 20)
# plt.axhline(20 * np.log10(1 + delta), linestyle='--', linewidth=1, color='r')
# plt.axhline(20 * np.log10(1 - delta), linestyle='--', linewidth=1, color='r')
#
# plt.xlabel('Frequency (Hz)')
# plt.ylabel('Gain (dB)')
# plt.title('Frequency Response')
# plt.ylim(-20, 5)
# plt.xlim(0.1, 50)
# plt.xscale('log')
# plt.grid(True)
# plt.tight_layout()
# plt.savefig("plot/HighPassFilter.svg", dpi=1200)
return filter_hf
def LowPassFilter():
# minima frequenza, detta frequenza di Nyquist (o anche cadenza di Nyquist), necessaria per campionare un segnale
# analogico senza perdere informazioni
nyq_rate = 1000 / 2
# regione di transizione del filtro dove il filtro passa da attenuare a non attenuare
width = 1 / nyq_rate
# oscillazione del filtro
ripple_db = 12
# ordine del filtro, parametro del filtro ottenuto con pigreco * shape of the window
num_of_taps, beta = kaiserord(ripple_db, width)
if num_of_taps % 2 == 0:
num_of_taps = num_of_taps + 1
# frequenza al di sopra della quale attenuare
cutoff_hz = 100
# filtro
filter_lp = firwin(num_of_taps, cutoff_hz / nyq_rate, window=('kaiser', beta), pass_zero='lowpass')
# valutare il filtro
# w, h = freqz(filter_lp, worN=40000)
#
# plt.plot((w / np.pi) * nyq_rate, 20 * np.log10(np.abs(h)), linewidth=2)
#
# plt.axvline(cutoff_hz + width * nyq_rate, linestyle='--', linewidth=1, color='g')
# plt.axvline(cutoff_hz - width * nyq_rate, linestyle='--', linewidth=1, color='g')
# plt.axhline(-ripple_db, linestyle='--', linewidth=1, color='c')
# delta = 10 ** (-ripple_db / 20)
# plt.axhline(20 * np.log10(1 + delta), linestyle='--', linewidth=1, color='r')
# plt.axhline(20 * np.log10(1 - delta), linestyle='--', linewidth=1, color='r')
#
# plt.xlabel('Frequency (Hz)')
# plt.ylabel('Gain (dB)')
# plt.title('Frequency Response')
# plt.ylim(-15, 5)
# plt.xlim(10, 200)
# plt.xscale('log')
# plt.grid(True)
# plt.tight_layout()
# plt.savefig("plot/LowPassFilter.svg", dpi=1200)
return filter_lp
def SmoothSignal(signal):
smoothed_signal = savgol_filter(signal, 29, 3) # window size 29, polynomial order 3
# plt.plot(signal[:int(len(signal) / 16)], label="raw signal")
# plt.plot(smoothed_signal[:int(len(smoothed_signal) / 16)], label="smoothed signal")
# plt.grid(True)
# plt.legend()
# plt.tight_layout()
# plt.savefig("plot/SmoothSignal.svg", dpi=1200)
return smoothed_signal