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: 2011-11-30 00022 ** Authors: 00023 ** Marc Wathelet (ISTerre, Grenoble, France) 00024 ** 00025 ***************************************************************************/ 00026 00027 #ifndef TAPERDELEGATE_H 00028 #define TAPERDELEGATE_H 00029 00030 #include "GeopsyCoreDLLExport.h" 00031 #include "TimeRange.h" 00032 #include "TaperParameters.h" 00033 00034 namespace GeopsyCore { 00035 00036 class TaperParameters; 00037 00038 class GEOPSYCORE_EXPORT TaperDelegate 00039 { 00040 public: 00041 TaperDelegate(const TaperParameters * p); 00042 ~TaperDelegate(); 00043 00044 template<typename sampleType> 00045 void apply(sampleType samples, int signalFirstSample, int signalLastSample, 00046 int windowFirstSample, int windowLastSample) const; 00047 00048 template<typename sampleType> 00049 void apply(sampleType samples, int nSamples, const TimeRange& range, 00050 double samplingFrequency) const; 00051 00052 class Sample 00053 { 00054 public: 00055 Sample(double * samples) {_samples=samples;} 00056 ~Sample() {} 00057 protected: 00058 double * _samples; 00059 }; 00060 00061 class RealSample: public Sample 00062 { 00063 public: 00064 RealSample(double * samples) : Sample(samples) {} 00065 00066 void setNull(int index) {_samples[index]=0;} 00067 void multiply(int index, double factor) {_samples[index]*=factor;} 00068 }; 00069 00070 class ComplexSample: public Sample 00071 { 00072 public: 00073 ComplexSample(int nSamples, double * samples) 00074 : Sample(samples) {_nSamples=nSamples; _nyquistIndex=_nSamples >> 1;} 00075 void setNull(int index) 00076 { 00077 _samples[index]=0.0; 00078 // If index==_nyquistIndex and even number of samples, _nSamples-index==index 00079 // Useless affectation to 0 but it is more effecient than condition testing for all samples. 00080 if(index>0) { 00081 _samples[_nSamples-index]=0; 00082 } 00083 } 00084 void multiply(int index, double factor) 00085 { 00086 _samples[index]*=factor; 00087 // Above shortcut cannot be done without multiplying twice 00088 if(index>0 && (index<_nyquistIndex || (_nSamples & 0x00000001))) { 00089 _samples[_nSamples-index]*=factor; 00090 } 00091 } 00092 protected: 00093 int _nSamples; 00094 int _nyquistIndex; 00095 }; 00096 private: 00097 const TaperParameters * _parameters; 00098 }; 00099 00103 template<typename sampleType> 00104 void TaperDelegate::apply(sampleType samples, int nSamples, const TimeRange& range, 00105 double samplingFrequency) const 00106 { 00107 TRACE; 00108 apply(samples, 0, nSamples-1, (int)round(range.start()*samplingFrequency), 00109 (int)round(range.end()*samplingFrequency)); 00110 } 00111 00122 template<typename sampleType> 00123 void TaperDelegate::apply(sampleType samples, int signalFirstSample, int signalLastSample, 00124 int windowFirstSample, int windowLastSample) const 00125 { 00126 TRACE; 00127 double width=windowLastSample-windowFirstSample; 00128 double widthFactor=1.0/width; 00129 00130 ASSERT(windowFirstSample<=windowLastSample); 00131 if(windowFirstSample>signalLastSample) { 00132 windowFirstSample=signalLastSample; 00133 } 00134 if(windowLastSample>signalLastSample) { 00135 windowLastSample=signalLastSample; 00136 } 00137 00138 int i=signalFirstSample; 00139 00140 // Sets signal eventually to null before the window 00141 if(!_parameters->reversed()) { 00142 for(; i<=windowFirstSample; i++) { 00143 samples.setNull(i); 00144 } 00145 } 00146 00147 // Transforms signal according to window type within window range 00148 switch(_parameters->window()) { 00149 case TaperParameters::Rectangular: 00150 if(_parameters->reversed()) { 00151 for(; i<=windowLastSample; i++) { 00152 samples.setNull(i); 00153 } 00154 } else { 00155 i=windowLastSample; 00156 } 00157 break; 00158 case TaperParameters::Hann: 00159 widthFactor*=2.0*M_PI; 00160 if(_parameters->reversed()) { 00161 for(; i<=windowLastSample; i++) { 00162 samples.multiply(i, 0.5+0.5*cos((i-windowFirstSample)*widthFactor)); 00163 } 00164 } else { 00165 for(; i<=windowLastSample; i++) { 00166 samples.multiply(i, 0.5-0.5*cos((i-windowFirstSample)*widthFactor)); 00167 } 00168 } 00169 break; 00170 case TaperParameters::Hamming: 00171 widthFactor*=2.0*M_PI; 00172 if(_parameters->reversed()) { 00173 for(; i<=windowLastSample; i++) { 00174 samples.multiply(i, 0.46+0.46*cos((i-windowFirstSample)*widthFactor)); 00175 } 00176 } else { 00177 for(; i<=windowLastSample; i++) { 00178 samples.multiply(i, 0.54-0.46*cos((i-windowFirstSample)*widthFactor)); 00179 } 00180 } 00181 break; 00182 case TaperParameters::Tukey: { 00183 ASSERT(_parameters->alpha()<=1.0); 00184 ASSERT(_parameters->alpha()>0.0); 00185 double term=2.0/_parameters->alpha(); 00186 double factor=term*widthFactor*M_PI; 00187 term=M_PI*(1.0-term); 00188 00189 int windowLastSampleLeft=(int)floor(_parameters->alpha()*width*0.5); 00190 if(windowLastSampleLeft>signalLastSample) { 00191 windowLastSampleLeft=signalLastSample; 00192 } 00193 int windowFirstSampleRight=(int)ceil(width*(1.0-_parameters->alpha()*0.5)); 00194 if(windowFirstSampleRight>signalLastSample) { 00195 windowFirstSampleRight=signalLastSample; 00196 } 00197 00198 if(_parameters->reversed()) { 00199 for(; i<=windowLastSampleLeft; i++) { 00200 samples.multiply(i, 0.5+0.5*cos(factor*(i-windowFirstSample))); 00201 } 00202 for(; i<=windowFirstSampleRight; i++) { 00203 samples.setNull(i); 00204 } 00205 for(; i<=windowLastSample; i++) { 00206 samples.multiply(i, 0.5-0.5*cos(factor*(i-windowFirstSample)+term)); 00207 } 00208 } else { 00209 for(; i<=windowLastSampleLeft; i++) { 00210 samples.multiply(i, 0.5-0.5*cos(factor*(i-windowFirstSample))); 00211 } 00212 for(i=windowFirstSampleRight; i<=windowLastSample; i++) { 00213 samples.multiply(i, 0.5+0.5*cos(factor*(i-windowFirstSample)+term)); 00214 } 00215 } 00216 } 00217 break; 00218 case TaperParameters::Cosine: 00219 widthFactor*=M_PI; 00220 if(_parameters->reversed()) { 00221 for(; i<=windowLastSample; i++) { 00222 samples.multiply(i, 1.0-sin((i-windowFirstSample)*widthFactor)); 00223 } 00224 } else { 00225 for(; i<=windowLastSample; i++) { 00226 samples.multiply(i, sin((i-windowFirstSample)*widthFactor)); 00227 } 00228 } 00229 break; 00230 case TaperParameters::Lanczos: { 00231 double factor=M_PI*2.0*widthFactor; 00232 if(_parameters->reversed()) { 00233 for(; i<=windowLastSample; i++) { 00234 double x=factor*(i-windowFirstSample)-M_PI; 00235 samples.multiply(i, 1.0-sin(x)/x); 00236 } 00237 } else { 00238 for(; i<=windowLastSample; i++) { 00239 double x=factor*(i-windowFirstSample)-M_PI; 00240 samples.multiply(i, sin(x)/x); 00241 } 00242 } 00243 } 00244 break; 00245 case TaperParameters::Triangular: { 00246 double factor=2.0*widthFactor; 00247 double term=0.5*width; 00248 if(_parameters->reversed()) { 00249 for(; i<=windowLastSample; i++) { 00250 samples.multiply(i, factor*fabs(i-windowFirstSample-term)); 00251 } 00252 } else { 00253 for(; i<=windowLastSample; i++) { 00254 samples.multiply(i, 1.0-factor*fabs(i-windowFirstSample-term)); 00255 } 00256 } 00257 } 00258 break; 00259 case TaperParameters::Bartlett: { 00260 double factor=2.0/(width+1.0); 00261 double term=0.5*width; 00262 if(_parameters->reversed()) { 00263 for(; i<=windowLastSample; i++) { 00264 samples.multiply(i, factor*fabs(i-windowFirstSample-term)); 00265 } 00266 } else { 00267 for(; i<=windowLastSample; i++) { 00268 samples.multiply(i, 1.0-factor*fabs(i-windowFirstSample-term)); 00269 } 00270 } 00271 } 00272 break; 00273 case TaperParameters::Gaussian: { 00274 ASSERT(_parameters->sigma()<=0.5); 00275 double factor=2.0*widthFactor; 00276 double invSigma=1.0/_parameters->sigma(); 00277 if(_parameters->reversed()) { 00278 for(; i<=windowLastSample; i++) { 00279 double x=(factor*(i-windowFirstSample)-1.0)*invSigma; 00280 samples.multiply(i, 1.0-exp(-0.5*x*x)); 00281 } 00282 } else { 00283 for(; i<=windowLastSample; i++) { 00284 double x=(factor*(i-windowFirstSample)-1.0)*invSigma; 00285 samples.multiply(i, exp(-0.5*x*x)); 00286 } 00287 } 00288 } 00289 break; 00290 case TaperParameters::BartlettHann: 00291 if(_parameters->reversed()) { 00292 for(; i<=windowLastSample; i++) { 00293 double x=(i-windowFirstSample)*widthFactor; 00294 samples.multiply(i, 0.38+0.48*fabs(x-0.5)+0.38*cos(2.0*M_PI*x)); 00295 } 00296 } else { 00297 for(; i<=windowLastSample; i++) { 00298 double x=(i-windowFirstSample)*widthFactor; 00299 samples.multiply(i, 0.62-0.48*fabs(x-0.5)-0.38*cos(2.0*M_PI*x)); 00300 } 00301 } 00302 break; 00303 case TaperParameters::Blackman: { 00304 double a0=(1.0-_parameters->alpha())*0.5; 00305 double a1=0.5; 00306 double a2=_parameters->alpha()*0.5; 00307 double f1=2.0*M_PI*widthFactor; 00308 double f2=f1*2.0; 00309 if(_parameters->reversed()) { 00310 for(; i<=windowLastSample; i++) { 00311 double x=i-windowFirstSample; 00312 samples.multiply(i, 1.0-(a0-a1*cos(x*f1)+a2*cos(x*f2))); 00313 } 00314 } else { 00315 for(; i<=windowLastSample; i++) { 00316 double x=i-windowFirstSample; 00317 samples.multiply(i, a0-a1*cos(x*f1)+a2*cos(x*f2)); 00318 } 00319 } 00320 } 00321 break; 00322 case TaperParameters::Nuttall: { 00323 double f1=2.0*M_PI*widthFactor; 00324 double f2=f1*2.0; 00325 double f3=f1*3.0; 00326 if(_parameters->reversed()) { 00327 for(; i<=windowLastSample; i++) { 00328 double x=i-windowFirstSample; 00329 samples.multiply(i, 1.0-(0.355768-0.487396*cos(x*f1)+0.144232*cos(x*f2)-0.012604*cos(x*f3))); 00330 } 00331 } else { 00332 for(; i<=windowLastSample; i++) { 00333 double x=i-windowFirstSample; 00334 samples.multiply(i, 0.355768-0.487396*cos(x*f1)+0.144232*cos(x*f2)-0.012604*cos(x*f3)); 00335 } 00336 } 00337 } 00338 break; 00339 case TaperParameters::BlackmanHarris: { 00340 double f1=2.0*M_PI*widthFactor; 00341 double f2=f1*2.0; 00342 double f3=f1*3.0; 00343 if(_parameters->reversed()) { 00344 for(; i<=windowLastSample; i++) { 00345 double x=i-windowFirstSample; 00346 samples.multiply(i, 1.0-(0.35875-0.48829*cos(x*f1)+0.14128*cos(x*f2)-0.01168*cos(x*f3))); 00347 } 00348 } else { 00349 for(; i<=windowLastSample; i++) { 00350 double x=i-windowFirstSample; 00351 samples.multiply(i, 0.35875-0.48829*cos(x*f1)+0.14128*cos(x*f2)-0.01168*cos(x*f3)); 00352 } 00353 } 00354 } 00355 break; 00356 case TaperParameters::BlackmanNuttall: { 00357 double f1=2.0*M_PI*widthFactor; 00358 double f2=f1*2.0; 00359 double f3=f1*3.0; 00360 if(_parameters->reversed()) { 00361 for(; i<=windowLastSample; i++) { 00362 double x=i-windowFirstSample; 00363 samples.multiply(i, 1.0-(0.3635819-0.4891775*cos(x*f1)+0.1365995*cos(x*f2)-0.0106411*cos(x*f3))); 00364 } 00365 } else { 00366 for(; i<=windowLastSample; i++) { 00367 double x=i-windowFirstSample; 00368 samples.multiply(i, 0.3635819-0.4891775*cos(x*f1)+0.1365995*cos(x*f2)-0.0106411*cos(x*f3)); 00369 } 00370 } 00371 } 00372 break; 00373 case TaperParameters::FlatTop: { 00374 double f1=2.0*M_PI*widthFactor; 00375 double f2=f1*2.0; 00376 double f3=f1*3.0; 00377 double f4=f1*4.0; 00378 if(_parameters->reversed()) { 00379 for(; i<=windowLastSample; i++) { 00380 double x=i-windowFirstSample; 00381 samples.multiply(i, 1.0-(1.0-1.93*cos(x*f1)+1.29*cos(x*f2)-0.388*cos(x*f3)+0.032*cos(x*f4))); 00382 } 00383 } else { 00384 for(; i<=windowLastSample; i++) { 00385 double x=i-windowFirstSample; 00386 samples.multiply(i, 1.0-1.93*cos(x*f1)+1.29*cos(x*f2)-0.388*cos(x*f3)+0.032*cos(x*f4)); 00387 } 00388 } 00389 } 00390 break; 00391 } 00392 00393 // Sets signal eventually to null after the window 00394 if(!_parameters->reversed()) { 00395 for(; i<=signalLastSample; i++) { 00396 samples.setNull(i); 00397 } 00398 } 00399 } 00400 00401 } // namespace GeopsyCore 00402 00403 #endif // TAPERDELEGATE_H
