- Email: [email protected]
Matched DDWT ROI Compression Engine for the Imaging Particle Detector
Matched DDWT ROI Compression Engine for the Imaging Particle Detector Wojciech Półchłopek, Roman Rumian Institute of Electronics, AGH University of Science and Technology, al. Mickiewicza 30, 30-059 Cracow, Poland phone: + (48) 12 617 27 00, fax: + (48) 12 617 30 45, email: [email protected], [email protected] Abstract: A new multirate DDWT (Double Density Wavelet Transform) matched filter ROI (region of interest) processor for the ICARUS (Imaging Cosmic And Rare Underground Signals) particle detector [Amerio(2004)] has been designed. The ROI extraction engine is based on matched and unmatched Wiener filtering using coupled DWT and DDWT processing. High through-put image-like ICARUS detector data (160MBps data rate per crate) is able to be compressed 600 times using multi-stage compression with ROI extraction. The 32-channel (160MBps data rate) processor has been fully implemented in a low cost FPGA device thanks to Fast Integer Arithmetic Wavelet Transform (FIAWT) [Półchłopek(2006)] algorithm implementation. Keywords: signal detection, signal processing algorithms, Wiener filters, spectral transformations, discrete-time systems 1. INTRODUCTION The T600 ICARUS Detector [Amerio(2004), Ankowski(2010)] is the large scale (dual 300 ton liquid Argon) electronic bubble-chamber with the total data throughput of several GB per second. The signals (see fig. 1) gathered in the chamber are of low S/N ratio because of the special signal conditions (charge collected in a very large volume form the ionized liquid Argon caused by the low energy particles as neutrinos). The important data occupies as low as 0.5% of the total collected volume. The typical low amplitude signal is shown in figure 1. The typical ICARUS signals can be fitted with 5 parameter function f(t) [Amerio(2004)]:
Amplitude
e
f (t ) = B + A
The time domain signal detector DAEDALUS [Arneodo et al. (1998)] designed for the 50l prototype was unable to handle with low SNR (close to 0dB) signals from 600 ton ICARUS module (T600), because was designed to operate with at least 20dB SNR signals. The only solution was to design time-frequency signal detector with time resolution at least of 32 samples which was the maximum decimation factor due to ROI window of 64 samples. 2. T600 SIGNAL AND NOISE ANALYSIS (WIENER SOLUTION) The ICARUS detector test runs revealed much lower than expected SNR of the signals and presence of various stationary and non-stationary noises and interferences. The noise analysis was presented in [Gibin(2003)] is shown in figure 2 (power spectral density function of signals were obtained with fitting procedure of detected signals with function (1.1)). The analysis will lead to the statistically optimal Wiener solution using PSD of signals Ss(f) and noise Sn(f) [Wiener(1950)],[ Miller(2004)]: 2
H( f ) =
− ( t − t0 ) τ 1 +τ 2
Ss ( f ) Ss ( f ) + Sn ( f )
(2.1)
− ( t −t 0 )
1+ e
τ2
which could be well fitted with exponent function of three parameters [10]: Time
⎛ f ⎞ ⎜ 1 c ⎟ H fit ( f ) = ⎜ e f −b ⎟ ⎜ ⎟ ⎝1 + e a ⎠ 2
Fig. 1. Typical signal from ICARUS detector (green) and fitted signal (red) f (t ) = B + A
e
− ( t − t0 ) τ 1 +τ 2 − ( t − t0 )
1+ e
τ2
(1.1)
2
(2.2)
The fitting function was chosen because of CDF [Cohen(1992), Donoho(1992)] spectral response similarity. The Hfit(f) function with parameters a=0.33·b; c=2.5·b well approximates CDF(2,2) spectral response (see fig. 2). The
procedure assumes linear phase filters, which can be used in lifting application biorthogonal wavelets transform CDF.
computationally complex and thus impossible to implement in the real-time ICARUS DAQ system.
H ( e jω )
4. DDWT DETECTION BASED ON WIENER SOLUTION FILTERS
f [Hz] x 102
Fig. 2. Wiener filter frequency response for ICARUS signals (eq. 1.3) – green, fitted with function Hfit(f) [Gibin(2003)] and CDF(2,2) synthesis wavelet (scale 32) amplitude spectrum.
3. WAVELET SCHEME DETECTION OF THE SIGNALS The procedure of detection described in [Batko (1999)] which uses Continuous Wavelet Transform can be well adapted to the ICARUS signals. The Continuous Wavelet Transform: +∞
Ca , b =
∫
−∞
⎛t −b⎞ f (t ) ⋅ ψ ⎜ ⎟ dt ⎝ a ⎠
(3.1)
of the signal is checked for local maxima and found maximum index bmax produces one dimensional function (see fig. 3) with one parameter bmax. The shape of this function is then checked for correlation with noise free signal shapes. This procedure uses CWT approximation which is 300
250
Ca,bmax 200
The real-time solution is possibile only with DDWT approximation of CWT. The shape of Ca,bmax function is then reduced to a three point approximation DWT (and for better time resolution DDWT). The DDWT approximation of CWT is efficient only when the middle filter is matched (scale 32 see figures 2 and 3) with optimal Wiener filter, and the two other filters are unmatched (scales 16 and 64 - see figure 3). This solution need six stages of DDWT decomposition which has been applied using DDWT FIAWT [Półchłopek (2002a),(2002b),(2006)] algorithm. The matched filter output (details on 5-th stage of decomposition) is then followed by over-threshold checking procedure. The two unmatched DDWT outputs are compared with dynamic threshold (their amplitude must be lower than matched filter output). 5. ROI ENCODING AND COMPRESSION RESULTS The initial state simulations were done in the Matlab® environment. Compression results were obtained via C++ data analysis program called Qscan with compression and ROI extraction interface. This application is designed for processing the ICARUS raw data files. Table 1. ROI extraction and compression results comparison. Time analysis ROI blocks detected (left 16017,15 chamber – low SNR) /2,80 total/per channel ROI blocks detected (right 5257,55 chamber – better SNR) /0,92 total/per channel Compression ratio (64 34,46 samples window) Compression ratio with additional lossless FIAWT 172 wavelet compression Compression ratio with additional high quality 344 lossy FIAWT wavelet compression False detection ratio 67,18% (left chamber - low SNR)
DDWT DDWT unmatched matched 5526,83 /0,96
5918,07 /1,03
4768,64 /0,83
5885,82 /1,03
71,21
62,11
356
310
712
620
13,72%
0,54%
150
100
50
0
0
20
40
60
80
100
120
140
a (scale)
Fig. 3. Maximum (bmax) coefficients of CWT.
Final ROI extraction and compression results are shown in the table 1. DDWT matched filter procedure shows high robustness to the noise and low false detection grade. Two dimensional image-like views (fig. 4) show algorithms sensitivities with small amplitude and low SNR (0dB and less) and confirms the best matched DDWT algorithm performance.
Figure 4. Signal recognition (ROI extraction trigger) results: time domain algorithm (a), DDWT matched filters algorithm (b) and DDWT algorithm (filters not matched) (c) RWB
S TROBE ADDR THRS
RESET SYNC_IN CLOCK
DATA_IN<9:0>
CONTROL LOGIC
DEMUX 1:32
PARAMETER REGISTERS
3 stage DWT processor ch0
3 stage DDWT processor ch0
3 stage DWT processor ch1
3 stage DDWT processor ch1
. . . . 3 stage DWT processor Ch31
. .
Windowing Block DWT and MUX 32:1
3 stage DDWT processor Ch31 Threshold and TRIGGER LOGIC
Figure 5. THESEUS FPGA module block diagram.
DATA_OUT<9:0>
OFL PEAKM<31:0>
6. HARDWARE IMPLEMENTATION The real-time ROI extraction engine application is based on DWT and DDWT (reverse CDF(2,2) wavelet) processor and threshold comparison block using oversampled FIAWT (according to a’trous algorithm) applied into the FPGA. The block diagram of the module THESEUS is shown in figure 5. The DWT and DDWT modules operate in a synchro-nous way, and fully implements FIAWT algorithm [Półchłopek (2006,2007)]. The synthesized project was implemented into the spe-cific FPGA device - Xilinx XC3S1500FG456-4 (speed grade 4 3,328CLBs). The final FPGA design consists of three stage DWT (CDF (2,2)) and three simplified DDWT (CDF (2,2)) used for data recognition in each of the 32 data channels sampled at 2,5MHz (multiplexed into 2 data streams of 40 MHz each). Each DWT stage resulted in utilization of 54 FPGA slices and could operate at frequency up to 50MHz. 7. CONCLUSIONS The new fully integer processing of the matched DDWT ROI extraction/compression enables very fast application and thus it can be very useful in application in real-time systems. The most important property of this concept is the possibility of simple and fast application into FPGA or ASIC chip. Until now the online compression has been applied to the part of the DAQ system only for testing purposes. The final application is foreseen in a scattered multi-processing unit DAQ architecture in mixed DSP/Xilinx-FPGA system. 8. ACKNOWLEGEMENTS Special thanks for the ICARUS collaboration especially for Padova Group (Univesity of Padova, Physics Department) and Polish ICARUS Group (IFJ, Cracow). 9. REFERENCES Amerio S., [et al.], Półchłopek W., [et al.], (ICARUS Collaboration) (2004). “Design, construction and tests of the ICARUS T600 detector”, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Volume: 527, Issue: 3, July 21, pp. 329-410. Ankowski A., [et al.], W. Półchłopek, [et al.] (ICARUS Collaboration) (2010) „Energy reconstruction of
electromagnetic showers from π0 decays with the ICARUS T600 Liquid Argon TPC” Acta Physica Polonica B; ISSN 0587-4254. — 2010 vol. 41 nr 1 s. 103–125. Arneodo, F.; Benetti, ; [et. al.] (ICARUS Collaboration) (1998). “Performance evaluation of a hit finding algorithm for the ICARUS detector” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment Volume: 412, Issue: 2-3, August 1, pp. 440453. Batko W., Mikulski A., (1999). Wavelet Transform of Impulse Signals, Machine Dynamics Problems, vol. 23, no. 2, pp. 139-146, Warsaw. Cohen A., Daubechies I., Feauveau J. C., (1992). Biorthogonal bases of compactly supported wavelets, Com-mun. on Pure and Appl. Math., vol. 45, pp. 485–560. Donoho D. L., (1992). “Interpolating wavelet transforms” Preprint, Department of Statistics, Stanford Univer-sity. Gibin D., (2003). “Wire signals : a first look for their usage in a local trigger” ICARUS Granada Meeting internal note. Miller S. L., Childers D. G., (2004). Probability and Random Processes, Elsevier Academic Press. Półchłopek W., Ziółko M., (2002a). “Wavelet Transform Compression and Denoising in Real-Time System” Proceedings of CNDSP Conferrence, Stafford, pp. 141148, Półchłopek W., Ventura S., Pietropaolo F., (2002b). “Wavelet Transform Compression and Denoising in Real-Time System (Proposal for the ICARUS DAQ System)” ICARUS TM2002/12 , Padova - ICARUS collaboration internal note: for pdf copy write to author Półchłopek W., Maj W., Padee W., (2006). “Fast Integer Arithmetic Wavelet Transform. Properties and Application in FPGA/DSP System” Proceedings of EUSIPCO 2006 Conference, Florence. Półchłopek W., (2007). “Discrete Time Wavelet Transform Applications in Real-Time Systems”, PhD Thesis, AGH University of Science and Technology, Cracow. Wiener N., (1950). Extrapolation, Interpolation and Smoothing of Stationary Time Series, The Technology Press and Wiley, New York.
Amplitude
e
f (t ) = B + A
The time domain signal detector DAEDALUS [Arneodo et al. (1998)] designed for the 50l prototype was unable to handle with low SNR (close to 0dB) signals from 600 ton ICARUS module (T600), because was designed to operate with at least 20dB SNR signals. The only solution was to design time-frequency signal detector with time resolution at least of 32 samples which was the maximum decimation factor due to ROI window of 64 samples. 2. T600 SIGNAL AND NOISE ANALYSIS (WIENER SOLUTION) The ICARUS detector test runs revealed much lower than expected SNR of the signals and presence of various stationary and non-stationary noises and interferences. The noise analysis was presented in [Gibin(2003)] is shown in figure 2 (power spectral density function of signals were obtained with fitting procedure of detected signals with function (1.1)). The analysis will lead to the statistically optimal Wiener solution using PSD of signals Ss(f) and noise Sn(f) [Wiener(1950)],[ Miller(2004)]: 2
H( f ) =
− ( t − t0 ) τ 1 +τ 2
Ss ( f ) Ss ( f ) + Sn ( f )
(2.1)
− ( t −t 0 )
1+ e
τ2
which could be well fitted with exponent function of three parameters [10]: Time
⎛ f ⎞ ⎜ 1 c ⎟ H fit ( f ) = ⎜ e f −b ⎟ ⎜ ⎟ ⎝1 + e a ⎠ 2
Fig. 1. Typical signal from ICARUS detector (green) and fitted signal (red) f (t ) = B + A
e
− ( t − t0 ) τ 1 +τ 2 − ( t − t0 )
1+ e
τ2
(1.1)
2
(2.2)
The fitting function was chosen because of CDF [Cohen(1992), Donoho(1992)] spectral response similarity. The Hfit(f) function with parameters a=0.33·b; c=2.5·b well approximates CDF(2,2) spectral response (see fig. 2). The
procedure assumes linear phase filters, which can be used in lifting application biorthogonal wavelets transform CDF.
computationally complex and thus impossible to implement in the real-time ICARUS DAQ system.
H ( e jω )
4. DDWT DETECTION BASED ON WIENER SOLUTION FILTERS
f [Hz] x 102
Fig. 2. Wiener filter frequency response for ICARUS signals (eq. 1.3) – green, fitted with function Hfit(f) [Gibin(2003)] and CDF(2,2) synthesis wavelet (scale 32) amplitude spectrum.
3. WAVELET SCHEME DETECTION OF THE SIGNALS The procedure of detection described in [Batko (1999)] which uses Continuous Wavelet Transform can be well adapted to the ICARUS signals. The Continuous Wavelet Transform: +∞
Ca , b =
∫
−∞
⎛t −b⎞ f (t ) ⋅ ψ ⎜ ⎟ dt ⎝ a ⎠
(3.1)
of the signal is checked for local maxima and found maximum index bmax produces one dimensional function (see fig. 3) with one parameter bmax. The shape of this function is then checked for correlation with noise free signal shapes. This procedure uses CWT approximation which is 300
250
Ca,bmax 200
The real-time solution is possibile only with DDWT approximation of CWT. The shape of Ca,bmax function is then reduced to a three point approximation DWT (and for better time resolution DDWT). The DDWT approximation of CWT is efficient only when the middle filter is matched (scale 32 see figures 2 and 3) with optimal Wiener filter, and the two other filters are unmatched (scales 16 and 64 - see figure 3). This solution need six stages of DDWT decomposition which has been applied using DDWT FIAWT [Półchłopek (2002a),(2002b),(2006)] algorithm. The matched filter output (details on 5-th stage of decomposition) is then followed by over-threshold checking procedure. The two unmatched DDWT outputs are compared with dynamic threshold (their amplitude must be lower than matched filter output). 5. ROI ENCODING AND COMPRESSION RESULTS The initial state simulations were done in the Matlab® environment. Compression results were obtained via C++ data analysis program called Qscan with compression and ROI extraction interface. This application is designed for processing the ICARUS raw data files. Table 1. ROI extraction and compression results comparison. Time analysis ROI blocks detected (left 16017,15 chamber – low SNR) /2,80 total/per channel ROI blocks detected (right 5257,55 chamber – better SNR) /0,92 total/per channel Compression ratio (64 34,46 samples window) Compression ratio with additional lossless FIAWT 172 wavelet compression Compression ratio with additional high quality 344 lossy FIAWT wavelet compression False detection ratio 67,18% (left chamber - low SNR)
DDWT DDWT unmatched matched 5526,83 /0,96
5918,07 /1,03
4768,64 /0,83
5885,82 /1,03
71,21
62,11
356
310
712
620
13,72%
0,54%
150
100
50
0
0
20
40
60
80
100
120
140
a (scale)
Fig. 3. Maximum (bmax) coefficients of CWT.
Final ROI extraction and compression results are shown in the table 1. DDWT matched filter procedure shows high robustness to the noise and low false detection grade. Two dimensional image-like views (fig. 4) show algorithms sensitivities with small amplitude and low SNR (0dB and less) and confirms the best matched DDWT algorithm performance.
Figure 4. Signal recognition (ROI extraction trigger) results: time domain algorithm (a), DDWT matched filters algorithm (b) and DDWT algorithm (filters not matched) (c) RWB
S TROBE ADDR THRS
RESET SYNC_IN CLOCK
DATA_IN<9:0>
CONTROL LOGIC
DEMUX 1:32
PARAMETER REGISTERS
3 stage DWT processor ch0
3 stage DDWT processor ch0
3 stage DWT processor ch1
3 stage DDWT processor ch1
. . . . 3 stage DWT processor Ch31
. .
Windowing Block DWT and MUX 32:1
3 stage DDWT processor Ch31 Threshold and TRIGGER LOGIC
Figure 5. THESEUS FPGA module block diagram.
DATA_OUT<9:0>
OFL PEAKM<31:0>
6. HARDWARE IMPLEMENTATION The real-time ROI extraction engine application is based on DWT and DDWT (reverse CDF(2,2) wavelet) processor and threshold comparison block using oversampled FIAWT (according to a’trous algorithm) applied into the FPGA. The block diagram of the module THESEUS is shown in figure 5. The DWT and DDWT modules operate in a synchro-nous way, and fully implements FIAWT algorithm [Półchłopek (2006,2007)]. The synthesized project was implemented into the spe-cific FPGA device - Xilinx XC3S1500FG456-4 (speed grade 4 3,328CLBs). The final FPGA design consists of three stage DWT (CDF (2,2)) and three simplified DDWT (CDF (2,2)) used for data recognition in each of the 32 data channels sampled at 2,5MHz (multiplexed into 2 data streams of 40 MHz each). Each DWT stage resulted in utilization of 54 FPGA slices and could operate at frequency up to 50MHz. 7. CONCLUSIONS The new fully integer processing of the matched DDWT ROI extraction/compression enables very fast application and thus it can be very useful in application in real-time systems. The most important property of this concept is the possibility of simple and fast application into FPGA or ASIC chip. Until now the online compression has been applied to the part of the DAQ system only for testing purposes. The final application is foreseen in a scattered multi-processing unit DAQ architecture in mixed DSP/Xilinx-FPGA system. 8. ACKNOWLEGEMENTS Special thanks for the ICARUS collaboration especially for Padova Group (Univesity of Padova, Physics Department) and Polish ICARUS Group (IFJ, Cracow). 9. REFERENCES Amerio S., [et al.], Półchłopek W., [et al.], (ICARUS Collaboration) (2004). “Design, construction and tests of the ICARUS T600 detector”, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Volume: 527, Issue: 3, July 21, pp. 329-410. Ankowski A., [et al.], W. Półchłopek, [et al.] (ICARUS Collaboration) (2010) „Energy reconstruction of
electromagnetic showers from π0 decays with the ICARUS T600 Liquid Argon TPC” Acta Physica Polonica B; ISSN 0587-4254. — 2010 vol. 41 nr 1 s. 103–125. Arneodo, F.; Benetti, ; [et. al.] (ICARUS Collaboration) (1998). “Performance evaluation of a hit finding algorithm for the ICARUS detector” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment Volume: 412, Issue: 2-3, August 1, pp. 440453. Batko W., Mikulski A., (1999). Wavelet Transform of Impulse Signals, Machine Dynamics Problems, vol. 23, no. 2, pp. 139-146, Warsaw. Cohen A., Daubechies I., Feauveau J. C., (1992). Biorthogonal bases of compactly supported wavelets, Com-mun. on Pure and Appl. Math., vol. 45, pp. 485–560. Donoho D. L., (1992). “Interpolating wavelet transforms” Preprint, Department of Statistics, Stanford Univer-sity. Gibin D., (2003). “Wire signals : a first look for their usage in a local trigger” ICARUS Granada Meeting internal note. Miller S. L., Childers D. G., (2004). Probability and Random Processes, Elsevier Academic Press. Półchłopek W., Ziółko M., (2002a). “Wavelet Transform Compression and Denoising in Real-Time System” Proceedings of CNDSP Conferrence, Stafford, pp. 141148, Półchłopek W., Ventura S., Pietropaolo F., (2002b). “Wavelet Transform Compression and Denoising in Real-Time System (Proposal for the ICARUS DAQ System)” ICARUS TM2002/12 , Padova - ICARUS collaboration internal note: for pdf copy write to author Półchłopek W., Maj W., Padee W., (2006). “Fast Integer Arithmetic Wavelet Transform. Properties and Application in FPGA/DSP System” Proceedings of EUSIPCO 2006 Conference, Florence. Półchłopek W., (2007). “Discrete Time Wavelet Transform Applications in Real-Time Systems”, PhD Thesis, AGH University of Science and Technology, Cracow. Wiener N., (1950). Extrapolation, Interpolation and Smoothing of Stationary Time Series, The Technology Press and Wiley, New York.