00001 /*************************************************************************** 00002 ** 00003 ** This file is part of GeopsyCore. 00004 ** 00005 ** This library is free software; you can redistribute it and/or 00006 ** modify it under the terms of the GNU Lesser General Public 00007 ** License as published by the Free Software Foundation; either 00008 ** version 2.1 of the License, or (at your option) any later version. 00009 ** 00010 ** This file is distributed in the hope that it will be useful, but WITHOUT 00011 ** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 ** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 00013 ** License for more details. 00014 ** 00015 ** You should have received a copy of the GNU Lesser General Public 00016 ** License along with this library; if not, write to the Free Software 00017 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00018 ** 00019 ** See http://www.geopsy.org for more information. 00020 ** 00021 ** Created : 2003-11-22 00022 ** Authors : 00023 ** Marc Wathelet 00024 ** Marc Wathelet (ULg, Liège, Belgium) 00025 ** Marc Wathelet (LGIT, Grenoble, France) 00026 ** 00027 ***************************************************************************/ 00028 00029 #ifndef DOUBLESIGNAL_H 00030 #define DOUBLESIGNAL_H 00031 00032 #include "GeopsyCoreDLLExport.h" 00033 #include "SignalTemplate.h" 00034 #include "ComplexSignal.h" 00035 #include "WindowingParameters.h" 00036 00037 namespace GeopsyCore { 00038 00039 class KeepSignal; 00040 class MorletParameters; 00041 00042 class GEOPSYCORE_EXPORT DoubleSignal: public SignalTemplate<double> 00043 { 00044 TRANSLATIONS("DoubleSignal") 00045 public: 00046 enum SignalType { 00047 Waveform=0, 00048 Spectrum=1, 00049 Phase=2, 00050 RealSpectrum=3, // For real spectrum and their fft's 00051 CAmpWaveform=4, 00052 CPhaseWaveform=5, 00053 ArrivalTime=6, 00054 UndefinedSignalType=7 00055 }; 00056 00057 DoubleSignal( ); 00058 DoubleSignal(int n); 00059 DoubleSignal(const DoubleSignal& o); 00060 00061 void lockProperties() const {_propertiesMutex.lockForWrite();} 00062 void constLockProperties() const {_propertiesMutex.lockForRead();} 00063 void unlockProperies() const {_propertiesMutex.unlock();} 00064 00065 // Properties 00066 double deltaT() const {return _deltaT;} 00067 double deltaF() const {return _deltaT;} 00068 virtual void setDeltaT(double newval) {_deltaT=newval;} 00069 void setDeltaF(double newval) {setDeltaT(newval);} 00070 double samplingFrequency() const {return 1.0/_deltaT;} 00071 double duration() const {return _deltaT * _nSamples;} 00072 inline virtual void setNSamples(int n); 00073 SignalType type() const {return _type;} 00074 void setType(SignalType t) {if(t!=UndefinedSignalType) _type=t;} 00075 inline bool isReal() const; 00076 inline bool isComplex() const; 00077 inline bool isSpectrum() const; 00078 void copyBasicProperties(const DoubleSignal& o); 00079 00080 static SignalType type(QString t); 00081 static QString typeString(SignalType st); 00082 00083 inline void copySamplesFrom(const DoubleSignal * sig); 00084 void copySamplesFrom(const DoubleSignal * sig, double sigStart, double thisStart, double tToCopy); 00085 void copySamplesFrom(const DoubleSignal * sig, int isigStart, int ithisStart, int nToCopy); 00086 void copyAmplitudeFrom(const DoubleSignal * sig); 00087 int add(const DoubleSignal * sig, double sigStart, double thisStart, 00088 double tToCopy, bool checkInvalid=false, double invalidValue=1e99); 00089 int add(const DoubleSignal * sig); 00090 00091 // Funtion to perform signal processing 00092 void fastFourierTransform(SignalType st); 00093 00094 SignalType saveType(SignalType st); 00095 void restoreType(SignalType st); 00096 const DoubleSignal * saveType(SignalType st) const; 00097 void restoreType(const DoubleSignal * sig) const; 00098 00099 bool subtractValue(double val=0.0); 00100 bool filter(const FilterParameters& param); 00101 bool whiten(); 00102 bool stddevClip(double factor); 00103 bool agc(double width); 00104 bool taper(const TimeRange& t, const TaperParameters& param); 00105 bool taper(const TaperParameters& param); 00106 bool abs(); 00107 bool square(); 00108 bool sqrt(); 00109 bool decimateAmplitude(double delta); 00110 bool shift(double dt); 00111 bool unglitch(double threshold); 00112 bool smooth(const SmoothingParameters& param, double minFreq, double maxFreq); 00113 bool setValue(double val); 00114 00115 bool discreteFourierTransform(); 00116 bool overSample(double factor, const DoubleSignal * sig); 00117 bool decimateTime(int factor, const DoubleSignal * sig); 00118 bool hv(const DoubleSignal * zSig, const DoubleSignal * hSig); 00119 00120 ComplexSignal * morletWavelet(const MorletParameters& param) const; 00121 ComplexSignal * positiveComplexSpectrum(bool fullSize=false) const; 00122 ComplexSignal * complexSpectrum() const; 00123 00124 void stalta(double tsta, double tlta, DoubleSignal * ratio) const; 00125 void staltaToKeep(const WindowingParameters::ShortLongTermAverage& param, KeepSignal * keep, double delay) const; 00126 void clipMaxToKeep(double percentage, KeepSignal * keep, double delay) const; 00127 bool correlation(const DoubleSignal * s1, const DoubleSignal * s2, double maxDelay); 00128 bool normalizedCorrelation(const DoubleSignal * s1, const DoubleSignal * s2, double maxDelay); 00129 bool timeCorrelation(const DoubleSignal * s1, const DoubleSignal * s2, double maxDelay); 00130 double correlation(const DoubleSignal * sig, double sigStart, double thisStart, double length) const; 00131 bool convolution(const DoubleSignal * s1, const DoubleSignal * s2); 00132 00133 double amplitudeAt(double abscissa) const; 00134 double interpoleAt(double abscissa) const; 00135 double maximumAmplitude(int itmin, int itmax) const; 00136 int maximumAmplitudeAt(int itmin, int itmax) const; 00137 double averageAmplitude(int itmin, int itmax) const; 00138 double energy(double minF, double maxF) const; 00139 double power(double excludeBorderRatio=0.0) const; 00140 00141 inline Complex complex(const double * samples, int i) const; 00142 inline void setComplex(double * samples, int i, const Complex& value) const; 00143 inline double amplitude2(const double * samples, int i) const; 00144 double amplitude(const double * samples, int i) const; 00145 void setAmplitude(double * samples, int i, double val) const; 00146 inline double phase(const double * samples, int i) const; 00147 inline double re(const double * samples, int i) const; 00148 inline double im(const double * samples, int i) const; 00149 inline void setRe(double * samples, int i, double value) const; 00150 inline void setIm(double * samples, int i, double value) const; 00151 00152 void debugPrint() const; 00153 protected: 00154 mutable QReadWriteLock _propertiesMutex; 00155 double _deltaT; 00156 SignalType _type; 00157 int _nyquistIndex; 00158 private: 00159 bool correlationCheck(const DoubleSignal * s1, const DoubleSignal * s2, double& maxDelay); 00160 bool correlationProcess(const DoubleSignal * s1, const DoubleSignal * s2, 00161 DoubleSignal * s1t, DoubleSignal * s2t); 00162 bool correlationNormalization(const DoubleSignal * s, DoubleSignal * norm) const; 00163 bool smoothFixedTriang(double width, double minF, double maxF); 00164 bool smoothFreqDepTriang(double ratio, double minF, double maxF); 00165 bool smoothKonno(double constant, double minF, double maxF); 00166 static void filb3(double fc, int ak, double pas, double *a, double *b); 00167 static void frecur(int la, double *a, double *xa, int lb, double *b, double *yb, 00168 int ndeb, int nfin, int npas, int isens, double *x, double *y); 00169 static void filrec(double fc, int type, int ndeb, int nfin, int npas, double *x, 00170 double *y, int nfil, double pas); 00171 int commonSamples(double delay, int otherNSamples, int& thisStart, int& otherStart) const; 00172 }; 00173 00174 inline void DoubleSignal::setNSamples(int n) 00175 { 00176 TRACE; 00177 SignalTemplate<double>::setNSamples(n); 00178 _nyquistIndex=_nSamples >> 1; 00179 } 00180 00181 inline double DoubleSignal::re(const double * samples, int i) const 00182 { 00183 TRACE; 00184 if(i<=_nyquistIndex) { 00185 return samples[i]; 00186 } else { 00187 return samples[_nSamples-i]; 00188 } 00189 } 00190 00191 inline double DoubleSignal::im(const double * samples, int i) const 00192 { 00193 TRACE; 00194 if(i==0) { 00195 return 0.0; 00196 } else if(i < _nyquistIndex) { 00197 return samples[ _nSamples -i ]; 00198 } else if(i==_nyquistIndex && !(_nSamples & 0x00000001)) { 00199 return 0.0; // For even sequency, sample at nyquist frequency 00200 } else { 00201 return -samples[ i ]; 00202 } 00203 } 00204 00205 inline void DoubleSignal::setRe(double * samples, int i, double value) const 00206 { 00207 TRACE; 00208 if(i<=_nyquistIndex) { 00209 samples[i]=value; 00210 } else { 00211 samples[_nSamples-i]=value; 00212 } 00213 } 00214 00215 inline void DoubleSignal::setIm(double * samples, int i, double value) const 00216 { 00217 TRACE; 00218 if(i<_nyquistIndex && i>0) { 00219 samples[_nSamples-i]=value; 00220 } 00221 // Not allowed to change the value for all other cases 00222 } 00223 00224 inline Complex DoubleSignal::complex(const double * samples, int i) const 00225 { 00226 TRACE; 00227 if(i==0) { 00228 return samples[i]; 00229 } else if(i<_nyquistIndex) { 00230 return Complex(samples[i], samples[_nSamples-i]); 00231 } else if(i==_nyquistIndex && !(_nSamples & 0x00000001) ) { 00232 return samples[i]; // For even sequency, sample at nyquist frequency 00233 } else { 00234 return Complex(samples[_nSamples-i], -samples[i]); 00235 } 00236 } 00237 00238 inline void DoubleSignal::setComplex(double * samples, int i, const Complex& value) const 00239 { 00240 TRACE; 00241 if(i>0 && i<_nyquistIndex) { 00242 samples[i]=value.re(); 00243 samples[_nSamples-i]=value.im(); 00244 } else { 00245 samples[i]=value.re(); 00246 } 00247 } 00248 00249 inline double DoubleSignal::amplitude2(const double * samples, int i) const 00250 { 00251 TRACE; 00252 if(i==0) { 00253 return samples[i]*samples[i]; 00254 } else if(i!=_nyquistIndex || (_nSamples & 0x00000001)) { 00255 return samples[i]*samples[i] + 00256 samples[_nSamples-i]*samples[_nSamples-i]; 00257 } else { // For even sequency, sample at nyquist frequency 00258 return samples[i]*samples[i]; 00259 } 00260 } 00261 00262 inline double DoubleSignal::amplitude(const double * samples, int i) const 00263 { 00264 return ::sqrt(amplitude2(samples, i)); 00265 } 00266 00267 inline double DoubleSignal::phase(const double * samples, int i) const 00268 { 00269 return complex(samples, i).phase(); 00270 } 00271 00272 inline void DoubleSignal::setAmplitude(double * samples, int i, double val) const 00273 { 00274 if(i==0) { 00275 samples[i]=val; 00276 } else if(i!=_nyquistIndex || (_nSamples & 0x00000001)) { 00277 double ampl=::sqrt(samples[i]*samples[i]+ 00278 samples[_nSamples-i]*samples[_nSamples-i]); 00279 if(ampl==0.0) return; 00280 double factor=val/ampl; 00281 samples[i ]*=factor; 00282 samples[_nSamples-i]*=factor; 00283 } else { // For even sequency, sample at nyquist frequency 00284 samples[i]=val; 00285 } 00286 } 00287 00288 inline bool DoubleSignal::isComplex() const 00289 { 00290 TRACE; 00291 switch (_type) { 00292 case Spectrum: 00293 case Phase: 00294 case CAmpWaveform: 00295 case CPhaseWaveform: 00296 return true; 00297 default: 00298 return false; 00299 } 00300 } 00301 00302 inline bool DoubleSignal::isReal() const 00303 { 00304 TRACE; 00305 switch (_type) { 00306 case Waveform: 00307 case RealSpectrum: 00308 return true; 00309 default: 00310 return false; 00311 } 00312 } 00313 00314 inline bool DoubleSignal::isSpectrum() const 00315 { 00316 TRACE; 00317 switch (_type) { 00318 case Spectrum: 00319 case Phase: 00320 case RealSpectrum: 00321 return true; 00322 default: 00323 return false; 00324 } 00325 } 00326 00327 inline void DoubleSignal::copySamplesFrom(const DoubleSignal * sig) 00328 { 00329 SignalTemplate<double>::copySamplesFrom(sig); 00330 } 00331 00332 } // namespace GeopsyCore 00333 00334 #endif // DOUBLESIGNAL_H
