Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

updating to use the correct sign for the iir 'a' coefficients

  • Loading branch information...
commit 3aa2720244c20cd35ae9f0b6beea265738faf2a0 1 parent 3c6404e
Adam Cox authored September 12, 2011
74  runTestBpAndCorr.py
@@ -6,12 +6,15 @@
6 6
 
7 7
 import matplotlib.pyplot as plt
8 8
 import numpy as np
  9
+
9 10
 noise = np.array(r['noise'])
10 11
 noisepower = np.array(r['noisepower'])
11 12
 template = np.array(r['template'])
12 13
 templatepower = np.array(r['templatepower'])
13 14
 signal = np.array(r['signal'])
14 15
 signalpower = np.array(r['signalpower'])
  16
+windowsignal = np.array(r['window_signal'])
  17
+windowtemplate = np.array(r['window_template'])
15 18
 bp_noise = np.array(r['bp_noise'])
16 19
 bp_noisepower = np.array(r['bp_noisepower'])
17 20
 bp_template = np.array(r['bp_template'])
@@ -19,33 +22,58 @@
19 22
 bp_signal = np.array(r['bp_signal'])
20 23
 bp_signalpower = np.array(r['bp_signalpower'])
21 24
 
  25
+#lets make the horizontal axis of the power spectrums show the correct frequency
  26
+hzaxis = np.array(range(0,len(signalpower)))
  27
+hzperpoint = r['nyquist']/len(signalpower)
  28
+hzaxis *= hzperpoint
  29
+
22 30
 h, w = sig.freqz(r['b'],r['a'])
23 31
 
24  
-corr = np.array(r['corr'])
  32
+corrL = np.array(r['corrLeft'])
  33
+corrM = np.array(r['corrMiddle'])
  34
+corrR = np.array(r['corrRight'])
  35
+
  36
+bp_tempLeft = np.array(r['bp_tempLeft'])
  37
+bp_tempMiddle = np.array(r['bp_tempMiddle'])
  38
+bp_tempRight = np.array(r['bp_tempRight'])
  39
+
25 40
 
26  
-f, axarr = plt.subplots(5, 3)
  41
+f, axarr = plt.subplots(7, 3)
27 42
 axarr[0,0].plot(noise)
28 43
 axarr[0,1].plot(template)
29 44
 axarr[0,2].plot(signal)
30  
-axarr[1,0].semilogy(noisepower)
31  
-axarr[1,1].semilogy(templatepower)
32  
-axarr[1,2].semilogy(signalpower)
33  
-axarr[2,0].plot(bp_noise)
34  
-axarr[2,1].plot(bp_template)
35  
-axarr[2,2].plot(bp_signal)
36  
-axarr[3,0].semilogy(bp_noisepower)
37  
-axarr[3,1].semilogy(bp_templatepower)
38  
-axarr[3,2].semilogy(bp_signalpower)
39  
-axarr[4,0].title('Digital filter frequency response')
40  
-axarr[4,0].title('Digital filter frequency response')
41  
-axarr[4,0].semilogy(h, np.abs(w), 'b')
42  
-axarr[4,0].ylabel('Amplitude (dB)', color='b')
43  
-#axarr[4,0].xlabel('Frequency (rad/sample)')
44  
-axarr[4,0].grid()
45  
-axarr[4,0].legend()
46  
-
47  
-ax2 = ax1.twinx()
  45
+axarr[1,0].loglog(hzaxis, noisepower)
  46
+axarr[1,1].loglog(hzaxis, templatepower)
  47
+axarr[1,2].loglog(hzaxis, signalpower)
  48
+#axarr[2,0].title('Digital filter frequency response')
  49
+#axarr[2,0].title('Digital filter frequency response')
  50
+
  51
+#axarr[2,0].ylabel('Amplitude (dB)', color='b')
  52
+#axarr[2,0].xlabel('Frequency (rad/sample)')
  53
+#axarr[2,0].grid()
  54
+#axarr[2,0].legend()
  55
+#ax1 = axarr[2,0].add_subplot(111)
  56
