Infrasound Technology Workshop Bermuda 2008 Enhancements to the

Infrasound Technology Workshop Bermuda, 2008 Enhancements to the CTBTO operational automatic infrasound processing system David Brown, Nicolas Brachet, and Ronan Le Bras International Data Centre Software Applications Section Preparatory Commission for the Comprehensive Nuclear-Test-Ban Treaty Organization Provisional Technical Secretariat Vienna International Centre P. O. Box 1200, A-1400 Vienna, Austria E-mail: David. Brown@ctbto. org E-mail: Nicolas. Brachet@ctbto. org

Overview Two enhancements to the IDC automatic infrasound processing system are being tested: q Station Noise characterization Ø purpose Ø requirements Ø method Ø testing Ø real data Ø still to be done q Infrasound amplitude determination Ø purpose Ø procedure Ø testing Ø still to be done Infrasound Technology Workshop, Bermuda November 2008 Page 2

Station Noise Characterization Purpose Ø Provide a method for users of IDC external products to determine station ‘noise’ levels. • will be useful in determining network detection capability Ø Provide an additional utility for internal PTS station performance monitoring. Requirements Ø Noise data for each station would need to reveal: • seasonal or monthly variations • day and night time variations Ø External/Internal Users need to access the binary data and plot files. Method Ø Employ Power Spectral Density methods • along the lines of the analysis of Bowman et al. • will run automatically on incoming data and update seasonal or monthly averages. Infrasound Technology Workshop, Bermuda November 2008 Page 3

Station Noise Characterization Choose 4 hourly periods for each station: Ø 03: 30 -04: 30, 09: 30 -10: 30, 15: 30 -16: 30, 21: 30 -22: 30 local time Choose a FFT window function: Ø References: • Harris 1978: On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform Proc. IEEE Vol 66, No. 1, 1978. • Heinzel 2002: https: //www. lisa. unihannover. de/? page=publikationen&sub=publikationen&type=publi&publitype=3&lang=en Ø Require moderate frequency resolution Ø Require good amplitude resolution Infrasound Technology Workshop, Bermuda November 2008 Page 4

Station Noise Characterization Windowing: Ø a significant part of the PSD measurement process no window equivalent to High Frequencies Infrasound Technology Workshop, Bermuda November 2008 Page 5

Station Noise Characterization Windowing: Ø non-periodic nature of the windowing process leads to spectral leakage Ø frequency content and spectral leakage controlled by the transform of the window function: ‘Dirac Comb’ spectral leakage Infrasound Technology Workshop, Bermuda November 2008 Page 6

Station Noise Characterization Windowing: Ø Ensure the side lobes of the Window Transform are below the noise floor of the sensor system (microbarometer + filter): • MB 2005: approx. 10 -7 Pa 2/Hz → Nutall 4 a window (Heinzel et al, 2002) • Chaparral : approx. 10 -9 Pa 2/Hz → Nutall 4 c window (Heinzel et al, 2002) Infrasound Technology Workshop, Bermuda November 2008 Page 7

Station Noise Characterization ‘Welch’ PSD method Ø Divide hourly data into overlapping 3 minute windows one hour … • without averaging, the variance of the PSD for a individual frequency picket is of the order of the PSD at that frequency. • 20 Hz sample-rate → 3600 samples or 1800 frequency pickets: • 0. 0000000 Hz→ DC • 0. 0015625 Hz → 640 sec period • 0. 0031250 Hz → 320 sec period • … • 9. 9984375 Hz Infrasound Technology Workshop, Bermuda November 2008 Page 8

