- Timestamp:
- 09/24/14 12:50:19 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/pyscripts/new_pyfact/pyfact.py
r17972 r17973 8 8 import numpy as np 9 9 import ROOT 10 import pylab11 12 from scipy.signal import firwin, iirdesign, lfilter13 10 14 11 ########## BUILDING OF THE SHARED OBJECT FILES ############################### … … 248 245 f.next(3) 249 246 print "Event Id:", f.cols['EventNum'] 250 251 import matplotlib.pyplot as plt252 plt.ion()253 """254 for i in range(f.header['NPIX']):255 plt.cla()256 plt.title("Event %d"%(f.cols['EventNum']))257 plt.plot(f.cols['CalibData2D'][i], '.:', label='pixel %d'%i)258 plt.legend()259 answer = raw_input('anykey for next pixel; "q" to quit. :')260 if 'q' in answer.lower():261 break262 """263 264 265 266 from scipy import signal267 from scipy.signal import freqs, iirfilter268 #b, a = iirfilter(4, [1, 10], 1, 60, analog=True, ftype='cheby1')269 #b, a = signal.iirfilter(17, [50, 200], rs=60, btype='band',270 # analog=True, ftype='cheby2')271 #w, h = signal.freqs(b, a, 1000)272 #w, h = freqs(b, a, worN=np.logspace(-1, 2, 1000))273 274 b, a = iirfilter(N=4, # the order of the filter275 Wn = [0.15, 0.25], # critical frequencies276 rp = 1, # maximum ripple in passband277 rs = 60, # minimum attenuation in stopband278 btype='bandstop', # {'bandpass', 'lowpass', 'highpass', 'bandstop'}279 #analog=True,280 ftype='cheby1')281 #w, h = freqs(b, a, worN=np.logspace(-1, 5, 1000))282 w, h = freqs(b, a)283 plt.figure()284 plt.semilogx(w, 20 * np.log10(abs(h)))285 plt.xlabel('Frequency')286 plt.ylabel('Amplitude response [dB]')287 plt.grid()288 plt.show()289 290 291 292 293 294 fig1, (ax11,ax12) = plt.subplots(nrows=2)295 fig2, (ax21,ax22) = plt.subplots(nrows=2)296 297 flt = firwin(13, [0.15,0.25])298 t = np.arange(f.cols['CalibData2D'].shape[1])/2.*1e-9299 fig3, (ax31, ax32) = plt.subplots(nrows=2)300 ax31.plot(flt)301 #w, h = freqs(flt, [1.0])302 w, h = freqs(flt, [1.0], worN=np.logspace(-1, 5, 1000))303 ax32.semilogx(w, 20 * np.log10(abs(h)))304 ax32.set_xlabel('Frequency')305 ax32.set_ylabel('Amplitude response [dB]')306 ax32.grid()307 308 309 310 for i in range(f.header['NPIX']):311 data = f.cols['CalibData2D'][i]312 fil_data = lfilter(flt, [1.0], data)313 314 ax11.cla()315 ax12.cla()316 ax11.set_title("Event %d"%(f.cols['EventNum']))317 ax11.plot(t, data, '.:', label='pixel %d'%i)318 ax11.legend()319 ax11.set_ylim(-20,40)320 ax12.plot(t[:fil_data.size], fil_data, '.:', label='pixel %d'%i)321 ax12.legend()322 ax12.set_ylim(-20,40)323 324 """325 Pxx, freqs, t, plot = pylab.specgram( f.cols['CalibData2D'][i],326 NFFT=32, Fs=2000000000,327 detrend=pylab.detrend_none,328 window=pylab.window_hanning,329 noverlap=16)330 """331 ax21.cla()332 ax22.cla()333 dat, freqs, bins, im = ax21.specgram(data,334 NFFT=64, Fs=2000000000,335 detrend=pylab.detrend_none,336 window=pylab.window_hanning,337 noverlap=60,338 interpolation='nearest')339 ax21.axis('tight')340 ax21.set_ylabel("frequency [Hz]")341 ax21.set_xlabel("time [s]")342 343 dat2, freqs2, bins2, im2 = ax22.specgram(data,344 NFFT=64, Fs=2000000000,345 detrend=pylab.detrend_none,346 window=pylab.window_hanning,347 noverlap=60,348 interpolation='nearest')349 ax22.axis('tight')350 ax22.set_ylabel("frequency [Hz]")351 ax22.set_xlabel("time [s]")352 353 plt.show()354 355 answer = raw_input('anykey for next pixel; "q" to quit. :')356 if 'q' in answer.lower():357 break
Note:
See TracChangeset
for help on using the changeset viewer.