+#ax2 = ax1.twinx()
  57
+# hzperpoint = tr['nyquist']/len(signalpower)
  58
+hzperpoint = r['nyquist']/h[len(h)-1]
  59
+h *= hzperpoint
  60
+axarr[2,0].loglog(h, np.abs(w), 'b')
48 61
 angles = np.unwrap(np.angle(w))
49  
-axarr[4,0].ylabel('Angle (radians)', color='g')
50  
-axarr[4,0].plot(h, angles, 'g')
51  
-axarr[4,2].plot(corr)
  62
+#axarr[2,1].ylabel('Angle (radians)', color='g')
  63
+#axarr[2,1].plot(h, angles, 'g')
  64
+axarr[2,1].plot(windowtemplate)
  65
+axarr[2,2].plot(windowsignal)
  66
+axarr[3,0].plot(bp_noise)
  67
+axarr[3,1].plot(bp_template)
  68
+axarr[3,2].plot(bp_signal)
  69
+axarr[4,0].loglog(hzaxis, bp_noisepower)
  70
+axarr[4,1].loglog(hzaxis, bp_templatepower)
  71
+axarr[4,2].loglog(hzaxis, bp_signalpower)
  72
+
  73
+axarr[5,0].plot(corrL, 'r')
  74
+axarr[5,1].plot(corrM, 'r')
  75
+axarr[5,2].plot(corrR, 'r')
  76
+
  77
+axarr[6,0].plot(bp_tempLeft, 'r')
  78
+axarr[6,1].plot(bp_tempMiddle, 'r')
  79
+axarr[6,2].plot(bp_tempRight, 'r')
2  simpleOptFilterTest2.py
@@ -21,7 +21,7 @@ def getNoisePulse(pl):
21 21
   #now filter the white noise with an iir filter to make it more realistic
22 22
   order = 4
23 23
   (b,a) = sig.iirfilter(order,0.1, btype='lowpass')
24  
-  iir4 = KIIRFourthOrder(-a[1], -a[2], -a[3], -a[4], b[0], b[1], b[2], b[3], b[4])
  24
+  iir4 = KIIRFourthOrder(a[1], a[2], a[3], a[4], b[0], b[1], b[2], b[3], b[4])
25 25
   iir4.SetInputPulse(whitenoise)
26 26
   print 'filter white noise', iir4.RunProcess()
27 27
 
3  simpleOptFilterTest3.py
@@ -22,7 +22,7 @@ def getNoisePulse(pl):
22 22
   #now filter the white noise with an iir filter to make it more realistic
23 23
   order = 4
24 24
   (b,a) = sig.iirfilter(order,0.1, btype='lowpass')
25  
-  iir4 = KIIRFourthOrder(-a[1], -a[2], -a[3], -a[4], b[0], b[1], b[2], b[3], b[4])
  25
+  iir4 = KIIRFourthOrder(a[1], a[2], a[3], a[4], b[0], b[1], b[2], b[3], b[4])
26 26
   iir4.SetInputPulse(whitenoise)
27 27
   print 'filter white noise', iir4.RunProcess()
28 28
 
@@ -98,6 +98,7 @@ def main(*args):
98 98
 
99 99
   tr = {}
100 100
   
  101
+  tr['nyquist'] = 50000  #50 kHz is the nyquist frequncy since we have a sample rate of 100 kHz
101 102
   
102 103
   ### create the pulse template and save it to the return
103 104
   (template, python_template) = getTemplate()
113  testBandPassFilterAndCorrelation.py
@@ -19,8 +19,8 @@ def getNoisePulse(pl):
19 19
   
20 20
   #now filter the white noise with an iir filter to make it more realistic
21 21
   order = 4
22  
-  (b,a) = sig.iirfilter(order,0.1, btype='lowpass')
23  
-  iir4 = KIIRFourthOrder(-a[1], -a[2], -a[3], -a[4], b[0], b[1], b[2], b[3], b[4])
  22
+  (b,a) = sig.iirfilter(order,0.5, btype='lowpass')
  23
