1 #ifndef CLITKSIGNAL_CXX
2 #define CLITKSIGNAL_CXX
4 #include "clitkSignal.h"
8 //---------------------------------------------------------------------
9 void Signal::Read(string fileName){
10 ifstream signalStream(fileName.c_str());
11 SignalValueType point;
12 if(!signalStream.is_open()){
13 std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
16 skipComment(signalStream);
17 while(!signalStream.eof()) {
18 skipComment(signalStream);
19 signalStream >> point;
20 skipComment(signalStream);
21 m_Data.push_back(point);
22 // cout << point << endl;
25 SetSamplingPeriod(1.);
27 //---------------------------------------------------------------------
30 //---------------------------------------------------------------------
31 void Signal::Read(string fileName, int col){
32 ifstream signalStream(fileName.c_str());
33 SignalValueType point;
34 if(!signalStream.is_open()){
35 std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
38 skipComment(signalStream);
39 while(!signalStream.eof()) {
40 skipComment(signalStream);
44 std::getline(signalStream, line);
47 istringstream iss(line);
48 for(int i=0; i<col; i++) { // skip col-1 first column
54 std::cerr << "ERROR: no col n" << col << " in the line '" << line << "' ?" << std::endl;
57 skipComment(signalStream);
58 m_Data.push_back(point);
59 // cout << point << endl;
62 SetSamplingPeriod(1.);
64 //---------------------------------------------------------------------
68 //---------------------------------------------------------------------
69 //Convert 1D image to signal
70 Signal Signal::ConvertImageToSignal( Signal::ImageType::Pointer image)
75 //make an image iterator
76 itk::ImageRegionConstIterator<ImageType> it(image,image->GetLargestPossibleRegion());
82 signal.push_back(it.Get());
87 signal.SetSamplingPeriod(image->GetSpacing()[0]);
91 //---------------------------------------------------------------------
94 //---------------------------------------------------------------------
95 //Convert 1D signal to image
96 Signal::ImageType::Pointer Signal::ConvertSignalToImage(Signal signal)
99 ImageType::Pointer image =ImageType::New();
100 ImageType::RegionType region;
101 ImageType::RegionType::IndexType index;
103 ImageType::RegionType::SizeType size;
104 size[0]=signal.size();
106 region.SetIndex(index);
107 region.SetSize(size);
109 image->SetRegions(region);
112 //make an image iterator
113 itk::ImageRegionIterator<ImageType> mIt(image,image->GetLargestPossibleRegion());
116 //make a signal iterator
117 Signal::const_iterator sIt=signal.begin();
120 while(sIt!=signal.end())
128 ImageType::SpacingType spacing;
129 spacing[0]=signal.GetSamplingPeriod();
130 image->SetSpacing(spacing);
135 //---------------------------------------------------------------------
139 //---------------------------------------------------------------------
140 void Signal::Write(const string fileName){
141 ofstream signalStream(fileName.c_str());
142 if(!signalStream.is_open()){
143 cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Write" << endl;
149 signalStream << *it << endl;
152 signalStream.close();
154 //---------------------------------------------------------------------
157 //---------------------------------------------------------------------
158 Signal & Signal::operator/=(Signal & div){
159 if(size()!= div.size())
161 std::cerr << "Error: signal size must be the same!" << std::endl;
172 cerr << "Division by 0 in operator/= skipped" << endl;
177 //---------------------------------------------------------------------
180 //---------------------------------------------------------------------
181 Signal & Signal::operator*=(Signal & mul){
182 if(size()!= mul.size())
184 std::cerr << "Error: signal size must be the same!" << std::endl;
196 //---------------------------------------------------------------------
199 //---------------------------------------------------------------------
200 Signal Signal::Normalize(double newMin,double newMax){
201 Signal temp (m_SamplingPeriod);
202 vector<double> extrema=GetGlobalMinMax();
203 iterator itSig=begin();
205 temp.push_back( ((*itSig)-extrema[0])*(newMax-newMin)/(extrema[1]-extrema[0]) + newMin );
210 //---------------------------------------------------------------------
213 //---------------------------------------------------------------------
214 vector<double> Signal::GetGlobalMinMax() const {
215 vector<double> extrema(2);
217 cerr << "ERROR: GetExtrema / No signal" << endl;
220 extrema[0]=m_Data[0];
221 extrema[1]=m_Data[0];
222 for(unsigned int i=1;i<m_Data.size();i++){
223 if(extrema[0]>m_Data[i]) extrema[0]=m_Data[i];
224 if(extrema[1]<m_Data[i]) extrema[1]=m_Data[i];
228 //---------------------------------------------------------------------
231 //---------------------------------------------------------------------
232 double Signal::GetGlobalMean() const {
234 for(unsigned int i=0;i<m_Data.size();i++){
237 return m/(double)m_Data.size();
239 //---------------------------------------------------------------------
242 //---------------------------------------------------------------------
243 Signal Signal::MovingAverageFilter ( unsigned int length) {
245 Signal temp(m_SamplingPeriod);
247 for (unsigned int i=0; i <size(); i++)
249 double accumulator=0.;
250 unsigned int scale=0;
251 for (unsigned int j= max(0,static_cast<int>(i)-static_cast<int>(length)); j<=min(size(), i+length); j++)
253 accumulator+=m_Data[j];
256 temp.push_back(accumulator/scale);
260 //---------------------------------------------------------------------
263 //---------------------------------------------------------------------
264 Signal Signal::GaussLikeFilter ( ) {
266 Signal temp(m_SamplingPeriod);
271 //first sample: mirrorring BC
272 temp.push_back((2.*m_Data[0]+2*m_Data[1])/4.);
275 for (unsigned int i=1; i <size()-1; i++)
276 temp.push_back((2.*m_Data[i]+m_Data[i-1]+m_Data[i+1])/4.);
278 //end sample: mirrorring BC
279 temp.push_back((2.*m_Data[size()-1]+2*m_Data[size()-2])/4.);
283 //---------------------------------------------------------------------
286 //---------------------------------------------------------------------
287 Signal Signal::NormalizeMeanStdDev(double newMean,double newStdDev)
289 iterator itSig=begin();
290 double sum=0, sum2=0;
291 Signal temp(m_SamplingPeriod);
293 //Calculate the mean and the std dev
297 sum2 += (*itSig) * (*itSig);
300 double oldMean=sum/size();
301 double oldStdDev=sqrt(sum2/size()-oldMean*oldMean);
304 double a = newStdDev/oldStdDev;
305 double b = newMean - a * oldMean;
309 temp.push_back( a *(*itSig) + b );
314 //---------------------------------------------------------------------
317 //---------------------------------------------------------------------
318 // Signal Signal::HighPassFilter (double sampPeriod, double cutOffFrequency )
321 // Signal temp(m_SamplingPeriod);
322 // temp.resize(size());
325 // SIGNAL_FFT_TYPE fft;
327 // //calculate the cut off frequency
328 // unsigned int samp=lrint(cutOffFrequency*static_cast<double>(size())*sampPeriod);
330 // //forward fft with empty fft
331 // if(fft.size()==0) OneDForwardFourier(*this, fft);
333 // //remove the frequencies
334 // for(unsigned int i=0;i<samp && i<fft.size();i++)
335 // fft[i]=complex<double>(0.,0.);
337 // //backward with remaining frequencies
338 // OneDBackwardFourier(fft,temp);
341 //---------------------------------------------------------------------
344 //---------------------------------------------------------------------
345 // Signal Signal::LowPassFilter (double sampPeriod, double cutOffFrequency )
348 // Signal temp(m_SamplingPeriod);
349 // temp.resize(size());
352 // SIGNAL_FFT_TYPE fft;
354 // //calculate the cut off frequency
355 // unsigned int samp=lrint(cutOffFrequency*static_cast<double>(size())*sampPeriod);
357 // //forward fft with empty fft
358 // if(fft.size()==0) OneDForwardFourier(*this, fft);
359 // unsigned int fsize=fft.size();
361 // //remove the frequencies
362 // unsigned int limit=min (samp, fsize);
363 // for(unsigned int i=limit;i<fft.size();i++)
364 // fft[i]=complex<double>(0.,0.);
366 // //backward with remaining frequencies
367 // OneDBackwardFourier(fft,temp);
370 //---------------------------------------------------------------------
373 //---------------------------------------------------------------------
374 // void Signal::OneDForwardFourier(const Signal& input, SIGNAL_FFT_TYPE & fft)
376 // //Create output array
377 // fft.resize(input.size()/2+1);
379 // double *tempCopy=new double[size()];
380 // copy(begin(), end(), tempCopy);
382 // //Forward Fourier Transform
384 // p=fftw_plan_dft_r2c_1d(size(),tempCopy,reinterpret_cast<fftw_complex*>(&(fft[0])),FFTW_ESTIMATE);
386 // fftw_destroy_plan(p);
387 // //delete tempCopy;
390 //---------------------------------------------------------------------
393 //---------------------------------------------------------------------
394 // void Signal::OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output)
399 // p=fftw_plan_dft_c2r_1d(output.size(),reinterpret_cast<fftw_complex*>(&(fft[0])),&(output[0]),FFTW_ESTIMATE);
401 // fftw_destroy_plan(p);
403 // vector<double>::iterator it=output.begin();
404 // while(it!=output.end()){
405 // *it /= (double)output.size();
410 //---------------------------------------------------------------------
413 //---------------------------------------------------------------------
414 // double Signal::MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft)
417 // if(fft.size()==0) OneDForwardFourier(sig,fft);
419 // double amplitude, amplitudeMax=abs(fft[1]);
420 // for(unsigned int i=1;i<fft.size();i++){
421 // amplitude=abs(fft[i]);
422 // if(amplitude>amplitudeMax){
424 // amplitudeMax=amplitude;
427 // return ((double)(posMax)/((double)sig.size()*sig.GetSamplingPeriod()));
429 //---------------------------------------------------------------------
432 //---------------------------------------------------------------------
433 Signal Signal::DetectLocalExtrema(unsigned int width)
435 Signal temp(m_SamplingPeriod);
437 unsigned int upper, lower;
439 //has to be at least 1
440 width=max(static_cast<int>(width), 1);
442 for(unsigned int i=0 ; i < size() ; i++){
446 for(unsigned int j=1; j< width+1; j++)
449 upper = min( size(), i+j);
450 lower = max( static_cast<int>(0), (int)i-(int)j);
453 if( ! (m_Data[i] >= m_Data[lower] && m_Data[i] >= m_Data[upper]))isMax=false;
456 if( ! (m_Data[i]<= m_Data[lower] && m_Data[i] <= m_Data[upper])) isMin=false;
458 //if neither, go to the next value
459 if( (!isMax) && (!isMin)) break;
462 if (isMax) temp.push_back(1);
463 else if (isMin) temp.push_back(0);
464 else temp.push_back(0.5);
468 //---------------------------------------------------------------------
471 //---------------------------------------------------------------------
472 Signal Signal::LimPhase()
475 //phase is defined as going from 0 to 1 linearly between end expiration
476 Signal phase(m_SamplingPeriod);
477 phase.resize(size());
479 unsigned int firstBeginCycle=0;
480 unsigned int firstEndCycle=0;
481 unsigned int beginCycle=0;
483 //=========================================
484 //Search the first value in extrema not 0.5
485 while(m_Data[beginCycle]!=1)
490 //We search the corresponding end
491 unsigned int endCycle=beginCycle+1;
492 while(endCycle != size() && m_Data[endCycle]!=1){
497 //============================================================
498 //Calculate phase at the beginning (before the first extremum)
499 for(unsigned int i=0 ; i<beginCycle ; i++)
502 //if before first cycle is longer than first cycle: let it go from 0 to 1
503 if((double)beginCycle/(double)(endCycle-beginCycle) > 1)
505 phase[i] = (double)(i-0)/(double)(beginCycle-0);
507 //copy the phase values later
510 firstBeginCycle=beginCycle;
511 firstEndCycle=endCycle;
516 //===================================================================
518 while(endCycle != size()){
520 //fill between extrema
521 for(unsigned int i=beginCycle ; i<endCycle ; i++){
522 phase[i] = (double)(i-beginCycle)/(double)(endCycle-beginCycle);
526 beginCycle = endCycle++;
527 while(endCycle != size() && m_Data[endCycle]!=1)
531 //===================================================================
532 //Calculate phase at the end (after the last extremum)
533 endCycle = beginCycle--;
535 //count backwards till the previous
536 while(m_Data[beginCycle]!=1) beginCycle--;
537 for(unsigned int i=endCycle ; i<size() ; i++)
540 //after last extrema is longer then last cycle
541 if((double)(size()-endCycle)/(double)(endCycle-beginCycle) > 1){
543 //make the last part go till 1
544 phase[i] = (double)(i-endCycle)/(double)(size()-endCycle);
547 //the last part is shorter, copy the last cycle values
549 phase[i] = phase[i -(endCycle-beginCycle)];
553 //===================================================================
554 //check it some remains to be copied in the beginning
555 if (firstBeginCycle!=0)
557 for(unsigned int i=0 ; i<firstBeginCycle ; i++)
558 phase[i] =phase[i+firstEndCycle-firstBeginCycle];
563 //---------------------------------------------------------------------
567 //---------------------------------------------------------------------
568 Signal Signal::MonPhase()
571 //phase is defined as going from 0 to 1 linearly between end expiration
572 Signal phase(m_SamplingPeriod);
573 phase.resize(size());
575 unsigned int firstBeginCycle=0;
576 unsigned int firstEndCycle=0;
577 unsigned int cycleCounter=0;
578 unsigned int beginCycle=0;
580 //===================================================================
581 //Search the first value in extrema not 0.5
582 while(m_Data[beginCycle]!=1)
587 //We search the corresponding end
588 unsigned int endCycle=beginCycle+1;
589 while(endCycle != size() && m_Data[endCycle]!=1){
594 //===================================================================
595 //Calculate phase at the beginning (before the first extremum)
596 for(unsigned int i=0 ; i<beginCycle ; i++)
599 //if before first cycle is longer than first cycle: let it go from 0 to 1
600 if((double)beginCycle/(double)(endCycle-beginCycle) > 1)
602 phase[i] = (double)(i-0)/(double)(beginCycle-0);
604 //copy the phase values later
607 firstBeginCycle=beginCycle;
608 firstEndCycle=endCycle;
613 //===================================================================
615 while(endCycle != size()){
618 //fill between extrema
619 for(unsigned int i=beginCycle ; i<endCycle ; i++)
620 phase[i] = (double)cycleCounter+(double)(i-beginCycle)/(double)(endCycle-beginCycle);
624 beginCycle = endCycle++;
625 while(endCycle != size() && m_Data[endCycle]!=1)
629 //===================================================================
630 //Calculate phase at the end (after the last extremum)
631 endCycle = beginCycle--;
633 //count backwards till the previous
634 while(m_Data[beginCycle]!=1) beginCycle--;
635 for(unsigned int i=endCycle ; i<size() ; i++)
638 //after last extrema is longer then last cycle
639 if((double)(size()-endCycle)/(double)(endCycle-beginCycle) > 1){
641 //make the last part go till 1
642 phase[i] = (double)cycleCounter+(double)(i-endCycle)/(double)(size()-endCycle);
645 //the last part is shorter, copy the last cycle values
647 phase[i] = phase[i -(endCycle-beginCycle)]+1;
651 //===================================================================
652 //check it some remains to be copied in the beginning
653 if (firstBeginCycle!=0)
655 for(unsigned int i=0 ; i<firstBeginCycle ; i++)
656 phase[i] =phase[i+firstEndCycle-firstBeginCycle]-1;
661 //---------------------------------------------------------------------
664 //---------------------------------------------------------------------
665 Signal Signal::MonPhaseDE(double eEPhaseValue, double eIPhaseValue)
667 //First calculate monPhase
668 Signal monPhase=this->MonPhase();
670 //Create an empty signal
671 Signal phase(size(), -1);
672 phase.SetSamplingPeriod(m_SamplingPeriod);
674 //Fill in the values at the extrema position
675 iterator phaseIt=phase.begin();
676 iterator monPhaseIt=monPhase.begin();
677 iterator extremaIt =begin();
678 while (extremaIt!= end())
680 if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
681 else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
682 extremaIt++; phaseIt++;monPhaseIt++;
688 //---------------------------------------------------------------------
691 //---------------------------------------------------------------------
692 Signal Signal::LinearlyInterpolateScatteredValues()
694 //Linearly interpolate the values in between
696 unsigned int beginCycle=0;
697 unsigned int endCycle=0;
699 //Create a new signal
700 Signal temp(size(),-1);
701 temp.SetSamplingPeriod(m_SamplingPeriod);
703 //start from the first value
704 while (m_Data[i]==-1)i++;
709 while ( (m_Data[i]==-1) && (i<size()) )i++;
715 for (unsigned int k=beginCycle;k<endCycle+1; k++)
716 temp[k]=m_Data[beginCycle]+(double)(k-beginCycle)*
717 (m_Data[endCycle]-m_Data[beginCycle])/(double)(endCycle-beginCycle);
722 while( (i< size()) && (m_Data[i]==-1) ) i++;
726 if (beginCycle!= size()-1)
728 //continue with same slope
729 for (unsigned int k=beginCycle+1; k<size();k++)
730 temp[k]=temp[beginCycle]+(double)(k-beginCycle)*(temp[beginCycle]-temp[beginCycle-1]);
738 while (temp[i]==-1)i++;
741 while (m_Data[i]==-1)i++;
744 //if the first filled half cycle is longer, copy it
745 if(beginCycle<(endCycle-beginCycle))
747 for (unsigned int k=0; k< beginCycle;k++)
748 temp[k]=temp[k+endCycle-beginCycle];
751 //if the first filled half cycle is longer, start from zero
754 for (unsigned int k=0; k< beginCycle;k++)
755 temp[k]=k*temp[beginCycle]/(beginCycle);
761 //---------------------------------------------------------------------
764 // //---------------------------------------------------------------------
765 // Signal Signal::ApproximateScatteredValuesWithBSplines(unsigned int splineOrder, unsigned int numberOfControlPoints)
767 // //filter requires a vector pixelType
768 // typedef itk::PointSet<VectorType, 1> PointSetType;
769 // PointSetType::Pointer pointSet=PointSetType::New();
770 // typedef PointSetType::PointType PointType;
773 // unsigned int pointIndex=0;
779 // p[0]= i;//JV spacing is allways 1
780 // pointSet->SetPoint( pointIndex, p );
781 // pointSet->SetPointData( pointIndex, m_Data[i] );
787 // //define the output signal properties
788 // ImageType::RegionType::SizeType outputSize;
789 // outputSize[0]= size();
790 // ImageType::PointType outputOrigin;
791 // outputOrigin[0]=0.0;//JV may need to be changed
792 // ImageType::SpacingType outputSpacing;
793 // outputSpacing[0]=1; //JV add relation to the original signal spacing
796 // typedef itk::BSplineScatteredDataPointSetToImageFilter< PointSetType, VectorImageType > PointSetToImageFilterType;
797 // PointSetToImageFilterType::Pointer pointSetToImageFilter= PointSetToImageFilterType::New();
798 // pointSetToImageFilter->SetInput(pointSet);
799 // pointSetToImageFilter->SetSplineOrder(splineOrder);//JV
800 // pointSetToImageFilter->SetSize(outputSize);
801 // pointSetToImageFilter->SetOrigin(outputOrigin);
802 // pointSetToImageFilter->SetSpacing(outputSpacing);
805 // itk::FixedArray<unsigned int,1> num;
806 // num[0]=numberOfControlPoints;
807 // pointSetToImageFilter->SetNumberOfControlPoints(num);//JV
808 // pointSetToImageFilter->Update();
809 // VectorImageType::Pointer approximatedSignal=pointSetToImageFilter->GetOutput();
811 // //Convert and return
812 // return ConvertVectorImageToSignal(approximatedSignal);
814 // //---------------------------------------------------------------------
817 //---------------------------------------------------------------------
818 Signal Signal::ConvertVectorImageToSignal(VectorImageType::Pointer image)
823 //make an image iterator
824 itk::ImageRegionConstIterator<VectorImageType> it(image,image->GetLargestPossibleRegion());
830 signal.push_back(it.Get()[0]);
835 signal.SetSamplingPeriod(image->GetSpacing()[0]);
839 //---------------------------------------------------------------------
842 //---------------------------------------------------------------------
843 Signal Signal::LimitSignalRange()
846 Signal signal(m_SamplingPeriod);
850 signal.push_back(*it-floor(*it));
857 //---------------------------------------------------------------------
858 double Signal::SSD(const Signal &sig2) const{
859 if(sig2.size() != size()){
860 cerr << "ERROR in Signal::SSD: signals don't have the same size" << endl;
864 for(unsigned int i=0;i<size();i++){
865 result+=pow(sig2[i]-m_Data[i],2);
871 //---------------------------------------------------------------------
874 //---------------------------------------------------------------------
875 void Signal::AddValue(double v) {
876 for(unsigned int i=0;i<size();i++) {
880 //---------------------------------------------------------------------
883 //---------------------------------------------------------------------
884 void Signal::ComputeAugmentedSpace(clitk::Signal & outputX,
885 clitk::Signal & outputY,
886 unsigned int delay) const {
887 if (size() <= delay) {
888 std::cerr << "Error in signal length is " << size()
889 << " while delay is " << delay << " : too short. Abort." << std::endl;
892 outputX.resize(size()-delay);
893 outputY.resize(size()-delay);
894 for(unsigned int i=0; i<outputX.size(); i++) {
895 outputX[i] = m_Data[i+delay];
896 outputY[i] = m_Data[i];
899 //---------------------------------------------------------------------
914 // double Signal::Compare(Signal & sigRef) {
915 // double coeffCorrParam[5]={0.,0.,0.,0.,0.};
917 // const_iterator itSig=begin();
918 // const_iterator itSigRef=sigRef.begin();//+offset;
919 // while(itSig!=end()){
920 // coeffCorrParam[0]+=*itSig;
921 // coeffCorrParam[1]+=*itSigRef;
922 // coeffCorrParam[2]+=(*itSig)*(*itSigRef);
923 // coeffCorrParam[3]+=(*itSig)*(*itSig);
924 // coeffCorrParam[4]+=(*itSigRef)*(*itSigRef);
925 // itSig++;itSigRef++;
928 // double coeffCorr = pow(size()*coeffCorrParam[2]-coeffCorrParam[0]*coeffCorrParam[1],2)/
929 // ((size()*coeffCorrParam[3]-pow(coeffCorrParam[0],2))*
930 // (size()*coeffCorrParam[4]-pow(coeffCorrParam[1],2)));
935 // int Signal::DerivateSigne( const_iterator & it) const{
939 // if((*it)==(*(it+pos)))
941 // else if((*it)<(*(it+pos)))
943 // else // Case : ((*it)>(*(it+pos)))
947 // void Signal::CenteredFiniteDifferences(Signal & result,int order,int* weights){
948 // const_iterator itSig=begin()+order;
949 // result.resize(size());
950 // iterator itDer=result.begin()+order;
951 // while(itSig!=end()-order){
953 // for(int i=-order;i<=order;i++){
954 // *itDer+=*(itSig-i)*weights[i+order];
960 // void Signal::FirstDerivate(Signal & result,int order){
962 // int weights[3]={-1,0,1};
963 // CenteredFiniteDifferences(result,order,weights);
965 // else if(order==2){
966 // int weights[5]={1,-8,0,8,-1};
967 // CenteredFiniteDifferences(result,order,weights);
971 // void Signal::SecondDerivate(Signal & result,int order){
973 // int weights[3]={1,-2,1};
974 // CenteredFiniteDifferences(result,order,weights);
976 // else if(order==2){
977 // int weights[5]={-1,16,-30,16,-1};
978 // CenteredFiniteDifferences(result,order,weights);
984 // void Signal::NormalizeMeanStdDev(double newMean,double newStdDev){
985 // iterator itSig=begin();
986 // double sum=0, sum2=0;
987 // while(itSig!=end()){
989 // sum2 += (*itSig) * (*itSig);
992 // double oldMean=sum/size();
993 // double oldStdDev=sqrt(sum2/size()-oldMean*oldMean);
995 // double a = newStdDev/oldStdDev;
996 // double b = newMean - a * oldMean;
998 // while(itSig!=end()){
999 // *itSig = a *(*itSig) + b;
1006 // void Signal::print(ostream & os, const int level) const {
1007 // os << "Size:" << m_Data.size() << endl;
1008 // const_iterator it=m_Data.begin();
1009 // while(it!=m_Data.end()){
1010 // os << *it << endl;
1018 // // istream& Signal::get(istream& is) {
1019 // // ERROR << "Signal::get NOT IMPLEMENTED";
1024 // // /** @b GridBase::put os
1027 // // ***************************************************/
1028 // // ostream& Signal::put(ostream& os) const {
1035 // void Signal::Crop(unsigned int posmin, unsigned int posmax){
1036 // if(posmin >= m_Data.size()) return;
1037 // if(posmax >= m_Data.size()) posmax=m_Data.size();
1038 // m_Data.erase(m_Data.begin()+posmax+1,m_Data.end());
1039 // m_Data.erase(m_Data.begin(),m_Data.begin()+posmin);
1042 // void Signal::LinearResample(const unsigned int newSize){
1044 // newData.push_back(front());
1045 // double posInOld,leftWeight,rightWeight;
1046 // int leftPos, rightPos;
1047 // for(unsigned int i=1 ; i < newSize-1 ; i++){
1048 // posInOld = (double)(i * (size()-1)) / (double)(newSize-1);
1049 // leftPos = (int)floor(posInOld);
1050 // rightPos = leftPos+1;
1051 // leftWeight = (double)rightPos - posInOld;
1052 // rightWeight = posInOld - (double)leftPos;
1053 // newData.push_back(m_Data[leftPos] * leftWeight + m_Data[rightPos] * rightWeight );
1056 // newData.push_back(back());
1061 // int Signal::FreqToSamp(double freq){
1062 // if(m_SamplingPeriod==-1.)
1063 // cerr << "ERROR: you did not initialize the sampling period" << endl;
1064 // return lrint(freq*(double)size()*m_SamplingPeriod);
1066 // double Signal::SampToFreq(int samp){
1067 // if(m_SamplingPeriod==-1.)
1068 // cerr << "ERROR: you did not initialize the sampling period" << endl;
1069 // return ((double)(samp)/((double)size()*m_SamplingPeriod));
1072 // //---------------------------------------------------------------------
1073 // Signal Signal::limPhaseDE(eIPhaseValue, eEPhaseValue)
1076 // //Create an empty signal
1077 // phase.resize(size());
1079 // iterator phaseIt=initialPhaseValues.begin();
1080 // iterator monPhaseIt=monPhase.begin();
1081 // iterator extremaIt =begin();
1083 // while (extremaIt!= end())
1085 // if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
1086 // else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
1087 // extremaIt++; phaseIt++;monPhaseIt++;
1095 #endif //#define CLITKSIGNAL