Skip to content

Adding ability to ignore bad channels, add a PCA based "select vals" option #24

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

Adding ability to ignore bad channels, add a PCA based "select vals" option #24

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
76ee4b6
This implements the PCA approach to signal protection
SFBayLaser Apr 23, 2025
3af54bb
If an entire channel is considered "bad" then zero the waveform (not …
SFBayLaser Apr 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 37 additions & 45 deletions icarus_signal_processing/Denoising.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,56 +85,38 @@ void icarus_signal_processing::Denoising::getSelectValsWPCA(ArrayFloat::const_it
{
auto nTicks = morphedWaveformsItr[0].size();

// Set a protection width
int halfWidth = std::max(int(window/8),4);

// std::cout << "***** getSelectVals ***** numChannels: " << numChannels << ", grouping: " << grouping << std::endl;

for (size_t i=0; i<numChannels; ++i)
{
VectorFloat localVec = morphedWaveformsItr[i];
Eigen::Vector<float,2> meanPos;
Eigen::Matrix<float,2,2> eigenVectors;
Eigen::Vector<float,2> eigenValues;
float nSigma;

float median = getMedian(localVec.begin(), localVec.end());
// need to create a "point cloud"
PointCloud<float> morphedCloud(nTicks);

std::vector<float> baseVec;
baseVec.resize(localVec.size());
for (size_t channelIdx=0; channelIdx<numChannels; ++channelIdx)
{
VectorFloat localVec = morphedWaveformsItr[channelIdx];

for (size_t j=0; j<baseVec.size(); ++j) baseVec[j] = morphedWaveformsItr[i][j] - median;
for(size_t idx = 0; idx < nTicks; idx++) morphedCloud[idx] = WavePoint<float>(idx,localVec[idx]);

float threshold = thresholdVec[i];
// Now use the PCA to compute correction
getIteratedPCAAxis(morphedCloud,meanPos,eigenVectors,eigenValues,nSigma);

// make sure the roi vector is initialized
std::fill(roiItr[i].begin(),roiItr[i].end(),false);
for (size_t j=0; j<nTicks; ++j) morphedCloud[j][1] = morphedWaveformsItr[channelIdx][j] - meanPos[1];

int lb(-2 * halfWidth);
meanPos[1] = 0.;

for (size_t j=0; j<nTicks; ++j)
//float threshold = std::max(float(5.)*std::sqrt(eigenValues[0]),thresholdVec[channelIdx]);
float threshold = thresholdVec[channelIdx] * std::sqrt(eigenValues[0]);

for(size_t tickIdx = 0; tickIdx < nTicks; tickIdx++)
{
// Are we over threshold?
if (std::fabs(baseVec[j]) > threshold)
{
// Check Bounds
selectValsItr[i][j] = true;

// Check ROI limits
if (lb < -halfWidth) lb = j - halfWidth;
}
else
{
selectValsItr[i][j] = false;

// Check if we had gone over threshold and need to set our ROI
if (lb > -(halfWidth+1))
{
int ub = j + halfWidth;
size_t lowerBound = std::max(lb, 0);
size_t upperBound = std::min(ub, (int) nTicks);
float cutVal = (eigenVectors * (morphedCloud[tickIdx] - meanPos))[0] + threshold;

for (size_t k=lowerBound; k<upperBound; ++k) roiItr[i][k] = true;
}

lb = -2 * halfWidth;
}
if (std::abs(morphedCloud[tickIdx][1]) > cutVal) selectValsItr[channelIdx][tickIdx] = true;
}
}

Expand Down Expand Up @@ -190,7 +172,9 @@ void icarus_signal_processing::Denoiser1D::operator()(ArrayFloat::iterator
selStart = morphStop;
}

getSelectVals(morphedWaveformsItr, selectValsItr, roiItr, thresholdVec, numChannels, grouping, window);
// Attempt to protect signal
//getSelectVals(morphedWaveformsItr, selectValsItr, roiItr, thresholdVec, numChannels, grouping, window);
getSelectValsWPCA(morphedWaveformsItr, selectValsItr, roiItr, thresholdVec, numChannels, grouping, window);

if (fOutputStats)
{
Expand Down Expand Up @@ -1862,6 +1846,12 @@ void icarus_signal_processing::Denoising::removeCoherentNoise(ArrayFloat::iterat
// get an instance of the waveform tools
icarus_signal_processing::WaveformTools<float> waveformTools;

// identify bad channels
std::vector<bool> isGoodChannel(numChannels);

for(size_t channelIdx = 0; channelIdx < numChannels; channelIdx++)
isGoodChannel[channelIdx] = std::count(selectValsItr[channelIdx].begin(),selectValsItr[channelIdx].end(),true) < 4096 ? true : false;

VectorFloat vl(grouping/2);

VectorFloat vu(grouping/2);
Expand All @@ -1884,11 +1874,11 @@ void icarus_signal_processing::Denoising::removeCoherentNoise(ArrayFloat::iterat
// Allow for signal protection?
// When subtracting in groups of 64 this does not seem to be working as one might expect
// So removing for now (12/6/2021)
// if (!selectValsItr[c][i])
// {
if (!selectValsItr[c][i])
{
if (c < group_mid) vl[idxL++] = filteredWaveformsItr[c][i];
else vu[idxU++] = filteredWaveformsItr[c][i];
// }
}
}

float median(0.);
Expand Down Expand Up @@ -1960,7 +1950,8 @@ void icarus_signal_processing::Denoising::removeCoherentNoise(ArrayFloat::iterat
for (auto k=group_start; k<group_end; ++k)
{
correctedMediansItr[k][i] = median;
waveLessCoherentItr[k][i] = filteredWaveformsItr[k][i] - median;
if (isGoodChannel[k]) waveLessCoherentItr[k][i] = filteredWaveformsItr[k][i] - median;
else waveLessCoherentItr[k][i] = 0.;
}
}
}
Expand Down Expand Up @@ -2312,8 +2303,9 @@ void icarus_signal_processing::Denoising::removeCoherentNoise(ArrayFloat::iterat
{
float median = correctedMediansItr[k][j];

waveLessCoherentItr[k][j] = filteredWaveformsItr[k][j] - median;
v[idxV++] = waveLessCoherentItr[k][j];
// Don't apply correction if channel labeled "bad"
if (!selectValsItr[k][j]) waveLessCoherentItr[k][j] = filteredWaveformsItr[k][j] - median;
v[idxV++] = waveLessCoherentItr[k][j];
}

intrinsicRMSItr[i][j] = std::sqrt(std::inner_product(v.begin(), v.begin()+idxV, v.begin(), 0.) / float(v.size()));
Expand Down