+  iir4 = KIIRFourthOrder(a[1], a[2], a[3], a[4], b[0], b[1], b[2], b[3], b[4])
24 24
   iir4.SetInputPulse(whitenoise)
25 25
   print 'filter white noise', iir4.RunProcess()
26 26
 
@@ -68,19 +68,56 @@ def getTemplate():
68 68
     
69 69
   return (t, pt)
70 70
 
71  
-def createSignal(pl, template, noisepulse, amplitude = 1., delay = 100):
  71
+def createSignal(pl, template, noisepulse, amplitude = 1., delay = 2095):
72 72
   signal = std.vector("double")()
73 73
   for i in range(pl):
74  
-    signal.push_back( amplitude*template[i-delay] + noisepulse[i])
  74
+    if i >= delay:
  75
+      signal.push_back( amplitude*template[i-delay] + noisepulse[i])
  76
+    else:
  77
+      signal.push_back(noisepulse[i])
75 78
     
76 79
   return signal
77 80
   
78 81
 
  82
+def shiftSignalLeft(signal, shift):
  83
+  newsignal = copy.copy(signal)
  84
+  for i in range(len(newsignal)):
  85
+    if i + shift < len(signal):
  86
+      newsignal[i] = signal[i+shift]
  87
+    else:
  88
+      newsignal[i] = 0
  89
+  return newsignal
  90
+  
  91
+def shiftSignalRight(signal, shift):
  92
+  newsignal = copy.copy(signal)
  93
+  for i in range(len(newsignal)):
  94
+    if i > shift:
  95
+      newsignal[i] = signal[i-shift]
  96
+    else:
  97
+      newsignal[i] = 0
  98
+  
  99
+  return newsignal
  100
+  
  101
+def addWindowFunction(signal, window, center):
  102
+  #this assumes that outside of the range of the window function, the window = 0
  103
+  #the middle position of the window function is aligned with the value 'center'
  104
+  firstPos = center - len(window)/2
  105
+  j = 0
  106
+  for i in range(len(signal)):
  107
+    if i < firstPos:
  108
+      signal[i] = 0
  109
+    elif i >= firstPos and i < firstPos + len(window):
  110
+      signal[i] = signal[i] * window[j]
  111
+      j+=1
  112
+    elif i >= firstPos + len(window):
  113
+      signal[i] = 0
  114
+      
79 115
 
80 116
 def main(*args):
81 117
 
82 118
   tr = {}  #the object that is returned by this function - used to analyze the results
83 119
   
  120
+  tr['nyquist'] = 50000  #50 kHz is the nyquist frequncy since we have a sample rate of 100 kHz
84 121
   
85 122
   ### create the pulse template and save it to the return
86 123
   (template, python_template) = getTemplate()
@@ -104,7 +141,7 @@ def main(*args):
104 141
     
105 142
   #create the signal by adding the template to the noise
106 143
   #save it to the return
107  
-  signal = createSignal(pulseLength, python_template, noisepulse, 10000.)
  144
+  signal = createSignal(pulseLength, python_template, noisepulse, 500.)
108 145
   signalpy = []
109 146
   for i in range(pulseLength):
110 147
     signalpy.append(signal[i])
@@ -113,10 +150,23 @@ def main(*args):
113 150
   tr['signalpower'] = calculatePower(signal) ### calculate the power 
114 151
   
115 152
   
  153
+  #add a windowing function
  154
+  x = sig.blackman(3200)
  155
+  tr['window'] = x.tolist()
  156
+  addWindowFunction(signal, x, 6556)
  157
+  addWindowFunction(template, x, np.argmin(template))
  158
+  
  159
+  tr['window_signal'] = []
  160
+  tr['window_template'] = []
  161
+  for i in range(len(signal)):
  162
+    tr['window_signal'].append(signal[i])
  163
+  
  164
+  for i in range(len(template)):
  165
+    tr['window_template'].append(template[i])
116 166
   #perform the filter
117 167
   
118  
-  (b,a) = sig.iirfilter(2,[0.01, 0.4])
119  
-  filter = KIIRFourthOrder(-a[1], -a[2], -a[3], -a[4], b[0], b[1], b[2], b[3], b[4])
  168