Station Noise Characterization Ø Number of overlapping windows is determined by the window ROV (Relative Overlap Value, see Heinzel et al, 2002), indicated in the following Table: window ROV % 3 -minute windows per hour ENBW (Equivalent Noise Band Width, bins) side lobe level emax % 3 db width (bins) (db) high frequency resolution windows Rectangular 0. 0 20 1. 0000 -13. 3 -36. 3 0. 9 general purpose windows Welch 29. 3 28 1. 2000 -21. 3 -22. 6 1. 2 Hanning 50. 0 39 1. 5000 -31. 5 -15. 1 1. 4 Hamming 50. 0 39 1. 3628 -42. 7 -18. 2 1. 3 Nutall 4 a 68. 0 60 2. 1253 -82. 6 -8. 1 2. 0 Nutall 4 c 65. 6 56 1. 9761 -98. 1 -9. 3 1. 9 high amplitude resolution windows Infrasound Technology Workshop, Bermuda FTSRS 75. 4 78 3. 7702 -76. 6 -0. 2 3. 7 HFT 248 D 84. 1 120 5. 6512 -284. 4 +0. 0 5. 6 November 2008 Page 9

Station Noise Characterization Ø Remove linear trend from each 3 minute data segment • otherwise DC component may generate spectral leakage Ø Multiply data segment with window function • remember to divide by the processing gain, ie, multiply by Ø Apply FFT Ø Determine PSD at frequency picket i according to the rule: N = samples; D = sample frequency; b = ENBW Accommodates the accumulation ability of the main lobe for white noise Ø average the PSD for all sub windows and take the log of the average • could consider taking the average of the log’s, but Bowman et al’s work suggests not log-normal distribution. Infrasound Technology Workshop, Bermuda November 2008 Page 10

Station Noise Characterization Testing Ø use artificial digitizer noise • Digitizing process generates (white) noise power spectral density: (see Lyons, 1997; Heinzel, 2002) is the Pa corresponding to one least-significant bit is the sampling frequency, Hz. Ø Example: with Infrasound Technology Workshop, Bermuda November 2008 Page 11

Station Noise Characterization Testing Ø use 4 different windows: q Rectangular (no window) q Hanning q Nutall 4 a (Heinzel 2002) q HFT 248 D (Heinzel 2002) -8. 079 Infrasound Technology Workshop, Bermuda November 2008 Page 12

Station Noise Characterization Testing Ø have initiated a series of ‘blind-tests’ with Lars Ceranna at BGR to assess the validity of the method. • using randomly chosen IMS array data • IS 07, IS 27, IS 32, IS 55 • 2007 -2008, 1 st day of the month 4 times per day. • will ‘manipulate’ the PSD methods until we get the same answers. Ø compare with known published results: Bowman et al. Infrasound Technology Workshop, Bermuda November 2008 Page 13

Station Noise Characterization Real Data Ø PSD noise data is being determined automatically for every IMS infrasound station. Ø Two outputs are generated: • Output 1: log-average PSD data for the given hours • binary data file + plot (. ps) • should be free of contamination from spectral leakage • have complete hourly and monthly-averaged PSD information for each station for September 2008 , October 2008 and the current part of November 2008. • for each 8 -element station the data growth rate is : 4 x (80000 + 130000 ) bytes = 840 Kb (. jpg Infrasound Technology Workshop, Bermuda binary) November 2008 Page 14

Station Noise Characterization Real Data Ø PSD noise data is being determined automatically for every IMS infrasound station. Ø Two outputs are generated: • Output 1: log-average PSD data for the given hours • binary data file + plot (. ps) • should be free of contamination from spectral leakage • have complete hourly and monthly-averaged PSD information for each station for September 2008 , October 2008 and the current part of November 2008. • for each 8 -element station the data growth rate is: 4 x (80000 + 130000 ) bytes = 520 Kb / day (approx. 5 Gb per year for 60 -stations) (generate graphic on the fly) Infrasound Technology Workshop, Bermuda November 2008 Page 15

Station Noise Characterization Real Data Ø PSD noise data is being determined automatically for every IMS infrasound station. Ø Two outputs are generated: • Output 2: monthly-average PSD data + variance smoothed using an 11 -point 6 -order Savitzky. Golay Filter for the specified times. • Savitzky-Golay filter o polynomial regression that preserves features of the distribution such as relative maxima, minima and width • have complete hourly and monthly-averaged PSD information for each station for September 2008 , October 2008 and the current part of November 2008. Infrasound Technology Workshop, Bermuda November 2008 Page 16

