1 | #include "SpikeRemoval.h"
|
---|
2 |
|
---|
3 | #include <iostream>
|
---|
4 |
|
---|
5 | // compute the mean of the left and right neighbors of a channel
|
---|
6 | vector<float> * computeN1mean( vector<float> &src)
|
---|
7 | {
|
---|
8 | vector<float> * dest = new vector<float>;
|
---|
9 |
|
---|
10 | dest->push_back(src[1]);
|
---|
11 |
|
---|
12 | for(unsigned int i = 1; i < src.size() - 1; i++){
|
---|
13 | /* if (i == 0){ // use right sample es mean
|
---|
14 | N1mean[i] = Ameas[i+1];
|
---|
15 | }
|
---|
16 | else if ( i == Samples-1 ){ //use left sample as mean
|
---|
17 | N1mean[i] = Ameas[i-1];
|
---|
18 | }
|
---|
19 | else{
|
---|
20 | N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
|
---|
21 | }
|
---|
22 | */
|
---|
23 | dest->push_back( ( src[i-1] + src[i+1] ) / 2. );
|
---|
24 | }
|
---|
25 |
|
---|
26 | dest->push_back(src[src.size() - 1]);
|
---|
27 | return dest;
|
---|
28 | } // end of computeN1mean computation
|
---|
29 |
|
---|
30 | void removeSpikes(
|
---|
31 | vector<float> &src,
|
---|
32 | vector<float> &dest,
|
---|
33 | const float CandidateTHR,
|
---|
34 | const float nextDiffTHR,
|
---|
35 | const float nextNextDiffTHR
|
---|
36 | ){
|
---|
37 | vector<float> * NextNeighborMean = computeN1mean(src);
|
---|
38 |
|
---|
39 | if (src.size() != NextNeighborMean->size())
|
---|
40 | //TODO .. if verbositylevel .. say something, use the return code to tell the caller..
|
---|
41 | dest.clear();
|
---|
42 | dest = src;
|
---|
43 |
|
---|
44 | float diff, nextDiff, nextNextDiff;
|
---|
45 |
|
---|
46 | // find the spike and replace it by mean value of neighbors
|
---|
47 | for (unsigned int i = 0; i < src.size(); i++) {
|
---|
48 |
|
---|
49 | diff = src[i] - (*NextNeighborMean)[i];
|
---|
50 |
|
---|
51 | if ( diff < CandidateTHR ){ // a spike candidate
|
---|
52 | // check consistency with a single channel spike
|
---|
53 |
|
---|
54 | if ( src[i+2] - ( src[i] + src[i+3] )/2. > 10. )
|
---|
55 | {
|
---|
56 | dest[i+1] = ( src[i] + src[i+3] )/2.;
|
---|
57 | dest[i+2] = ( src[i] + src[i+3] )/2.;
|
---|
58 | i = i + 3;
|
---|
59 | }
|
---|
60 | else
|
---|
61 | {
|
---|
62 | nextDiff = src[i+1] - (*NextNeighborMean)[i+1];
|
---|
63 | nextNextDiff = src[i+2] - (*NextNeighborMean)[i+2];
|
---|
64 |
|
---|
65 | if ( ( nextDiff > nextDiffTHR * diff ) && ( nextNextDiff < nextNextDiffTHR ) ){
|
---|
66 | dest[i+1] = (*NextNeighborMean)[i+1];
|
---|
67 | (*NextNeighborMean)[i+2] = (src[i+1] - src[i+3] / 2.);
|
---|
68 | i = i + 2;//do not care about the next sample it was the spike
|
---|
69 | }
|
---|
70 | // treatment for the end of the pipeline must be added !!!
|
---|
71 | }
|
---|
72 | }
|
---|
73 | } // end of spike search and correction
|
---|
74 | delete NextNeighborMean;
|
---|
75 | return;
|
---|
76 | }
|
---|