77import  matplotlib .pylab  as  plt 
88import  numpy  as  np 
99from  scipy  import  signal 
10+ import  iir 
1011
11- SAMPLE_RATE  =  2000   # Hertz 
12- DURATION  =  5   # Seconds 
13- NTAPS  =  1000 
14- SAMPLES  =  DURATION  *  SAMPLE_RATE 
15- fs  =  SAMPLE_RATE 
12+ fs  =  2000   # Hertz 
1613
1714def  generate_sine_wave (freq , sample_rate , duration , magnitude ):
18-     x  =  np .linspace (0 , duration , sample_rate  *  duration , endpoint = False )
15+     x  =  np .linspace (0 , duration , fs  *  duration , endpoint = False )
1916    frequencies  =  x  *  freq 
2017    y  =  magnitude  *  np .sin ((2  *  np .pi ) *  frequencies )
2118    return  x , y 
22- 
23- class  IIR2Filter (object ):
24-     """ 
25-     given a set of coefficients for the IIR filter this class creates an object that keeps the variables needed for the 
26-     IIR filtering as well as creating the function "filter(x)" with filters the input data. 
27-     Attributes: 
28-         @:param coefficients: input coefficients of the IIR filter as an array of 6 elements where the first three 
29-         coefficients are for the FIR filter part of the IIR filter and the last three coefficients are for the IIR part 
30-         of the IIR filter. 
31-     """ 
32- 
33-     def  __init__ (self , coefficients ):
34-         self .myCoefficients  =  coefficients 
35-         self .IIRcoeff  =  self .myCoefficients [3 :6 ]
36-         self .FIRcoeff  =  self .myCoefficients [0 :3 ]
37-         self .acc_input  =  0 
38-         self .acc_output  =  0 
39-         self .buffer1  =  0 
40-         self .buffer2  =  0 
41-         self .input  =  0 
42-         self .output  =  0 
43- 
44-     def  filter (self , input ):
45-         """ 
46-         :param input: input value to be processed. 
47-         :return: processed value. 
48-         """ 
49- 
50-         self .input  =  input 
51-         self .output  =  0 
52- 
53-         """ 
54-         IIR Part of the filter: 
55-             The accumulated input are the values of the IIR coefficients multiplied by the variables of the filter:  
56-             the input and the delay lines. 
57-         """ 
58-         self .acc_input  =  (self .input  +  self .buffer1 
59-                           *  - self .IIRcoeff [1 ] +  self .buffer2  *  - self .IIRcoeff [2 ])
60- 
61-         """ 
62-         FIR Part of the filter:      
63-             The accumulated output are the values of the FIR coefficients multiplied by the variables of the filter:  
64-             the input and the delay lines.  
65-          
66-         """ 
67- 
68-         self .acc_output  =  (self .acc_input  *  self .FIRcoeff [0 ]
69-                            +  self .buffer1  *  self .FIRcoeff [1 ] +  self .buffer2 
70-                            *  self .FIRcoeff [2 ])
71- 
72-         # Shifting the values on the delay line: acc_input->buffer1->buffer2 
73-         self .buffer2  =  self .buffer1 
74-         self .buffer1  =  self .acc_input 
75-         self .input  =  self .acc_output 
76-         self .output  =  self .acc_output 
77-         return  self .output 
7819
20+ # Generate 2 sinusoidal waves with a frequencies of 1 and 50 Hz for test 
7921f1  =  1 
8022f2  =  50 
8123magnitude_1  =  1 
8224magnitude_2  =  0.5 
83- x , sine1  =  generate_sine_wave (f1 , SAMPLE_RATE , DURATION , magnitude_1 )
84- _ , noise  =  generate_sine_wave (f2 , SAMPLE_RATE , DURATION , magnitude_2 )
25+ DURATION  =  5   # Seconds 
26+ x , sine1  =  generate_sine_wave (f1 , fs , DURATION , magnitude_1 )
27+ _ , noise  =  generate_sine_wave (f2 , fs , DURATION , magnitude_2 )
8528mysignal  =  sine1  +  noise 
29+ 
8630plt .figure (1 )
8731plt .plot (mysignal )
8832
8933# ----------- IIR filter --------------  
90- cutoff  =  [0.5 , 45 ]
9134order  =  20 
9235fs  =   2000 
93- for   i   in   range ( len ( cutoff )): 
94-      cutoff [ i ]  =   cutoff [ i ] / (fs  *  2 )
36+ f1   =   0.5   / ( fs   *   2 ) 
37+ f2   =   45   / (fs  *  2 )
9538
9639cutoff  = [0.000125 , 0.01125 ] 
40+ 
9741coeff  =  signal .butter (order , cutoff , 'bandpass' , output = 'sos' )
9842
9943# If the order of the filter is 1, is one IIR 2nd order filter otherwise, it is a chain of IIR filters. 
44+ SAMPLES  =  DURATION  *  fs 
10045
101- myFilter  =  IIR2Filter (coeff [0 ])
102- 
46+ myFilter  =  iir .IIR_filter (coeff )
10347y  =  np .zeros (SAMPLES )
10448for  i  in  range (SAMPLES ):
10549    y [i ] =  myFilter .filter (mysignal [i ])
10650
107- # np.savetxt('coefficients.dat', b, fmt='%f', delimiter='') 
108- 
51+ # Now plot original and filtered signal 
10952plt .figure (2 )
11053plt .subplot (211 )
11154plt .plot (mysignal )
@@ -114,20 +57,19 @@ def filter(self, input):
11457plt .plot (y )
11558# plt.ylim((-150, 300)) 
11659
60+ # Calculate the fourier trasnform  
11761unfilteredfft  =  np .fft .fft (mysignal )
11862IIRfilteredfft  =  np .fft .fft (y )
11963
120- T  =  1 / SAMPLE_RATE 
121- 
122- # 1/T = frequency 
123- f  =  np .linspace (0 , SAMPLE_RATE , SAMPLES )
64+ T  =  1 / fs 
65+ f  =  np .linspace (0 , fs , SAMPLES )
12466
12567plt .figure (3 )
12668plt .subplot (211 )
127- plt .plot (f [:SAMPLES  //  2 ], np .abs (unfilteredfft )[:SAMPLES  //  2 ] *  1  /  SAMPLES )  # 1 / N is a normalization factor 
128- # plt.plot(20.0*np.log10(unfilteredfft)) 
69+ plt .plot (f [:SAMPLES  //  2 ], np .abs (unfilteredfft )[:SAMPLES  //  2 ] *  1  /  SAMPLES )  
12970plt .subplot (212 )
130- plt .plot (f [:SAMPLES  //  2 ], np .abs (IIRfilteredfft )[:SAMPLES  //  2 ] *  1  /  SAMPLES )  # 1 / N is a normalization factor 
71+ plt .plot (f [:SAMPLES  //  2 ], np .abs (IIRfilteredfft )[:SAMPLES  //  2 ] *  1  /  SAMPLES ) # 1 / N is a normalization factor 
72+ 
73+ np .savetxt ('coefficients.dat' , b , fmt = '%f' , delimiter = '' )
13174
132- print (coeff )
13375plt .show ()
0 commit comments