1 #ifndef CLITKSIGNAL_CXX
2 #define CLITKSIGNAL_CXX
4 #include "clitkSignal.h"
9 //---------------------------------------------------------------------
10 Signal::Signal(const Signal & s) { // Copy
13 //---------------------------------------------------------------------
16 //---------------------------------------------------------------------
17 void Signal::CopyFrom(const Signal & s) {
18 SetSamplingPeriod(s.GetSamplingPeriod());
19 // DD(GetSamplingPeriod());
21 for(uint i=0; i<size(); i++) {
25 //---------------------------------------------------------------------
29 //---------------------------------------------------------------------
30 void Signal::Read(string fileName){
31 std::ifstream signalStream(fileName.c_str());
32 SignalValueType point;
33 if(!signalStream.is_open()){
34 std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
37 skipComment(signalStream);
38 while(!signalStream.eof()) {
39 skipComment(signalStream);
40 signalStream >> point;
41 skipComment(signalStream);
42 m_Data.push_back(point);
43 // cout << point << endl;
46 SetSamplingPeriod(1.);
48 //---------------------------------------------------------------------
51 //---------------------------------------------------------------------
52 bool Signal::Read(string fileName, int col){
53 ifstream signalStream(fileName.c_str());
54 SignalValueType point;
55 if(!signalStream.is_open()){
56 std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
59 skipComment(signalStream);
60 while(!signalStream.eof()) {
61 skipComment(signalStream);
65 std::getline(signalStream, line);
68 istringstream iss(line);
69 for(int i=0; i<col; i++) { // skip col-1 first column
76 std::cerr << "ERROR: no col n" << col << " in the line '" << line << "' ?" << std::endl;
77 return false; //exit(0);
81 skipComment(signalStream);
82 m_Data.push_back(point);
83 // cout << point << endl;
86 SetSamplingPeriod(1.);
89 //---------------------------------------------------------------------
93 //---------------------------------------------------------------------
94 //Convert 1D image to signal
95 Signal Signal::ConvertImageToSignal( Signal::ImageType::Pointer image)
100 //make an image iterator
101 itk::ImageRegionConstIterator<ImageType> it(image,image->GetLargestPossibleRegion());
107 signal.push_back(it.Get());
112 signal.SetSamplingPeriod(image->GetSpacing()[0]);
116 //---------------------------------------------------------------------
119 //---------------------------------------------------------------------
120 //Convert 1D signal to image
121 Signal::ImageType::Pointer Signal::ConvertSignalToImage(Signal signal)
124 ImageType::Pointer image =ImageType::New();
125 ImageType::RegionType region;
126 ImageType::RegionType::IndexType index;
128 ImageType::RegionType::SizeType size;
129 size[0]=signal.size();
131 region.SetIndex(index);
132 region.SetSize(size);
134 image->SetRegions(region);
137 //make an image iterator
138 itk::ImageRegionIterator<ImageType> mIt(image,image->GetLargestPossibleRegion());
141 //make a signal iterator
142 Signal::const_iterator sIt=signal.begin();
145 while(sIt!=signal.end())
153 ImageType::SpacingType spacing;
154 spacing[0]=signal.GetSamplingPeriod();
155 image->SetSpacing(spacing);
160 //---------------------------------------------------------------------
164 //---------------------------------------------------------------------
165 void Signal::Write(const string fileName){
166 ofstream signalStream(fileName.c_str());
167 if(!signalStream.is_open()){
168 cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Write" << endl;
174 signalStream << *it << endl;
177 signalStream.close();
179 //---------------------------------------------------------------------
182 //---------------------------------------------------------------------
183 Signal & Signal::operator/=(Signal & div){
184 if(size()!= div.size())
186 std::cerr << "Error: signal size must be the same!" << std::endl;
197 cerr << "Division by 0 in operator/= skipped" << endl;
202 //---------------------------------------------------------------------
205 //---------------------------------------------------------------------
206 Signal & Signal::operator*=(Signal & mul){
207 if(size()!= mul.size())
209 std::cerr << "Error: signal size must be the same!" << std::endl;
221 //---------------------------------------------------------------------
224 //---------------------------------------------------------------------
225 Signal Signal::Normalize(double newMin,double newMax){
226 Signal temp (m_SamplingPeriod);
227 vector<double> extrema=GetGlobalMinMax();
228 iterator itSig=begin();
230 temp.push_back( ((*itSig)-extrema[0])*(newMax-newMin)/(extrema[1]-extrema[0]) + newMin );
235 //---------------------------------------------------------------------
238 //---------------------------------------------------------------------
239 vector<double> Signal::GetGlobalMinMax() const {
240 vector<double> extrema(2);
242 cerr << "ERROR: GetExtrema / No signal" << endl;
245 extrema[0]=m_Data[0];
246 extrema[1]=m_Data[0];
247 for(unsigned int i=1;i<m_Data.size();i++){
248 if(extrema[0]>m_Data[i]) extrema[0]=m_Data[i];
249 if(extrema[1]<m_Data[i]) extrema[1]=m_Data[i];
253 //---------------------------------------------------------------------
256 //---------------------------------------------------------------------
257 double Signal::GetGlobalMean() const {
259 for(unsigned int i=0;i<m_Data.size();i++){
262 return m/(double)m_Data.size();
264 //---------------------------------------------------------------------
267 //---------------------------------------------------------------------
268 Signal Signal::MovingAverageFilter ( unsigned int length) {
270 Signal temp(m_SamplingPeriod);
272 for (unsigned int i=0; i <size(); i++)
274 double accumulator=0.;
275 unsigned int scale=0;
276 for (unsigned int j= max(0,static_cast<int>(i)-static_cast<int>(length)); j<=min(size(), i+length); j++)
278 accumulator+=m_Data[j];
281 temp.push_back(accumulator/scale);
285 //---------------------------------------------------------------------
288 //---------------------------------------------------------------------
289 Signal Signal::GaussLikeFilter ( ) {
291 Signal temp(m_SamplingPeriod);
296 //first sample: mirrorring BC
297 temp.push_back((2.*m_Data[0]+2*m_Data[1])/4.);
300 for (unsigned int i=1; i <size()-1; i++)
301 temp.push_back((2.*m_Data[i]+m_Data[i-1]+m_Data[i+1])/4.);
303 //end sample: mirrorring BC
304 temp.push_back((2.*m_Data[size()-1]+2*m_Data[size()-2])/4.);
308 //---------------------------------------------------------------------
311 //---------------------------------------------------------------------
312 Signal Signal::NormalizeMeanStdDev(double newMean,double newStdDev)
314 iterator itSig=begin();
315 double sum=0, sum2=0;
316 Signal temp(m_SamplingPeriod);
318 //Calculate the mean and the std dev
322 sum2 += (*itSig) * (*itSig);
325 double oldMean=sum/size();
326 double oldStdDev=sqrt(sum2/size()-oldMean*oldMean);
329 double a = newStdDev/oldStdDev;
330 double b = newMean - a * oldMean;
334 temp.push_back( a *(*itSig) + b );
339 //---------------------------------------------------------------------
342 //---------------------------------------------------------------------
343 // Signal Signal::HighPassFilter (double sampPeriod, double cutOffFrequency )
346 // Signal temp(m_SamplingPeriod);
347 // temp.resize(size());
350 // SIGNAL_FFT_TYPE fft;
352 // //calculate the cut off frequency
353 // unsigned int samp=lrint(cutOffFrequency*static_cast<double>(size())*sampPeriod);
355 // //forward fft with empty fft
356 // if(fft.size()==0) OneDForwardFourier(*this, fft);
358 // //remove the frequencies
359 // for(unsigned int i=0;i<samp && i<fft.size();i++)
360 // fft[i]=complex<double>(0.,0.);
362 // //backward with remaining frequencies
363 // OneDBackwardFourier(fft,temp);
366 //---------------------------------------------------------------------
369 //---------------------------------------------------------------------
370 // Signal Signal::LowPassFilter (double sampPeriod, double cutOffFrequency )
373 // Signal temp(m_SamplingPeriod);
374 // temp.resize(size());
377 // SIGNAL_FFT_TYPE fft;
379 // //calculate the cut off frequency
380 // unsigned int samp=lrint(cutOffFrequency*static_cast<double>(size())*sampPeriod);
382 // //forward fft with empty fft
383 // if(fft.size()==0) OneDForwardFourier(*this, fft);
384 // unsigned int fsize=fft.size();
386 // //remove the frequencies
387 // unsigned int limit=min (samp, fsize);
388 // for(unsigned int i=limit;i<fft.size();i++)
389 // fft[i]=complex<double>(0.,0.);
391 // //backward with remaining frequencies
392 // OneDBackwardFourier(fft,temp);
395 //---------------------------------------------------------------------
398 //---------------------------------------------------------------------
399 // void Signal::OneDForwardFourier(const Signal& input, SIGNAL_FFT_TYPE & fft)
401 // //Create output array
402 // fft.resize(input.size()/2+1);
404 // double *tempCopy=new double[size()];
405 // copy(begin(), end(), tempCopy);
407 // //Forward Fourier Transform
409 // p=fftw_plan_dft_r2c_1d(size(),tempCopy,reinterpret_cast<fftw_complex*>(&(fft[0])),FFTW_ESTIMATE);
411 // fftw_destroy_plan(p);
412 // //delete tempCopy;
415 //---------------------------------------------------------------------
418 //---------------------------------------------------------------------
419 // void Signal::OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output)
424 // p=fftw_plan_dft_c2r_1d(output.size(),reinterpret_cast<fftw_complex*>(&(fft[0])),&(output[0]),FFTW_ESTIMATE);
426 // fftw_destroy_plan(p);
428 // vector<double>::iterator it=output.begin();
429 // while(it!=output.end()){
430 // *it /= (double)output.size();
435 //---------------------------------------------------------------------
438 //---------------------------------------------------------------------
439 // double Signal::MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft)
442 // if(fft.size()==0) OneDForwardFourier(sig,fft);
444 // double amplitude, amplitudeMax=abs(fft[1]);
445 // for(unsigned int i=1;i<fft.size();i++){
446 // amplitude=abs(fft[i]);
447 // if(amplitude>amplitudeMax){
449 // amplitudeMax=amplitude;
452 // return ((double)(posMax)/((double)sig.size()*sig.GetSamplingPeriod()));
454 //---------------------------------------------------------------------
457 //---------------------------------------------------------------------
458 Signal Signal::DetectLocalExtrema(unsigned int width)
460 Signal temp(m_SamplingPeriod);
462 unsigned int upper, lower;
464 //has to be at least 1
465 width=max(static_cast<int>(width), 1);
467 for(unsigned int i=0 ; i < size() ; i++){
471 for(unsigned int j=1; j< width+1; j++)
474 upper = min( size(), i+j);
475 lower = max( static_cast<int>(0), (int)i-(int)j);
478 if( ! (m_Data[i] >= m_Data[lower] && m_Data[i] >= m_Data[upper]))isMax=false;
481 if( ! (m_Data[i]<= m_Data[lower] && m_Data[i] <= m_Data[upper])) isMin=false;
483 //if neither, go to the next value
484 if( (!isMax) && (!isMin)) break;
487 if (isMax) temp.push_back(1);
488 else if (isMin) temp.push_back(0);
489 else temp.push_back(0.5);
493 //---------------------------------------------------------------------
496 //---------------------------------------------------------------------
497 Signal Signal::LimPhase()
500 //phase is defined as going from 0 to 1 linearly between end expiration
501 Signal phase(m_SamplingPeriod);
502 phase.resize(size());
504 unsigned int firstBeginCycle=0;
505 unsigned int firstEndCycle=0;
506 unsigned int beginCycle=0;
508 //=========================================
509 //Search the first value in extrema not 0.5
510 while(m_Data[beginCycle]!=1)
515 //We search the corresponding end
516 unsigned int endCycle=beginCycle+1;
517 while(endCycle != size() && m_Data[endCycle]!=1){
522 //============================================================
523 //Calculate phase at the beginning (before the first extremum)
524 for(unsigned int i=0 ; i<beginCycle ; i++)
527 //if before first cycle is longer than first cycle: let it go from 0 to 1
528 if((double)beginCycle/(double)(endCycle-beginCycle) > 1)
530 phase[i] = (double)(i-0)/(double)(beginCycle-0);
532 //copy the phase values later
535 firstBeginCycle=beginCycle;
536 firstEndCycle=endCycle;
541 //===================================================================
543 while(endCycle != size()){
545 //fill between extrema
546 for(unsigned int i=beginCycle ; i<endCycle ; i++){
547 phase[i] = (double)(i-beginCycle)/(double)(endCycle-beginCycle);
551 beginCycle = endCycle++;
552 while(endCycle != size() && m_Data[endCycle]!=1)
556 //===================================================================
557 //Calculate phase at the end (after the last extremum)
558 endCycle = beginCycle--;
560 //count backwards till the previous
561 while(m_Data[beginCycle]!=1) beginCycle--;
562 for(unsigned int i=endCycle ; i<size() ; i++)
565 //after last extrema is longer then last cycle
566 if((double)(size()-endCycle)/(double)(endCycle-beginCycle) > 1){
568 //make the last part go till 1
569 phase[i] = (double)(i-endCycle)/(double)(size()-endCycle);
572 //the last part is shorter, copy the last cycle values
574 phase[i] = phase[i -(endCycle-beginCycle)];
578 //===================================================================
579 //check it some remains to be copied in the beginning
580 if (firstBeginCycle!=0)
582 for(unsigned int i=0 ; i<firstBeginCycle ; i++)
583 phase[i] =phase[i+firstEndCycle-firstBeginCycle];
588 //---------------------------------------------------------------------
592 //---------------------------------------------------------------------
593 Signal Signal::MonPhase()
596 //phase is defined as going from 0 to 1 linearly between end expiration
597 Signal phase(m_SamplingPeriod);
598 phase.resize(size());
600 unsigned int firstBeginCycle=0;
601 unsigned int firstEndCycle=0;
602 unsigned int cycleCounter=0;
603 unsigned int beginCycle=0;
605 //===================================================================
606 //Search the first value in extrema not 0.5
607 while(m_Data[beginCycle]!=1)
612 //We search the corresponding end
613 unsigned int endCycle=beginCycle+1;
614 while(endCycle != size() && m_Data[endCycle]!=1){
619 //===================================================================
620 //Calculate phase at the beginning (before the first extremum)
621 for(unsigned int i=0 ; i<beginCycle ; i++)
624 //if before first cycle is longer than first cycle: let it go from 0 to 1
625 if((double)beginCycle/(double)(endCycle-beginCycle) > 1)
627 phase[i] = (double)(i-0)/(double)(beginCycle-0);
629 //copy the phase values later
632 firstBeginCycle=beginCycle;
633 firstEndCycle=endCycle;
638 //===================================================================
640 while(endCycle != size()){
643 //fill between extrema
644 for(unsigned int i=beginCycle ; i<endCycle ; i++)
645 phase[i] = (double)cycleCounter+(double)(i-beginCycle)/(double)(endCycle-beginCycle);
649 beginCycle = endCycle++;
650 while(endCycle != size() && m_Data[endCycle]!=1)
654 //===================================================================
655 //Calculate phase at the end (after the last extremum)
656 endCycle = beginCycle--;
658 //count backwards till the previous
659 while(m_Data[beginCycle]!=1) beginCycle--;
660 for(unsigned int i=endCycle ; i<size() ; i++)
663 //after last extrema is longer then last cycle
664 if((double)(size()-endCycle)/(double)(endCycle-beginCycle) > 1){
666 //make the last part go till 1
667 phase[i] = (double)cycleCounter+(double)(i-endCycle)/(double)(size()-endCycle);
670 //the last part is shorter, copy the last cycle values
672 phase[i] = phase[i -(endCycle-beginCycle)]+1;
676 //===================================================================
677 //check it some remains to be copied in the beginning
678 if (firstBeginCycle!=0)
680 for(unsigned int i=0 ; i<firstBeginCycle ; i++)
681 phase[i] =phase[i+firstEndCycle-firstBeginCycle]-1;
686 //---------------------------------------------------------------------
689 //---------------------------------------------------------------------
690 Signal Signal::MonPhaseDE(double eEPhaseValue, double eIPhaseValue)
692 //First calculate monPhase
693 Signal monPhase=this->MonPhase();
695 //Create an empty signal
696 Signal phase(size(), -1);
697 phase.SetSamplingPeriod(m_SamplingPeriod);
699 //Fill in the values at the extrema position
700 iterator phaseIt=phase.begin();
701 iterator monPhaseIt=monPhase.begin();
702 iterator extremaIt =begin();
703 while (extremaIt!= end())
705 if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
706 else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
707 extremaIt++; phaseIt++;monPhaseIt++;
713 //---------------------------------------------------------------------
716 //---------------------------------------------------------------------
717 Signal Signal::LinearlyInterpolateScatteredValues()
719 //Linearly interpolate the values in between
721 unsigned int beginCycle=0;
722 unsigned int endCycle=0;
724 //Create a new signal
725 Signal temp(size(),-1);
726 temp.SetSamplingPeriod(m_SamplingPeriod);
728 //start from the first value
729 while (m_Data[i]==-1)i++;
734 while ( (m_Data[i]==-1) && (i<size()) )i++;
740 for (unsigned int k=beginCycle;k<endCycle+1; k++)
741 temp[k]=m_Data[beginCycle]+(double)(k-beginCycle)*
742 (m_Data[endCycle]-m_Data[beginCycle])/(double)(endCycle-beginCycle);
747 while( (i< size()) && (m_Data[i]==-1) ) i++;
751 if (beginCycle!= size()-1)
753 //continue with same slope
754 for (unsigned int k=beginCycle+1; k<size();k++)
755 temp[k]=temp[beginCycle]+(double)(k-beginCycle)*(temp[beginCycle]-temp[beginCycle-1]);
763 while (temp[i]==-1)i++;
766 while (m_Data[i]==-1)i++;
769 //if the first filled half cycle is longer, copy it
770 if(beginCycle<(endCycle-beginCycle))
772 for (unsigned int k=0; k< beginCycle;k++)
773 temp[k]=temp[k+endCycle-beginCycle];
776 //if the first filled half cycle is longer, start from zero
779 for (unsigned int k=0; k< beginCycle;k++)
780 temp[k]=k*temp[beginCycle]/(beginCycle);
786 //---------------------------------------------------------------------
790 // //---------------------------------------------------------------------
791 // Signal Signal::ApproximateScatteredValuesWithBSplines(unsigned int splineOrder, unsigned int numberOfControlPoints)
793 // //filter requires a vector pixelType
794 // typedef itk::PointSet<VectorType, 1> PointSetType;
795 // PointSetType::Pointer pointSet=PointSetType::New();
796 // typedef PointSetType::PointType PointType;
799 // unsigned int pointIndex=0;
805 // p[0]= i;//JV spacing is allways 1
806 // pointSet->SetPoint( pointIndex, p );
807 // pointSet->SetPointData( pointIndex, m_Data[i] );
813 // //---------------------------------------------------------------------
814 // Signal Signal::ApproximateScatteredValuesWithBSplines(unsigned int splineOrder, unsigned int numberOfControlPoints)
816 // //filter requires a vector pixelType
817 // typedef itk::PointSet<VectorType, 1> PointSetType;
818 // PointSetType::Pointer pointSet=PointSetType::New();
819 // typedef PointSetType::PointType PointType;
822 // unsigned int pointIndex=0;
828 // p[0]= i;//JV spacing is allways 1
829 // pointSet->SetPoint( pointIndex, p );
830 // pointSet->SetPointData( pointIndex, m_Data[i] );
836 // //define the output signal properties
837 // ImageType::RegionType::SizeType outputSize;
838 // outputSize[0]= size();
839 // ImageType::PointType outputOrigin;
840 // outputOrigin[0]=0.0;//JV may need to be changed
841 // ImageType::SpacingType outputSpacing;
842 // outputSpacing[0]=1; //JV add relation to the original signal spacing
845 // typedef itk::BSplineScatteredDataPointSetToImageFilter< PointSetType, VectorImageType > PointSetToImageFilterType;
846 // PointSetToImageFilterType::Pointer pointSetToImageFilter= PointSetToImageFilterType::New();
847 // pointSetToImageFilter->SetInput(pointSet);
848 // pointSetToImageFilter->SetSplineOrder(splineOrder);//JV
849 // pointSetToImageFilter->SetSize(outputSize);
850 // pointSetToImageFilter->SetOrigin(outputOrigin);
851 // pointSetToImageFilter->SetSpacing(outputSpacing);
854 // itk::FixedArray<unsigned int,1> num;
855 // num[0]=numberOfControlPoints;
856 // pointSetToImageFilter->SetNumberOfControlPoints(num);//JV
857 // pointSetToImageFilter->Update();
858 // VectorImageType::Pointer approximatedSignal=pointSetToImageFilter->GetOutput();
860 // //Convert and return
861 // return ConvertVectorImageToSignal(approximatedSignal);
863 // //---------------------------------------------------------------------
866 //---------------------------------------------------------------------
867 Signal Signal::ConvertVectorImageToSignal(VectorImageType::Pointer image)
872 //make an image iterator
873 itk::ImageRegionConstIterator<VectorImageType> it(image,image->GetLargestPossibleRegion());
879 signal.push_back(it.Get()[0]);
884 signal.SetSamplingPeriod(image->GetSpacing()[0]);
888 //---------------------------------------------------------------------
891 //---------------------------------------------------------------------
892 Signal Signal::LimitSignalRange()
895 Signal signal(m_SamplingPeriod);
899 signal.push_back(*it-floor(*it));
906 //---------------------------------------------------------------------
907 double Signal::SSD(const Signal &sig2) const{
908 if(sig2.size() != size()){
909 cerr << "ERROR in Signal::SSD: signals don't have the same size" << endl;
913 for(unsigned int i=0;i<size();i++){
914 result+=pow(sig2[i]-m_Data[i],2);
920 //---------------------------------------------------------------------
923 //---------------------------------------------------------------------
924 void Signal::AddValue(double v) {
925 for(unsigned int i=0;i<size();i++) {
929 //---------------------------------------------------------------------
932 //---------------------------------------------------------------------
933 void Signal::ComputeAugmentedSpace(clitk::Signal & outputX,
934 clitk::Signal & outputY,
935 unsigned int delay) const {
936 if (size() <= delay) {
937 std::cerr << "Error in signal length is " << size()
938 << " while delay is " << delay << " : too short. Abort." << std::endl;
941 outputX.resize(size()-delay);
942 outputY.resize(size()-delay);
943 for(unsigned int i=0; i<outputX.size(); i++) {
944 outputX[i] = m_Data[i+delay];
945 outputY[i] = m_Data[i];
948 //---------------------------------------------------------------------
951 //---------------------------------------------------------------------
952 void Signal::ResampleWithTimeWarpTo(clitk::Signal & output, bool lin) const {
953 double sp = output.GetSamplingPeriod()/output.GetTotalTimeDuration()*GetTotalTimeDuration();
954 // DD(GetTotalTimeDuration());
955 // DD(output.GetTotalTimeDuration());
958 for(uint i=0; i<output.size(); i++) {
959 output[i] = GetValueAtLin(i*sp);
963 for(uint i=0; i<output.size(); i++) {
964 output[i] = GetValueAtNN(i*sp);
968 //---------------------------------------------------------------------
971 //---------------------------------------------------------------------
972 double Signal::GetValueAtNN(double t) const {
973 int index = (int)lrint(t/GetSamplingPeriod()); // closest point
975 if ((index<0) || (index>=(int)size())) {
976 FATAL("out of bound for time " << t << "in this signal having "
977 << size() << " values with sampling " << GetSamplingPeriod() << std::endl);
979 return (*this)[index];
981 //---------------------------------------------------------------------
984 //---------------------------------------------------------------------
985 double Signal::GetValueAtLin(double t) const {
986 double ti = t/GetSamplingPeriod();
989 int index = (int)floor(ti); // floor point
992 FATAL("out of bound for time " << t << " in this signal having "
993 << size() << " values with sampling " << GetSamplingPeriod() << std::endl);
995 if (index>=(int)size()-1) { // last value
996 return (*this)[index];
1000 // DD((*this)[index]);
1001 // DD((*this)[index+1]);
1002 return ((ti-index)*((*this)[index+1]) + (1.0-(ti-index))*((*this)[index]));
1005 //---------------------------------------------------------------------
1008 //---------------------------------------------------------------------
1009 void clitk::Signal::Shift(int d, int length) {
1010 std::cout << "Shift d=" << d << " length = " << length << std::endl;
1011 // int d = (int)lrint(s/GetSamplingPeriod()); // closest integer delta
1014 temp.resize(length);//(uint)size());
1015 temp.SetSamplingPeriod(GetSamplingPeriod());
1016 for(uint i=0; i<temp.size(); i++) {
1017 int j = (i+d) % size();
1020 temp[i] = (*this)[j];
1024 //---------------------------------------------------------------------
1025 //---------------------------------------------------------------------
1026 double Signal::UnormalizedXcov(const Signal &input2, Signal &output) const {
1032 assert(size() == input2.size()); //les signaux à comparer doivent avoir la même taille
1035 //calcul de la moyenne
1036 double mean1=0.,mean2=0.;
1037 mean1=this->GetGlobalMean();
1038 mean2=input2.GetGlobalMean();
1040 for (unsigned int k= 0;k<N;k++){
1042 for (unsigned int n=0;n<=k;n++){
1043 value=value+((*this)[n] - mean1)*(input2[n+N-k-1] - mean2);
1047 for (unsigned int k=N;k<max-1;k++){
1049 for (unsigned int n=0;n<max-k-1;n++){
1050 value=value+((*this)[n+k-N+1] - mean1)*(input2[n] - mean2); //*input2??
1060 //---------------------------------------------------------------------
1061 //---------------------------------------------------------------------
1062 double Signal::UnormalizedCorrelation(const Signal &input2) const {
1067 assert(size() == input2.size()); //les signaux à comparer doivent avoir la même taille
1069 //calcul de la moyenne
1070 double mean1=0.,mean2=0.;
1071 mean1=this->GetGlobalMean();
1072 mean2=input2.GetGlobalMean();
1074 for (unsigned int n=0;n<N;n++){
1075 value=value+((*this)[n] - mean1)*(input2[n] - mean2);
1080 //---------------------------------------------------------------------
1081 //---------------------------------------------------------------------
1082 double Signal::Correlation(const Signal &input2) const {
1084 assert(size() == input2.size()); //les signaux à comparer doivent avoir la même taille
1086 double cc11,cc22,cc;
1087 Signal input1(*this);
1089 cc=this->UnormalizedCorrelation(input2);
1090 cc11=this->UnormalizedCorrelation(input1);
1091 cc22=input2.UnormalizedCorrelation (input2);
1092 cc=cc/sqrt(cc11*cc22);
1093 std::cout<<"coefficient de corrélation : "<<cc<<std::endl;
1101 //---------------------------------------------------------------------
1102 //---------------------------------------------------------------------
1104 Signal Signal::Xcov(const Signal &input2 , int &shift, double &coef, double &cc) {
1106 Signal temp_output,output;
1112 temp_output.resize(max);
1114 Signal input1(*this);
1116 cc=this->UnormalizedXcov(input2,output);
1117 cc11=this->UnormalizedCorrelation(input1);
1118 cc22=input2.UnormalizedCorrelation (input2);
1120 unsigned int tailleOutput;
1121 tailleOutput=output.size();
1122 for (unsigned int k= 0;k<tailleOutput;k++){
1123 output[k]=output[k]/sqrt(cc11 * cc22 );
1128 for (unsigned int k= 0;k<tailleOutput;k++){
1134 cc=this->Correlation(input2);
1136 std::cout<<"shift: "<<shift<<std::endl;
1137 std::cout<<"XcovMax: "<<coef<<std::endl;
1143 //---------------------------------------------------------------------
1144 //---------------------------------------------------------------------
1147 // double Signal::Compare(Signal & sigRef) {
1148 // double coeffCorrParam[5]={0.,0.,0.,0.,0.};
1150 // const_iterator itSig=begin();
1151 // const_iterator itSigRef=sigRef.begin();//+offset;
1152 // while(itSig!=end()){
1153 // coeffCorrParam[0]+=*itSig;
1154 // coeffCorrParam[1]+=*itSigRef;
1155 // coeffCorrParam[2]+=(*itSig)*(*itSigRef);
1156 // coeffCorrParam[3]+=(*itSig)*(*itSig);
1157 // coeffCorrParam[4]+=(*itSigRef)*(*itSigRef);
1158 // itSig++;itSigRef++;
1161 // double coeffCorr = pow(size()*coeffCorrParam[2]-coeffCorrParam[0]*coeffCorrParam[1],2)/
1162 // ((size()*coeffCorrParam[3]-pow(coeffCorrParam[0],2))*
1163 // (size()*coeffCorrParam[4]-pow(coeffCorrParam[1],2)));
1165 // return coeffCorr;
1168 // int Signal::DerivateSigne( const_iterator & it) const{
1172 // if((*it)==(*(it+pos)))
1174 // else if((*it)<(*(it+pos)))
1176 // else // Case : ((*it)>(*(it+pos)))
1180 // void Signal::CenteredFiniteDifferences(Signal & result,int order,int* weights){
1181 // const_iterator itSig=begin()+order;
1182 // result.resize(size());
1183 // iterator itDer=result.begin()+order;
1184 // while(itSig!=end()-order){
1186 // for(int i=-order;i<=order;i++){
1187 // *itDer+=*(itSig-i)*weights[i+order];
1193 // void Signal::FirstDerivate(Signal & result,int order){
1195 // int weights[3]={-1,0,1};
1196 // CenteredFiniteDifferences(result,order,weights);
1198 // else if(order==2){
1199 // int weights[5]={1,-8,0,8,-1};
1200 // CenteredFiniteDifferences(result,order,weights);
1204 // void Signal::SecondDerivate(Signal & result,int order){
1206 // int weights[3]={1,-2,1};
1207 // CenteredFiniteDifferences(result,order,weights);
1209 // else if(order==2){
1210 // int weights[5]={-1,16,-30,16,-1};
1211 // CenteredFiniteDifferences(result,order,weights);
1217 // void Signal::NormalizeMeanStdDev(double newMean,double newStdDev){
1218 // iterator itSig=begin();
1219 // double sum=0, sum2=0;
1220 // while(itSig!=end()){
1222 // sum2 += (*itSig) * (*itSig);
1225 // double oldMean=sum/size();
1226 // double oldStdDev=sqrt(sum2/size()-oldMean*oldMean);
1228 // double a = newStdDev/oldStdDev;
1229 // double b = newMean - a * oldMean;
1231 // while(itSig!=end()){
1232 // *itSig = a *(*itSig) + b;
1239 // void Signal::print(ostream & os, const int level) const {
1240 // os << "Size:" << m_Data.size() << endl;
1241 // const_iterator it=m_Data.begin();
1242 // while(it!=m_Data.end()){
1243 // os << *it << endl;
1251 // // istream& Signal::get(istream& is) {
1252 // // ERROR << "Signal::get NOT IMPLEMENTED";
1257 // // /** @b GridBase::put os
1260 // // ***************************************************/
1261 // // ostream& Signal::put(ostream& os) const {
1268 // void Signal::Crop(unsigned int posmin, unsigned int posmax){
1269 // if(posmin >= m_Data.size()) return;
1270 // if(posmax >= m_Data.size()) posmax=m_Data.size();
1271 // m_Data.erase(m_Data.begin()+posmax+1,m_Data.end());
1272 // m_Data.erase(m_Data.begin(),m_Data.begin()+posmin);
1275 // void Signal::LinearResample(const unsigned int newSize){
1277 // newData.push_back(front());
1278 // double posInOld,leftWeight,rightWeight;
1279 // int leftPos, rightPos;
1280 // for(unsigned int i=1 ; i < newSize-1 ; i++){
1281 // posInOld = (double)(i * (size()-1)) / (double)(newSize-1);
1282 // leftPos = (int)floor(posInOld);
1283 // rightPos = leftPos+1;
1284 // leftWeight = (double)rightPos - posInOld;
1285 // rightWeight = posInOld - (double)leftPos;
1286 // newData.push_back(m_Data[leftPos] * leftWeight + m_Data[rightPos] * rightWeight );
1289 // newData.push_back(back());
1294 // int Signal::FreqToSamp(double freq){
1295 // if(m_SamplingPeriod==-1.)
1296 // cerr << "ERROR: you did not initialize the sampling period" << endl;
1297 // return lrint(freq*(double)size()*m_SamplingPeriod);
1299 // double Signal::SampToFreq(int samp){
1300 // if(m_SamplingPeriod==-1.)
1301 // cerr << "ERROR: you did not initialize the sampling period" << endl;
1302 // return ((double)(samp)/((double)size()*m_SamplingPeriod));
1305 // //---------------------------------------------------------------------
1306 // Signal Signal::limPhaseDE(eIPhaseValue, eEPhaseValue)
1309 // //Create an empty signal
1310 // phase.resize(size());
1312 // iterator phaseIt=initialPhaseValues.begin();
1313 // iterator monPhaseIt=monPhase.begin();
1314 // iterator extremaIt =begin();
1316 // while (extremaIt!= end())
1318 // if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
1319 // else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
1320 // extremaIt++; phaseIt++;monPhaseIt++;
1328 #endif //#define CLITKSIGNAL