// Wavelet Analysis Tool // universal data container for wavelet transforms // used with DMT and ROOT // //$Id: wseries.hh,v 1.1 2005/05/16 06:05:10 igor Exp $ #ifndef WSERIES_HH #define WSERIES_HH #ifndef WAVEARRAY_HH #include "wavearray.hh" #endif #include #include #include "WaveDWT.hh" #include #include "wavecomplex.hh" #include "TNamed.h" typedef wavecomplex d_complex; typedef std::vector vector_int; template class WSeries : public wavearray { public: // constructors //: Default constructor WSeries(); //: Construct WSeries for specific wavelet type //+ default constructor explicit WSeries(const Wavelet &w); //: Construct from wavearray //!param: value - data to initialize the WSeries object explicit WSeries(const wavearray& value, const Wavelet &w); //: Copy constructor //!param: value - object to copy from WSeries(const WSeries& value); //: destructor virtual ~WSeries(); // operators WSeries& operator= (const wavearray &); WSeries& operator= (const WSeries &); WSeries& operator= (const DataType_t); // operator[](const slice &) sets the Slice object of the wavearray class virtual WSeries& operator[](const std::slice &); //: operators for WSeries objects, which can have different length //: it is required they have the same type of transform (standard or binary) //: and the same size of approximation levels. //: warning: there is no check that approximation levels have the same //: sampling rate. virtual WSeries& operator+=(WSeries &); virtual WSeries& operator-=(WSeries &); virtual WSeries& operator*=(WSeries &); // just to trick ANSI standard virtual WSeries& operator+=(wavearray &); virtual WSeries& operator-=(wavearray &); virtual WSeries& operator*=(wavearray &); virtual WSeries& operator+=(const DataType_t); virtual WSeries& operator-=(const DataType_t); virtual WSeries& operator*=(const DataType_t); // multiply layer by layer this and input wseries void mul(WSeries &); //: Dump data array to an ASCII file virtual void Dump(const char*, int=0); // accessors //: Get maximum possible level of decompostion int getMaxLevel(); //: Get level of decompostion inline int getLevel(){ return pWavelet->m_Level; } //: set level of decompostion inline void setLevel(size_t n){ pWavelet->m_Level = n; } //: Set black pixel probability inline void setbpp(double f) { bpp=f; return; } //: Get black pixel probability inline double getbpp() const { return bpp; } //: Set wavelet rate (0-level) inline void wrate(double r) {wRate=r; return;} //: Get wavelet rate inline double wrate() const { return wRate; } //: Set low frequency boundary inline void setlow(double f) { f_low=f>0. ? f : 0.; return; } //: Get low frequency boundary inline double getlow() const { return f_low; } //: Set high frequency boundary inline void sethigh(double f) { f_high = f; return; } //: get high frequency boundary inline double gethigh() const { return f_high; } //: Get max layer of decompostion inline int maxLayer(){ return pWavelet->BinaryTree() ? (1<getSlice(0).size();} // maximum sample index in 00 phase inline size_t maxIndex(){return sizeZero()*(maxLayer()+1)-1;} //: Get slice structure for specified layer inline std::slice getSlice(double n){ return pWavelet->getSlice(n); } //: Get wavelet layer frequency resolution inline double resolution(size_t n=0){ return frequency(1)-frequency(0); } //: Get central frequency of a layer // l - layer number. // TF map binary dyadic WDM // zero layer Fc dF/2 Fo 0 // non-zero Fc +n*dF Fn +n*dF double frequency(int l); //: Get layer index for input frequency // f - frequency Hz // TF map binary dyadic WDM // zero layer Fc dF/2 Fo 0 // non-zero Fc +n*dF Fn +n*dF int layer(double f); //: Extract wavelet coefficients from specified layer //!param: n - layer number int getLayer(wavearray &w, double n); //: replace wavelet data for specified layer with data from Sequence //!param: n - layer number void putLayer(wavearray &, double n); // extract wavelet amplitude for layer m and wavelet time index n // or sample(n,m) // !param: n - wavelet time index // !param: m - layer number inline DataType_t getSample(int n, double m) { std::slice S = this->getSlice(m); return (ndata[n*S.stride()+S.start()] : 0; } // extract wavelet amplitude for layer m and wavelet time index n // or sample(n,m) // !param: a - new amplitude // !param: n - wavelet time index // !param: m - layer number inline DataType_t putSample(DataType_t a, int n, double m) { std::slice S = this->getSlice(m); if(ndata[n*S.stride()+S.start()]=a; } // mutators virtual void resize(unsigned int); virtual void resample(double, int=6); //: initialize wavelet parameters from Wavelet object void setWavelet(const Wavelet &w); //: return true if WDM transform bool isWDM() {return pWavelet->m_WaveType==WDMT ? true : false;} //: Perform n steps of forward wavelet transform //!param: wavelet - n is number of steps (-1 means full decomposition) // WDM - n=-1 - otrhonormal, n=0 - power, n>0 - upsampled power void Forward(int n = -1); void Forward(wavearray &, int n = -1); void Forward(wavearray &, Wavelet &, int n = -1); //: Perform n steps of inverse wavelet transform //!param: n - number of steps (-1 means full reconstruction) void Inverse(int n = -1); // bandpass data and store in TF domain // ts - input time series // flow - low frequence boundary // fhigh - high frequency boundary // n - decomposition parameter void bandpass(wavearray &ts, double flow, double fhigh, int n=-1); // set wseries coefficients to a in layers between // flow - low frequence boundary // fhigh - high frequency boundary // a - value void bandpass(double flow, double fhigh, double a=0.); // maxEnergy: put maximum energy of delayed samples in this // param: wavearray - input time series // param: wavelet - wavelet used for the transformation // param: double - range of time delays // param: int - downsample factor to obtain coarse TD steps size_t maxEnergy(wavearray&ts, Wavelet &w, double t, int n=1); // produce wavescan from two up-sampled WDM TF data // par1 - high-frequency resolution WDM // par2 - low-frequency resolution WDM void scan(WSeries &, int=-1); // void scan(WSeries &, WSeries &); //++++++++++++++ wavelet data conditioning +++++++++++++++++++++++++++ //: calculate running medians with window of t seconds //: and subtract the median from this (false key) or //: calculate median for abs(this) and normalize this (true) virtual void median(double t, bool norm=false); //: apply linear predictor to each layer. // param: filter length in seconds // param: filter mode: -1/0/1 - backward/symmetric/forward // param: stride for filter training (0 - train on whole TS) // param: boundary offset to account for wavelet artifacts virtual void lprFilter(double,int=0,double=0.,double=0.); //: tracking of noise non-stationarity and whitening. // param 1 - time window dT. if = 0 - dT=T, where T is wavearray duration // param 2 - mode: 0 - no whitening, 1 - single whitening, >1 - double whitening // mode <0 - whitening using guadrature (WDM wavelet only) // param 3 - boundary offset // param 4 - noise sampling interval (window stride) // the number of measurements is k=int((T-2*offset)/stride) // if stride=0, then stride is set to dT // return: noise array if param2>0, median if param2=0 //!what it does: each wavelet layer is devided into k intervals. //!The data for each interval is sorted and the following parameters //!are calculated: median and the amplitude //!corresponding to 31% percentile (wp). Wavelet amplitudes (w) are //!normalized as w' = (w-median(t))/wp(t), where median(t) and wp(t) //!is a linear interpolation between (median,wp) measurements for //!each interval. virtual WSeries white(double,int,double=0.,double=0.); //: whiten TF map by using input noise array virtual bool white(WSeries ws, int mode=0); //: local whitening, works only for binary wavelets. //: returns array of noise rms for wavelet layers //!param: n - number of decomposition steps //!algorithm: //! 1) do forward wavelet transform with n decomposition steps //! 2) whiten wavelet layers and calculate noise rms as //! 1/Sum(1/var) //! 3) do inverse wavelet transform with n reconstruction steps virtual wavearray filter(size_t); //: works only for binary wavelets. //: calculates, corrects and returns noise variability //!param: first - time window to calculate normalization constants //! second - low frequency boundary for correction //! third - high frequency boundary for correction //!algorithm: //! 1) sort wavelet amplitudes with the same time stamp //! 2) calculate left(p) and right(p) amplitudes //! put (right(p)-left(p))/2 into output array //! 3) if first parameter >0 - devide WSeries by average //! variability virtual WSeries variability(double=0., double=-1., double=-1.); //: Selection of a fixed fraction of pixels //: reduced wavelet amplitudes are stored in this //: Returns fraction of non-zero coefficients. //!param: t - sub interval duration. If can not divide on integer // number of sub-intervals then add leftover to the last // one. //!param: f - black pixel fraction //!param: m - mode //!options: f = 0, m = 0 - returns black pixel occupancy //! m = 1 - set threshold f //! m = 2 - random policy //! m = 0 - random pixel selection virtual double fraction(double=0.,double=0.,int=0); //: calculate rank logarithic significance of wavelet pixels //: reduced wavelet amplitudes are stored in this //: Returns pixel occupancy for significance>0. //!param: n - sub-interval duration in seconds //!param: f - black pixel fraction //!options: f = 0 - returns black pixel occupancy virtual double significance(double, double=1.); //: calculate running rank logarithic significance of wavelet pixels //: reduced wavelet amplitudes are stored in this //: Returns pixel occupancy for significance>0. //!param: n - sub-interval duration in domain units //!param: f - black pixel fraction //!options: f = 0 - returns black pixel occupancy virtual double rsignificance(size_t=0, double=1.); //: calculate running rank logarithic significance of wavelet pixels //: reduced wavelet amplitudes are stored in this //: Returns pixel occupancy for significance>0. //!param: T - sliding window duration in seconds //!param: f - black pixel fraction //!param: t - sliding step in seconds //!options: f = 0 - returns black pixel occupancy //!options: t = 0 - sliding step = wavelet time resolution. virtual double rSignificance(double, double=1., double=0.); //: calculate running logarithic significance of wavelet pixels //: reduced wavelet amplitudes are stored in this //: Returns pixel occupancy for significance>0. //!param: T - sliding window duration in seconds //!param: f - black pixel fraction //!param: t - sliding step in seconds //!options: f = 0 - returns black pixel occupancy //!options: t = 0 - sliding step = wavelet time resolution. virtual double gSignificance(double, double=1., double=0.); //: Selection of a fixed fraction of pixels for each wavelet layer //: Returns fraction of non-zero coefficients. //!param: f - black pixel fraction //!param: m - mode //!options: f = 0 - returns black pixel occupancy //! m = 1 - set threshold f, returns percentile amplitudes //! m =-1 - set threshold f, returns wavelet amplitudes //! m > 1 - random policy,returns percentile amplitudes //! m <-1 - random policy,returns wavelet amplitudes //! m = 0 - random pixel selection //! if m<0 return wavelet amplitudes instead of the percentile amplitude virtual double percentile(double=0.,int=0, WSeries* = NULL); //: clean up single pixels //!param: S - threshold on pixel significance //!return pixel occupancy. virtual double pixclean(double=0.); //: select pixels from *this which satisfy a coincidence rule //: within specified window //!param: WSeries object used for coincidence //!param: coincidence window in seconds //!return pixel occupancy virtual double coincidence(WSeries &, int=0, int=0, double=0.); //: select pixels from *this which satisfy a coincidence rule //: within specified window w, above threshold T, //!param: WSeries object used for coincidence //!param: coincidence window in seconds //!param: threshold on significance //!return pixel occupancy virtual double Coincidence(WSeries &, double=0., double=0.); //: calculate calibration coefficients and apply energy calibration //: to wavelet data in this //: AS_Q calibration: R(f)=(1 + gamma*(R*C-1))/alpha*C //: DARM calibration: R(f)=(1 + gamma*(R*C-1))/gamma*C //: input //!param: number of samples in calibration arrays R & C //!param: frequency resolution //!param: pointer to response function R in Fourier domain //!param: pointer to sensing function C in Fourier domain //!param: time dependent calibration coefficient alpha //!param: time dependent calibration coefficient gamma //!param: 0/1 - AS_Q/DARM_ERR calibration, by default is 0 //!return array with calibration constants for each wavelet layer virtual WSeries calibrate(size_t, double, d_complex*, d_complex*, wavearray &, wavearray &, size_t ch=0); //: mask WSeries data with the pixel mask defined in wavecluster object //!param: int n //! if n<0, zero pixels defined in mask (regression) //! if n>=0, zero all pixels except ones defined in the mask //!param: bool - if true, set WSeries data to be positive //! if pMask.size()=0, mask(0,true) is equivalent to abs(data) //!return core pixel occupancy //virtual double mask(int=1, bool=false); //: calculate pixel occupancy for clusters passed selection cuts //!param: true - core pccupancy, false - total occupancy; //!return wavearray with occupancy. //virtual wavearray occupancy(bool=true); // print wseries parameters void print(); // *MENU* virtual void Browse(TBrowser *b) {print();} // data members //: parameters of wavelet transform WaveDWT* pWavelet; //: whitening mode ULong_t w_mode; //: black pixel probability double bpp; //: wavelet zero layer rate double wRate; //: low frequency boundary double f_low; //: high frequency boundary double f_high; ClassDef(WSeries,1) }; // class WSeries #endif // WSERIES_HH