+  (b,a) = sig.iirfilter(2,[0.001, 0.01])
  169
+  filter = KIIRFourthOrder(a[1], a[2], a[3], a[4], b[0], b[1], b[2], b[3], b[4])
120 170
   
121 171
   #here's the filter's frequency response function
122 172
   tr['b'] = b
@@ -164,17 +214,62 @@ def main(*args):
164 214
   correlation = KCorrelation()
165 215
 
166 216
   #set the template as the response to the correlation function
167  
-  correlation.SetResponse(bp_template)
  217
+  bp_tempLeft = shiftSignalRight(bp_template, 2090)
  218
+  bp_tempMiddle = shiftSignalRight(bp_template, 2095)
  219
+  bp_tempRight = shiftSignalRight(bp_template, 2100)
  220
+  
  221
+  bp_py = []
  222
+  for i in range(len(bp_tempLeft)):
  223
+    bp_py.append(bp_tempLeft[i])
  224
+  tr['bp_tempLeft'] = bp_py
  225
+  
  226
+  bp_py = []
  227
+  for i in range(len(bp_tempMiddle)):
  228
+    bp_py.append(bp_tempMiddle[i])
  229
+  tr['bp_tempMiddle'] = bp_py
  230
+  
  231
+  bp_py = []
  232
+  for i in range(len(bp_tempRight)):
  233
+    bp_py.append(bp_tempRight[i])
  234
+  tr['bp_tempRight'] = bp_py
  235
+  
  236
+  
  237
+  correlation.SetResponse(bp_tempLeft)
168 238
   correlation.SetInputPulse(bp_signal)
169 239
   correlation.RunProcess()
  240
+  corrOut = std.vector("double")()
  241
+  corrOutpy = []
  242
+  for i in range(correlation.GetOutputPulseSize()):
  243
+    corrOut.push_back(correlation.GetOutputPulse()[i])
  244
+    corrOutpy.append(correlation.GetOutputPulse()[i])
  245
+    
  246
+  tr['corrLeft'] = corrOutpy
170 247
   
  248
+  
  249
+  correlation.SetResponse(bp_tempMiddle)
  250
+  correlation.SetInputPulse(bp_signal)
  251
+  correlation.RunProcess()
171 252
   corrOut = std.vector("double")()
172 253
   corrOutpy = []
173 254
   for i in range(correlation.GetOutputPulseSize()):
174 255
     corrOut.push_back(correlation.GetOutputPulse()[i])
175 256
     corrOutpy.append(correlation.GetOutputPulse()[i])
176 257
     
177  
-  tr['corr'] = corrOutpy
  258
+  tr['corrMiddle'] = corrOutpy
  259
+  
  260
+  
  261
+  correlation.SetResponse(bp_tempRight)
  262
+  correlation.SetInputPulse(bp_signal)
  263
+  correlation.RunProcess()
  264
+  corrOut = std.vector("double")()
  265
+  corrOutpy = []
  266
+  for i in range(correlation.GetOutputPulseSize()):
  267
+    corrOut.push_back(correlation.GetOutputPulse()[i])
  268
+    corrOutpy.append(correlation.GetOutputPulse()[i])
  269
+    
  270
+  tr['corrRight'] = corrOutpy
  271
+  
  272
+  
178 273
     
179 274
   return tr
180 275
 
2  testiir4.py
@@ -42,7 +42,7 @@ def main(*args):
42 42
   tr['a'] = a
43 43
   tr['filter order'] =  order
44 44
   
45  
-  filter = KIIRFourthOrder(-a[1], -a[2], -a[3], -a[4], b[0], b[1], b[2], b[3], b[4])
  45
+  filter = KIIRFourthOrder(a[1], a[2], a[3], a[4], b[0], b[1], b[2], b[3], b[4])
46 46
   filter.SetInputPulse(whitenoise)
47 47
   print 'run filter', filter.RunProcess()
48 48
   

0 notes on commit 3aa2720

Please sign in to comment.
Something went wrong with that request. Please try again.