Station Noise Characterization Still to be done q incorporate as part of mainstream IDC processing q make data accessible to external users via a web-service interface q incorporate into Station Performance Monitoring tool for internal PTS use Infrasound Technology Workshop, Bermuda November 2008 Page 17

Infrasound Amplitude Determination Purpose Ø Provide a method for determining infrasound amplitudes in Pascals • useful information to write in the Bulletin • may be useful in Network processing Infrasound Technology Workshop, Bermuda November 2008 Page 18

Infrasound Amplitude Determination Procedure Ø Apply Window Function in the time domain • Use FTSRS ‘Flat Top’window • process sufficient data such that the window function will not drop by more than 50% from its peak value over the duration of the signal segment. • otherwise signal will get mangled by the window division processes window function signal 0. 5 extra data 0. 0 signal start time Infrasound Technology Workshop, Bermuda signal end time November 2008 Page 19

Infrasound Amplitude Determination Procedure Ø Apply FIR filter • two filter bands chosen: • 0. 5 to 3. 0 Hz • minfreq and maxfreq determined by pmcc detection algorithm Ø Remove window function from signal segment by dividing Ø Three amplitude measures applied to a time-aligned beam : • beam aligned according to azimuth + slowness in CSS detection table • Peak-to-Peak: used in several attenuation laws • LANL • DASE • RMS (user specifiable window length) • Amplitude of the Analytic Trace (see Olson, 2000) Ø For each detection: • 6 amplitude measures are being determined • the Peak-to-Peak 0. 5 -3. 0 Hz amplitude is being written to the CSS arrival table, all others are being written to the amplitude table. Infrasound Technology Workshop, Bermuda November 2008 Page 20

Infrasound Amplitude Determination Testing Ø A sinc function is used to test beamformer, filtering procedure and amplitude measurement • Peak-to-Peak (P-P) amplitude = 2. 434 ; Analytic amplitude = 2. 0 Nutall 4 a window 0. 001 - 9. 999 Hz 0. 15 - 0. 60 Hz Infrasound Technology Workshop, Bermuda Overlapping traces: • beamformer is working correctly • filter process is not generating a phase-shift • window division is not mangling the data non-overlapping traces: • filter process is not generating a phase-shift • window division is not mangling the data November 2008 Page 21

Infrasound Amplitude Determination Testing window high frequency resolution rectangle welch general purpose hanning Results: • for all windows the P-P amplitude measure 0. 001 to 9. 999 Hz 0. 15 to 0. 6 Hz measure provides the P-P 2. 434 2. 104 correct amplitude in the ‘unfiltered’ case RMS 1. 623 0. 900 • the HFT 248 D window Analytic 1. 987 1. 215 provides the best analytic P-P 2. 434 2. 104 amplitude in the unfiltered RMS 1. 618 0. 900 case. • all windows provided Analytic 1. 981 1. 215 reduced amplitudes for a 2 P-P 2. 434 2. 104 -octave band-pass filter RMS 1. 623 0. 900 centred on the dominant Analytic 1. 986 1. 215 frequency. P-P 2. 434 2. 104 RMS 1. 628 0. 900 Analytic 1. 991 1. 215 P-P 2. 434 2. 104 RMS 1. 632 0. 900 Analytic 1. 996 1. 215 nutall 4 a high amplitude resolution HFT 248 D Infrasound Technology Workshop, Bermuda November 2008 Page 22

Infrasound Amplitude Determination Still to be done q evaluate results Ø compare with analyst measured amplitudes q investigate utility of different amplitude measures Ø Analytic amplitude as a microbarom classifier ? Ø the RMS amplitude as a signal significance measure ? q incorporate as part of mainstream IDC processing Infrasound Technology Workshop, Bermuda November 2008 Page 23

Conclusions Summary q Two enhancements to the automatic infrasound processing system are currently under development: Ø station noise characterization Ø amplitude determination q Should be operational in around 3 -4 months. Infrasound Technology Workshop, Bermuda November 2008 Page 24