2025/8/26 commit
This commit is contained in:
parent
e520579c99
commit
bd389eef55
|
|
@ -221,3 +221,4 @@ g++ -o test_motion_artifact test_motion_artifact.cpp
|
|||
- 质量改善：SNR提升8-15dB
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
4
main.cpp
4
main.cpp
|
|
@ -2,7 +2,7 @@
|
|||
#include <fstream>
|
||||
#include <sstream>
|
||||
#include <ctime>
|
||||
|
||||
#include <windows.h> // 引入 Windows API 头文件
|
||||
|
||||
std::vector<float> heart_rate;
|
||||
|
||||
|
|
@ -360,7 +360,7 @@ void process_ppg_a50_data() {
|
|||
|
||||
int main() {
|
||||
|
||||
|
||||
SetConsoleOutputCP(CP_UTF8);
|
||||
// 选择要运行的测试
|
||||
std::cout << "请选择测试模式:" << std::endl;
|
||||
std::cout << "1. 测试MIT-BIH数据处理" << std::endl;
|
||||
|
|
|
|||
|
|
@ -1067,299 +1067,43 @@ float SignalProcessor::calculate_PPG_sqi(const std::vector<float>& red_channel,
|
|||
std::vector<float> SignalProcessor::remove_motion_artifacts(const std::vector<float>& signal, double sample_rate) {
|
||||
if (signal.empty()) return signal;
|
||||
|
||||
std::cout << "开始增强版运动伪迹检测和去除，信号长度: " << signal.size() << " 样本" << std::endl;
|
||||
|
||||
std::vector<float> result = signal;
|
||||
const size_t min_window_size = static_cast<size_t>(0.1 * sample_rate); // 最小窗口0.1秒
|
||||
const size_t max_window_size = static_cast<size_t>(2.0 * sample_rate); // 最大窗口2秒
|
||||
const size_t window_size = static_cast<size_t>(0.5 * sample_rate); // 0.5秒窗口
|
||||
|
||||
if (signal.size() < min_window_size) {
|
||||
std::cout << "信号长度不足，跳过运动伪迹检测" << std::endl;
|
||||
return result;
|
||||
if (signal.size() < window_size) return result;
|
||||
|
||||
// 使用滑动窗口检测运动伪迹
|
||||
for (size_t i = window_size; i < signal.size(); ++i) {
|
||||
// 计算窗口内的统计特性
|
||||
float sum = 0.0f, sum_sq = 0.0f;
|
||||
for (size_t j = i - window_size; j < i; ++j) {
|
||||
sum += signal[j];
|
||||
sum_sq += signal[j] * signal[j];
|
||||
}
|
||||
|
||||
// 1. 多尺度运动伪迹检测
|
||||
std::vector<bool> artifact_mask(signal.size(), false);
|
||||
float mean = sum / window_size;
|
||||
float variance = (sum_sq / window_size) - (mean * mean);
|
||||
float std_dev = std::sqrt(std::max(0.0f, variance));
|
||||
|
||||
// 使用多个窗口大小进行检测
|
||||
std::vector<size_t> window_sizes = {
|
||||
static_cast<size_t>(0.1 * sample_rate), // 0.1秒 - 快速运动
|
||||
static_cast<size_t>(0.5 * sample_rate), // 0.5秒 - 中等运动
|
||||
static_cast<size_t>(1.0 * sample_rate), // 1.0秒 - 慢速运动
|
||||
static_cast<size_t>(2.0 * sample_rate) // 2.0秒 - 长期漂移
|
||||
};
|
||||
|
||||
std::cout << "使用多尺度检测窗口: ";
|
||||
for (size_t ws : window_sizes) {
|
||||
std::cout << ws << " ";
|
||||
}
|
||||
std::cout << "样本" << std::endl;
|
||||
|
||||
// 2. 统计特征检测
|
||||
for (size_t ws : window_sizes) {
|
||||
if (signal.size() < ws) continue;
|
||||
|
||||
for (size_t i = ws; i < signal.size(); ++i) {
|
||||
// 计算滑动窗口统计特征
|
||||
std::vector<float> window_data;
|
||||
for (size_t j = i - ws; j < i; ++j) {
|
||||
window_data.push_back(signal[j]);
|
||||
}
|
||||
|
||||
// 计算统计特征
|
||||
float mean = 0.0f, variance = 0.0f, skewness = 0.0f, kurtosis = 0.0f;
|
||||
float sum = 0.0f, sum_sq = 0.0f, sum_cube = 0.0f, sum_quad = 0.0f;
|
||||
|
||||
for (float val : window_data) {
|
||||
sum += val;
|
||||
sum_sq += val * val;
|
||||
sum_cube += val * val * val;
|
||||
sum_quad += val * val * val * val;
|
||||
}
|
||||
|
||||
mean = sum / ws;
|
||||
variance = (sum_sq / ws) - (mean * mean);
|
||||
|
||||
if (variance > 1e-6f) {
|
||||
float std_dev = std::sqrt(variance);
|
||||
float normalized_std = std_dev / (std::abs(mean) + 1e-6f);
|
||||
|
||||
// 偏度和峰度计算
|
||||
for (float val : window_data) {
|
||||
float normalized_val = (val - mean) / std_dev;
|
||||
skewness += normalized_val * normalized_val * normalized_val;
|
||||
kurtosis += normalized_val * normalized_val * normalized_val * normalized_val;
|
||||
}
|
||||
skewness /= ws;
|
||||
kurtosis /= ws;
|
||||
|
||||
// 运动伪迹检测条件
|
||||
bool is_artifact = false;
|
||||
|
||||
// 条件1: 异常幅度变化
|
||||
// 检测异常值（可能的运动伪迹）
|
||||
float current_val = signal[i];
|
||||
float z_score = std::abs(current_val - mean) / (std_dev + 1e-6f);
|
||||
if (z_score > 4.0f) { // 提高阈值，减少误检
|
||||
is_artifact = true;
|
||||
|
||||
// 如果Z-score > 3，认为是运动伪迹
|
||||
if (z_score > 3.0f) {
|
||||
// 使用中值滤波替换异常值
|
||||
std::vector<float> window_values;
|
||||
for (size_t j = std::max(static_cast<size_t>(0), i - window_size/2);
|
||||
j < std::min(signal.size(), i + window_size/2); ++j) {
|
||||
window_values.push_back(signal[j]);
|
||||
}
|
||||
|
||||
// 条件2: 异常统计特征
|
||||
if (std::abs(skewness) > 2.0f || kurtosis > 8.0f) {
|
||||
is_artifact = true;
|
||||
}
|
||||
|
||||
// 条件3: 异常方差变化
|
||||
if (normalized_std > 0.5f) {
|
||||
is_artifact = true;
|
||||
}
|
||||
|
||||
// 条件4: 梯度异常检测
|
||||
if (i > 0) {
|
||||
float gradient = std::abs(current_val - signal[i-1]);
|
||||
float avg_gradient = 0.0f;
|
||||
for (size_t j = 1; j < std::min(ws, i); ++j) {
|
||||
avg_gradient += std::abs(signal[i-j] - signal[i-j-1]);
|
||||
}
|
||||
avg_gradient /= std::min(ws, i);
|
||||
|
||||
if (gradient > 3.0f * avg_gradient) {
|
||||
is_artifact = true;
|
||||
}
|
||||
}
|
||||
|
||||
if (is_artifact) {
|
||||
artifact_mask[i] = true;
|
||||
if (!window_values.empty()) {
|
||||
std::sort(window_values.begin(), window_values.end());
|
||||
result[i] = window_values[window_values.size() / 2]; // 中值
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 3. 频域特征检测
|
||||
std::cout << "开始频域特征检测..." << std::endl;
|
||||
|
||||
// 使用FFT检测异常频率成分
|
||||
const size_t fft_size = std::min(static_cast<size_t>(1024), signal.size());
|
||||
if (fft_size >= 64) {
|
||||
std::vector<float> fft_window(fft_size);
|
||||
std::vector<float> power_spectrum(fft_size / 2);
|
||||
|
||||
for (size_t i = 0; i < signal.size() - fft_size; i += fft_size / 4) { // 50%重叠
|
||||
// 提取窗口数据
|
||||
for (size_t j = 0; j < fft_size; ++j) {
|
||||
fft_window[j] = signal[i + j];
|
||||
}
|
||||
|
||||
// 应用窗函数（Hanning窗）
|
||||
for (size_t j = 0; j < fft_size; ++j) {
|
||||
float window_val = 0.5f * (1.0f - std::cos(2.0f * M_PI * j / (fft_size - 1)));
|
||||
fft_window[j] *= window_val;
|
||||
}
|
||||
|
||||
// 计算功率谱密度
|
||||
for (size_t j = 0; j < fft_size / 2; ++j) {
|
||||
power_spectrum[j] = 0.0f;
|
||||
for (size_t k = 0; k < fft_size; ++k) {
|
||||
float phase = 2.0f * M_PI * j * k / fft_size;
|
||||
float real_part = fft_window[k] * std::cos(phase);
|
||||
float imag_part = fft_window[k] * std::sin(phase);
|
||||
power_spectrum[j] += real_part * real_part + imag_part * imag_part;
|
||||
}
|
||||
power_spectrum[j] /= fft_size;
|
||||
}
|
||||
|
||||
// 检测异常频率成分
|
||||
float total_power = 0.0f, high_freq_power = 0.0f;
|
||||
for (size_t j = 0; j < fft_size / 2; ++j) {
|
||||
total_power += power_spectrum[j];
|
||||
if (j > fft_size / 4) { // 高频成分
|
||||
high_freq_power += power_spectrum[j];
|
||||
}
|
||||
}
|
||||
|
||||
float high_freq_ratio = high_freq_power / (total_power + 1e-6f);
|
||||
|
||||
// 如果高频成分比例异常，标记为运动伪迹
|
||||
if (high_freq_ratio > 0.3f) {
|
||||
for (size_t j = i; j < std::min(i + fft_size, signal.size()); ++j) {
|
||||
artifact_mask[j] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 4. 形态学检测
|
||||
std::cout << "开始形态学检测..." << std::endl;
|
||||
|
||||
// 检测尖峰和突变
|
||||
for (size_t i = 2; i < signal.size() - 2; ++i) {
|
||||
float current = signal[i];
|
||||
float prev = signal[i-1];
|
||||
float next = signal[i+1];
|
||||
float prev2 = signal[i-2];
|
||||
float next2 = signal[i+2];
|
||||
|
||||
// 尖峰检测
|
||||
if ((current > prev && current > next &&
|
||||
current - prev > 2.0f * std::abs(next - prev)) ||
|
||||
(current < prev && current < next &&
|
||||
prev - current > 2.0f * std::abs(next - prev))) {
|
||||
artifact_mask[i] = true;
|
||||
}
|
||||
|
||||
// 突变检测
|
||||
float local_std = std::sqrt((std::pow(prev2 - prev, 2) + std::pow(prev - current, 2) +
|
||||
std::pow(current - next, 2) + std::pow(next - next2, 2)) / 4.0f);
|
||||
if (std::abs(current - prev) > 3.0f * local_std) {
|
||||
artifact_mask[i] = true;
|
||||
}
|
||||
}
|
||||
|
||||
// 5. 智能修复策略
|
||||
std::cout << "开始智能修复..." << std::endl;
|
||||
|
||||
size_t artifact_count = 0;
|
||||
for (size_t i = 0; i < signal.size(); ++i) {
|
||||
if (artifact_mask[i]) {
|
||||
artifact_count++;
|
||||
|
||||
// 根据伪迹类型选择修复策略
|
||||
if (i > 0 && i < signal.size() - 1) {
|
||||
// 策略1: 中值插值
|
||||
std::vector<float> neighbors;
|
||||
for (int offset = -2; offset <= 2; ++offset) {
|
||||
int idx = static_cast<int>(i) + offset;
|
||||
if (idx >= 0 && idx < static_cast<int>(signal.size()) && !artifact_mask[idx]) {
|
||||
neighbors.push_back(signal[idx]);
|
||||
}
|
||||
}
|
||||
|
||||
if (!neighbors.empty()) {
|
||||
if (neighbors.size() >= 3) {
|
||||
// 使用中值
|
||||
std::sort(neighbors.begin(), neighbors.end());
|
||||
result[i] = neighbors[neighbors.size() / 2];
|
||||
} else {
|
||||
// 使用均值
|
||||
float sum = 0.0f;
|
||||
for (float val : neighbors) sum += val;
|
||||
result[i] = sum / neighbors.size();
|
||||
}
|
||||
} else {
|
||||
// 使用线性插值
|
||||
int left_idx = -1, right_idx = -1;
|
||||
for (int offset = 1; offset <= 5; ++offset) {
|
||||
if (i - offset >= 0 && !artifact_mask[i - offset]) {
|
||||
left_idx = i - offset;
|
||||
break;
|
||||
}
|
||||
}
|
||||
for (int offset = 1; offset <= 5; ++offset) {
|
||||
if (i + offset < signal.size() && !artifact_mask[i + offset]) {
|
||||
right_idx = i + offset;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (left_idx >= 0 && right_idx >= 0) {
|
||||
float weight = static_cast<float>(i - left_idx) / (right_idx - left_idx);
|
||||
result[i] = signal[left_idx] * (1.0f - weight) + signal[right_idx] * weight;
|
||||
} else if (left_idx >= 0) {
|
||||
result[i] = signal[left_idx];
|
||||
} else if (right_idx >= 0) {
|
||||
result[i] = signal[right_idx];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 6. 后处理：平滑修复后的信号
|
||||
if (artifact_count > 0) {
|
||||
std::cout << "检测到 " << artifact_count << " 个运动伪迹点，开始后处理..." << std::endl;
|
||||
|
||||
// 对修复区域进行轻微平滑
|
||||
std::vector<float> smoothed = result;
|
||||
const size_t smooth_window = static_cast<size_t>(0.05 * sample_rate); // 50ms窗口
|
||||
|
||||
for (size_t i = smooth_window; i < result.size() - smooth_window; ++i) {
|
||||
if (artifact_mask[i]) {
|
||||
float sum = 0.0f;
|
||||
size_t count = 0;
|
||||
|
||||
for (size_t j = i - smooth_window; j <= i + smooth_window; ++j) {
|
||||
if (j < result.size()) {
|
||||
sum += result[j];
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count > 0) {
|
||||
smoothed[i] = sum / count;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
result = smoothed;
|
||||
}
|
||||
|
||||
// 7. 质量评估
|
||||
float improvement_ratio = 0.0f;
|
||||
if (artifact_count > 0) {
|
||||
// 计算修复前后的信号质量改善
|
||||
float original_variance = 0.0f, repaired_variance = 0.0f;
|
||||
for (size_t i = 1; i < signal.size(); ++i) {
|
||||
original_variance += std::pow(signal[i] - signal[i-1], 2);
|
||||
repaired_variance += std::pow(result[i] - result[i-1], 2);
|
||||
}
|
||||
|
||||
if (original_variance > 0) {
|
||||
improvement_ratio = (original_variance - repaired_variance) / original_variance;
|
||||
}
|
||||
}
|
||||
|
||||
std::cout << "运动伪迹检测和去除完成" << std::endl;
|
||||
std::cout << "检测到的伪迹点数量: " << artifact_count << std::endl;
|
||||
std::cout << "信号质量改善比例: " << (improvement_ratio * 100) << "%" << std::endl;
|
||||
|
||||
return result;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,132 +0,0 @@
|
|||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <fstream>
|
||||
#include <random>
|
||||
|
||||
// 生成包含运动伪迹的测试信号
|
||||
std::vector<float> generate_test_signal(int sample_rate, float duration) {
|
||||
int num_samples = static_cast<int>(sample_rate * duration);
|
||||
std::vector<float> signal(num_samples);
|
||||
|
||||
std::random_device rd;
|
||||
std::mt19937 gen(rd());
|
||||
std::normal_distribution<float> noise(0.0f, 0.05f);
|
||||
|
||||
for (int i = 0; i < num_samples; ++i) {
|
||||
float t = static_cast<float>(i) / sample_rate;
|
||||
|
||||
// 基础信号
|
||||
float base = 0.5 * sin(2 * M_PI * 1.0 * t);
|
||||
|
||||
// 运动伪迹
|
||||
float artifact = 0.0f;
|
||||
if (i % 500 == 0) { // 每500个样本添加伪迹
|
||||
artifact = 0.8f * exp(-std::pow((t - static_cast<int>(t)) * 10, 2));
|
||||
}
|
||||
|
||||
signal[i] = base + artifact + noise(gen);
|
||||
}
|
||||
|
||||
return signal;
|
||||
}
|
||||
|
||||
// 简化的增强版运动伪迹检测
|
||||
std::vector<float> remove_motion_artifacts_simple(const std::vector<float>& signal, double sample_rate) {
|
||||
if (signal.empty()) return signal;
|
||||
|
||||
std::vector<float> result = signal;
|
||||
const size_t window_size = static_cast<size_t>(0.5 * sample_rate);
|
||||
|
||||
if (signal.size() < window_size) return result;
|
||||
|
||||
std::vector<bool> artifact_mask(signal.size(), false);
|
||||
|
||||
// 多尺度检测
|
||||
for (size_t i = window_size; i < signal.size(); ++i) {
|
||||
std::vector<float> window;
|
||||
for (size_t j = i - window_size; j < i; ++j) {
|
||||
window.push_back(signal[j]);
|
||||
}
|
||||
|
||||
float sum = 0.0f, sum_sq = 0.0f;
|
||||
for (float val : window) {
|
||||
sum += val;
|
||||
sum_sq += val * val;
|
||||
}
|
||||
|
||||
float mean = sum / window_size;
|
||||
float variance = (sum_sq / window_size) - (mean * mean);
|
||||
|
||||
if (variance > 1e-6f) {
|
||||
float std_dev = std::sqrt(variance);
|
||||
float current = signal[i];
|
||||
float z_score = std::abs(current - mean) / (std_dev + 1e-6f);
|
||||
|
||||
if (z_score > 4.0f) {
|
||||
artifact_mask[i] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 修复伪迹
|
||||
size_t count = 0;
|
||||
for (size_t i = 0; i < signal.size(); ++i) {
|
||||
if (artifact_mask[i]) {
|
||||
count++;
|
||||
|
||||
// 使用邻域中值
|
||||
std::vector<float> neighbors;
|
||||
for (int offset = -2; offset <= 2; ++offset) {
|
||||
int idx = static_cast<int>(i) + offset;
|
||||
if (idx >= 0 && idx < static_cast<int>(signal.size()) && !artifact_mask[idx]) {
|
||||
neighbors.push_back(signal[idx]);
|
||||
}
|
||||
}
|
||||
|
||||
if (!neighbors.empty()) {
|
||||
std::sort(neighbors.begin(), neighbors.end());
|
||||
result[i] = neighbors[neighbors.size() / 2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::cout << "检测到 " << count << " 个运动伪迹点" << std::endl;
|
||||
return result;
|
||||
}
|
||||
|
||||
// 保存信号到CSV
|
||||
void save_to_csv(const std::vector<float>& signal, const std::string& filename, double sample_rate) {
|
||||
std::ofstream file(filename);
|
||||
if (!file.is_open()) return;
|
||||
|
||||
file << "Time(s),Amplitude\n";
|
||||
for (size_t i = 0; i < signal.size(); ++i) {
|
||||
float time = static_cast<float>(i) / sample_rate;
|
||||
file << time << "," << signal[i] << "\n";
|
||||
}
|
||||
file.close();
|
||||
}
|
||||
|
||||
int main() {
|
||||
std::cout << "=== 增强版运动伪迹检测测试 ===" << std::endl;
|
||||
|
||||
const int sample_rate = 100;
|
||||
const float duration = 5.0f;
|
||||
|
||||
// 生成测试信号
|
||||
std::vector<float> original = generate_test_signal(sample_rate, duration);
|
||||
std::cout << "生成测试信号，长度: " << original.size() << " 样本" << std::endl;
|
||||
|
||||
// 应用运动伪迹检测和去除
|
||||
std::vector<float> cleaned = remove_motion_artifacts_simple(original, sample_rate);
|
||||
|
||||
// 保存结果
|
||||
save_to_csv(original, "original_with_artifacts.csv", sample_rate);
|
||||
save_to_csv(cleaned, "cleaned_signal.csv", sample_rate);
|
||||
|
||||
std::cout << "测试完成，结果已保存到CSV文件" << std::endl;
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,339 +0,0 @@
|
|||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <fstream>
|
||||
#include <random>
|
||||
#include <algorithm>
|
||||
|
||||
// 模拟PPG信号生成函数，包含运动伪迹
|
||||
std::vector<float> generate_test_ppg_with_artifacts(int sample_rate, float duration_seconds) {
|
||||
int num_samples = static_cast<int>(sample_rate * duration_seconds);
|
||||
std::vector<float> signal(num_samples);
|
||||
|
||||
std::random_device rd;
|
||||
std::mt19937 gen(rd());
|
||||
std::normal_distribution<float> noise_dist(0.0f, 0.05f);
|
||||
std::uniform_real_distribution<float> artifact_dist(0.0f, 1.0f);
|
||||
|
||||
for (int i = 0; i < num_samples; ++i) {
|
||||
float t = static_cast<float>(i) / sample_rate;
|
||||
|
||||
// 基础PPG信号（心率约60bpm）
|
||||
float ppg_component = 0.5 * sin(2 * M_PI * 1.0 * t);
|
||||
|
||||
// 呼吸调制
|
||||
float respiration = 0.1 * sin(2 * M_PI * 0.2 * t);
|
||||
|
||||
// 基线漂移
|
||||
float baseline_drift = 0.05 * sin(2 * M_PI * 0.01 * t);
|
||||
|
||||
// 高频噪声
|
||||
float high_freq_noise = 0.02 * sin(2 * M_PI * 10.0 * t) +
|
||||
0.02 * sin(2 * M_PI * 15.0 * t);
|
||||
|
||||
// 运动伪迹（随机出现）
|
||||
float motion_artifact = 0.0f;
|
||||
if (artifact_dist(gen) < 0.01f) { // 1%概率出现运动伪迹
|
||||
float artifact_type = artifact_dist(gen);
|
||||
if (artifact_type < 0.3f) {
|
||||
// 尖峰伪迹
|
||||
motion_artifact = 0.8f * exp(-std::pow((t - static_cast<int>(t)) * 10, 2));
|
||||
} else if (artifact_type < 0.6f) {
|
||||
// 基线跳跃
|
||||
motion_artifact = 0.3f * (artifact_dist(gen) - 0.5f);
|
||||
} else {
|
||||
// 高频振荡
|
||||
motion_artifact = 0.2f * sin(2 * M_PI * 25.0 * t) * exp(-std::pow((t - static_cast<int>(t)) * 5, 2));
|
||||
}
|
||||
}
|
||||
|
||||
// 组合所有成分
|
||||
signal[i] = ppg_component + respiration + baseline_drift + high_freq_noise + motion_artifact + noise_dist(gen);
|
||||
}
|
||||
|
||||
return signal;
|
||||
}
|
||||
|
||||
// 简化的增强版运动伪迹检测算法（用于测试）
|
||||
std::vector<float> enhanced_motion_artifact_removal(const std::vector<float>& signal, double sample_rate) {
|
||||
if (signal.empty()) return signal;
|
||||
|
||||
std::vector<float> result = signal;
|
||||
const size_t min_window_size = static_cast<size_t>(0.1 * sample_rate);
|
||||
|
||||
if (signal.size() < min_window_size) return result;
|
||||
|
||||
// 多尺度检测窗口
|
||||
std::vector<size_t> window_sizes = {
|
||||
static_cast<size_t>(0.1 * sample_rate),
|
||||
static_cast<size_t>(0.5 * sample_rate),
|
||||
static_cast<size_t>(1.0 * sample_rate)
|
||||
};
|
||||
|
||||
std::vector<bool> artifact_mask(signal.size(), false);
|
||||
|
||||
// 统计特征检测
|
||||
for (size_t ws : window_sizes) {
|
||||
if (signal.size() < ws) continue;
|
||||
|
||||
for (size_t i = ws; i < signal.size(); ++i) {
|
||||
std::vector<float> window_data;
|
||||
for (size_t j = i - ws; j < i; ++j) {
|
||||
window_data.push_back(signal[j]);
|
||||
}
|
||||
|
||||
// 计算统计特征
|
||||
float sum = 0.0f, sum_sq = 0.0f;
|
||||
for (float val : window_data) {
|
||||
sum += val;
|
||||
sum_sq += val * val;
|
||||
}
|
||||
|
||||
float mean = sum / ws;
|
||||
float variance = (sum_sq / ws) - (mean * mean);
|
||||
|
||||
if (variance > 1e-6f) {
|
||||
float std_dev = std::sqrt(variance);
|
||||
float current_val = signal[i];
|
||||
float z_score = std::abs(current_val - mean) / (std_dev + 1e-6f);
|
||||
|
||||
// 检测异常值
|
||||
if (z_score > 4.0f) {
|
||||
artifact_mask[i] = true;
|
||||
}
|
||||
|
||||
// 梯度异常检测
|
||||
if (i > 0) {
|
||||
float gradient = std::abs(current_val - signal[i-1]);
|
||||
float avg_gradient = 0.0f;
|
||||
for (size_t j = 1; j < std::min(ws, i); ++j) {
|
||||
avg_gradient += std::abs(signal[i-j] - signal[i-j-1]);
|
||||
}
|
||||
avg_gradient /= std::min(ws, i);
|
||||
|
||||
if (gradient > 3.0f * avg_gradient) {
|
||||
artifact_mask[i] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 形态学检测
|
||||
for (size_t i = 2; i < signal.size() - 2; ++i) {
|
||||
float current = signal[i];
|
||||
float prev = signal[i-1];
|
||||
float next = signal[i+1];
|
||||
|
||||
// 尖峰检测
|
||||
if ((current > prev && current > next &&
|
||||
current - prev > 2.0f * std::abs(next - prev)) ||
|
||||
(current < prev && current < next &&
|
||||
prev - current > 2.0f * std::abs(next - prev))) {
|
||||
artifact_mask[i] = true;
|
||||
}
|
||||
}
|
||||
|
||||
// 智能修复
|
||||
size_t artifact_count = 0;
|
||||
for (size_t i = 0; i < signal.size(); ++i) {
|
||||
if (artifact_mask[i]) {
|
||||
artifact_count++;
|
||||
|
||||
// 使用中值插值
|
||||
std::vector<float> neighbors;
|
||||
for (int offset = -2; offset <= 2; ++offset) {
|
||||
int idx = static_cast<int>(i) + offset;
|
||||
if (idx >= 0 && idx < static_cast<int>(signal.size()) && !artifact_mask[idx]) {
|
||||
neighbors.push_back(signal[idx]);
|
||||
}
|
||||
}
|
||||
|
||||
if (!neighbors.empty()) {
|
||||
if (neighbors.size() >= 3) {
|
||||
std::sort(neighbors.begin(), neighbors.end());
|
||||
result[i] = neighbors[neighbors.size() / 2];
|
||||
} else {
|
||||
float sum = 0.0f;
|
||||
for (float val : neighbors) sum += val;
|
||||
result[i] = sum / neighbors.size();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::cout << "检测到 " << artifact_count << " 个运动伪迹点" << std::endl;
|
||||
return result;
|
||||
}
|
||||
|
||||
// 计算信号质量指标
|
||||
struct SignalQualityMetrics {
|
||||
float snr_db;
|
||||
float variance;
|
||||
float peak_to_peak;
|
||||
float artifact_ratio;
|
||||
};
|
||||
|
||||
SignalQualityMetrics calculate_signal_quality(const std::vector<float>& signal, double sample_rate) {
|
||||
SignalQualityMetrics metrics;
|
||||
|
||||
if (signal.empty()) {
|
||||
metrics.snr_db = 0.0f;
|
||||
metrics.variance = 0.0f;
|
||||
metrics.peak_to_peak = 0.0f;
|
||||
metrics.artifact_ratio = 0.0f;
|
||||
return metrics;
|
||||
}
|
||||
|
||||
// 计算基本统计量
|
||||
float sum = 0.0f, sum_sq = 0.0f;
|
||||
float min_val = signal[0], max_val = signal[0];
|
||||
|
||||
for (float val : signal) {
|
||||
sum += val;
|
||||
sum_sq += val * val;
|
||||
min_val = std::min(min_val, val);
|
||||
max_val = std::max(max_val, val);
|
||||
}
|
||||
|
||||
float mean = sum / signal.size();
|
||||
metrics.variance = (sum_sq / signal.size()) - (mean * mean);
|
||||
metrics.peak_to_peak = max_val - min_val;
|
||||
|
||||
// 计算SNR
|
||||
float signal_power = 0.0f;
|
||||
float noise_power = 0.0f;
|
||||
|
||||
for (float val : signal) {
|
||||
signal_power += std::pow(val - mean, 2);
|
||||
}
|
||||
signal_power /= signal.size();
|
||||
|
||||
for (size_t i = 1; i < signal.size(); ++i) {
|
||||
float diff = signal[i] - signal[i-1];
|
||||
noise_power += diff * diff;
|
||||
}
|
||||
noise_power /= (signal.size() - 1);
|
||||
|
||||
if (noise_power > 1e-6f) {
|
||||
metrics.snr_db = 10.0f * std::log10(signal_power / noise_power);
|
||||
} else {
|
||||
metrics.snr_db = 0.0f;
|
||||
}
|
||||
|
||||
// 估算伪迹比例（通过异常梯度检测）
|
||||
size_t artifact_count = 0;
|
||||
for (size_t i = 1; i < signal.size(); ++i) {
|
||||
float gradient = std::abs(signal[i] - signal[i-1]);
|
||||
if (gradient > 3.0f * std::sqrt(metrics.variance)) {
|
||||
artifact_count++;
|
||||
}
|
||||
}
|
||||
|
||||
metrics.artifact_ratio = static_cast<float>(artifact_count) / signal.size();
|
||||
|
||||
return metrics;
|
||||
}
|
||||
|
||||
// 保存信号到CSV文件
|
||||
void save_signal_to_csv(const std::vector<float>& signal, const std::string& filename, double sample_rate) {
|
||||
std::ofstream file(filename);
|
||||
if (!file.is_open()) {
|
||||
std::cerr << "无法创建文件: " << filename << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
file << "Time(s),Amplitude\n";
|
||||
for (size_t i = 0; i < signal.size(); ++i) {
|
||||
float time = static_cast<float>(i) / sample_rate;
|
||||
file << time << "," << signal[i] << "\n";
|
||||
}
|
||||
file.close();
|
||||
std::cout << "信号已保存到: " << filename << std::endl;
|
||||
}
|
||||
|
||||
int main() {
|
||||
std::cout << "=== 增强版运动伪迹检测和去除算法测试 ===" << std::endl;
|
||||
|
||||
// 测试参数
|
||||
const int sample_rate = 100; // 100Hz采样率
|
||||
const float duration = 10.0f; // 10秒数据
|
||||
|
||||
std::cout << "采样率: " << sample_rate << "Hz" << std::endl;
|
||||
std::cout << "信号长度: " << duration << "秒" << std::endl;
|
||||
|
||||
// 1. 生成包含运动伪迹的测试信号
|
||||
std::cout << "\n1. 生成测试信号..." << std::endl;
|
||||
std::vector<float> original_signal = generate_test_ppg_with_artifacts(sample_rate, duration);
|
||||
std::cout << "测试信号生成完成，长度: " << original_signal.size() << " 样本" << std::endl;
|
||||
|
||||
// 2. 计算原始信号质量
|
||||
std::cout << "\n2. 分析原始信号质量..." << std::endl;
|
||||
SignalQualityMetrics original_metrics = calculate_signal_quality(original_signal, sample_rate);
|
||||
std::cout << "原始信号质量指标:" << std::endl;
|
||||
std::cout << " SNR: " << original_metrics.snr_db << " dB" << std::endl;
|
||||
std::cout << " 方差: " << original_metrics.variance << std::endl;
|
||||
std::cout << " 峰峰值: " << original_metrics.peak_to_peak << std::endl;
|
||||
std::cout << " 伪迹比例: " << (original_metrics.artifact_ratio * 100) << "%" << std::endl;
|
||||
|
||||
// 3. 应用运动伪迹检测和去除
|
||||
std::cout << "\n3. 应用增强版运动伪迹检测和去除..." << std::endl;
|
||||
std::vector<float> cleaned_signal = enhanced_motion_artifact_removal(original_signal, sample_rate);
|
||||
|
||||
// 4. 计算处理后信号质量
|
||||
std::cout << "\n4. 分析处理后信号质量..." << std::endl;
|
||||
SignalQualityMetrics cleaned_metrics = calculate_signal_quality(cleaned_signal, sample_rate);
|
||||
std::cout << "处理后信号质量指标:" << std::endl;
|
||||
std::cout << " SNR: " << cleaned_metrics.snr_db << " dB" << std::endl;
|
||||
std::cout << " 方差: " << cleaned_metrics.variance << std::endl;
|
||||
std::cout << " 峰峰值: " << cleaned_metrics.peak_to_peak << std::endl;
|
||||
std::cout << " 伪迹比例: " << (cleaned_metrics.artifact_ratio * 100) << "%" << std::endl;
|
||||
|
||||
// 5. 计算改善效果
|
||||
std::cout << "\n5. 改善效果分析..." << std::endl;
|
||||
float snr_improvement = cleaned_metrics.snr_db - original_metrics.snr_db;
|
||||
float variance_reduction = (original_metrics.variance - cleaned_metrics.variance) / original_metrics.variance * 100;
|
||||
float artifact_reduction = (original_metrics.artifact_ratio - cleaned_metrics.artifact_ratio) / original_metrics.artifact_ratio * 100;
|
||||
|
||||
std::cout << "改善效果:" << std::endl;
|
||||
std::cout << " SNR改善: " << snr_improvement << " dB" << std::endl;
|
||||
std::cout << " 方差减少: " << variance_reduction << "%" << std::endl;
|
||||
std::cout << " 伪迹减少: " << artifact_reduction << "%" << std::endl;
|
||||
|
||||
// 6. 保存结果
|
||||
std::cout << "\n6. 保存结果..." << std::endl;
|
||||
save_signal_to_csv(original_signal, "original_signal_with_artifacts.csv", sample_rate);
|
||||
save_signal_to_csv(cleaned_signal, "cleaned_signal.csv", sample_rate);
|
||||
|
||||
// 7. 生成对比报告
|
||||
std::ofstream report("motion_artifact_removal_report.txt");
|
||||
if (report.is_open()) {
|
||||
report << "=== 运动伪迹检测和去除算法测试报告 ===" << std::endl;
|
||||
report << "测试时间: " << std::time(nullptr) << std::endl;
|
||||
report << "采样率: " << sample_rate << "Hz" << std::endl;
|
||||
report << "信号长度: " << duration << "秒" << std::endl;
|
||||
report << "\n原始信号质量:" << std::endl;
|
||||
report << " SNR: " << original_metrics.snr_db << " dB" << std::endl;
|
||||
report << " 方差: " << original_metrics.variance << std::endl;
|
||||
report << " 峰峰值: " << original_metrics.peak_to_peak << std::endl;
|
||||
report << " 伪迹比例: " << (original_metrics.artifact_ratio * 100) << "%" << std::endl;
|
||||
report << "\n处理后信号质量:" << std::endl;
|
||||
report << " SNR: " << cleaned_metrics.snr_db << " dB" << std::endl;
|
||||
report << " 方差: " << cleaned_metrics.variance << std::endl;
|
||||
report << " 峰峰值: " << cleaned_metrics.peak_to_peak << std::endl;
|
||||
report << " 伪迹比例: " << (cleaned_metrics.artifact_ratio * 100) << "%" << std::endl;
|
||||
report << "\n改善效果:" << std::endl;
|
||||
report << " SNR改善: " << snr_improvement << " dB" << std::endl;
|
||||
report << " 方差减少: " << variance_reduction << "%" << std::endl;
|
||||
report << " 伪迹减少: " << artifact_reduction << "%" << std::endl;
|
||||
report.close();
|
||||
std::cout << "测试报告已保存到: motion_artifact_removal_report.txt" << std::endl;
|
||||
}
|
||||
|
||||
std::cout << "\n=== 测试完成 ===" << std::endl;
|
||||
std::cout << "请检查生成的CSV文件和测试报告" << std::endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,123 +0,0 @@
|
|||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <fstream>
|
||||
|
||||
// 模拟PPG信号生成函数
|
||||
std::vector<float> generate_test_ppg_signal(int sample_rate, float duration_seconds) {
|
||||
int num_samples = static_cast<int>(sample_rate * duration_seconds);
|
||||
std::vector<float> signal(num_samples);
|
||||
|
||||
for (int i = 0; i < num_samples; ++i) {
|
||||
float t = static_cast<float>(i) / sample_rate;
|
||||
|
||||
// 生成包含多个频率成分的测试信号
|
||||
float ppg_component = 0.5 * sin(2 * M_PI * 1.2 * t); // 1.2Hz PPG信号
|
||||
float noise_component = 0.1 * sin(2 * M_PI * 50.0 * t); // 50Hz工频干扰
|
||||
float high_freq_noise = 0.05 * sin(2 * M_PI * 100.0 * t); // 100Hz高频噪声
|
||||
|
||||
signal[i] = ppg_component + noise_component + high_freq_noise;
|
||||
}
|
||||
|
||||
return signal;
|
||||
}
|
||||
|
||||
// 简化的陷波滤波器实现（用于测试）
|
||||
std::vector<float> simple_notch_filter(const std::vector<float>& input, double sample_rate,
|
||||
double target_freq, double bandwidth) {
|
||||
if (input.empty()) return {};
|
||||
|
||||
// 检查采样率和目标频率的合理性
|
||||
if (sample_rate <= 0 || target_freq <= 0) {
|
||||
std::cout << "警告: 采样率或目标频率无效，返回原始信号" << std::endl;
|
||||
return input;
|
||||
}
|
||||
|
||||
// 检查目标频率是否接近奈奎斯特频率
|
||||
const double nyquist_freq = sample_rate / 2.0;
|
||||
if (target_freq >= nyquist_freq * 0.8) {
|
||||
std::cout << "警告: 目标频率 " << target_freq << "Hz 接近奈奎斯特频率 " << nyquist_freq << "Hz，跳过陷波滤波" << std::endl;
|
||||
return input;
|
||||
}
|
||||
|
||||
// 使用带阻滤波替代陷波滤波
|
||||
std::vector<float> output = input;
|
||||
const double low_cutoff = target_freq - bandwidth/2;
|
||||
const double high_cutoff = target_freq + bandwidth/2;
|
||||
|
||||
// 简单的带阻滤波实现
|
||||
for (size_t i = 2; i < output.size(); ++i) {
|
||||
// 简单的移动平均滤波
|
||||
output[i] = (input[i] + input[i-1] + input[i-2]) / 3.0f;
|
||||
}
|
||||
|
||||
return output;
|
||||
}
|
||||
|
||||
// 保存信号到CSV文件
|
||||
void save_signal_to_csv(const std::vector<float>& signal, const std::string& filename, double sample_rate) {
|
||||
std::ofstream file(filename);
|
||||
if (!file.is_open()) {
|
||||
std::cerr << "无法创建文件: " << filename << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
file << "Time(s),Amplitude\n";
|
||||
for (size_t i = 0; i < signal.size(); ++i) {
|
||||
float time = static_cast<float>(i) / sample_rate;
|
||||
file << time << "," << signal[i] << "\n";
|
||||
}
|
||||
file.close();
|
||||
std::cout << "信号已保存到: " << filename << std::endl;
|
||||
}
|
||||
|
||||
int main() {
|
||||
std::cout << "=== 50Hz陷波滤波器测试程序 ===" << std::endl;
|
||||
|
||||
// 测试不同的采样率
|
||||
std::vector<int> test_sample_rates = {50, 100, 250, 500};
|
||||
|
||||
for (int sample_rate : test_sample_rates) {
|
||||
std::cout << "\n--- 测试采样率: " << sample_rate << "Hz ---" << std::endl;
|
||||
|
||||
// 生成测试信号
|
||||
std::vector<float> test_signal = generate_test_ppg_signal(sample_rate, 2.0); // 2秒数据
|
||||
std::cout << "生成测试信号，长度: " << test_signal.size() << " 样本" << std::endl;
|
||||
|
||||
// 保存原始信号
|
||||
std::string original_filename = "original_signal_" + std::to_string(sample_rate) + "Hz.csv";
|
||||
save_signal_to_csv(test_signal, original_filename, sample_rate);
|
||||
|
||||
// 测试陷波滤波
|
||||
std::vector<float> filtered_signal;
|
||||
if (sample_rate > 100) {
|
||||
// 高采样率使用标准陷波滤波
|
||||
std::cout << "使用标准陷波滤波 (49.5-50.5Hz)" << std::endl;
|
||||
filtered_signal = simple_notch_filter(test_signal, sample_rate, 50.0, 1.0);
|
||||
} else {
|
||||
// 低采样率使用带阻滤波
|
||||
std::cout << "使用带阻滤波 (45-55Hz)，避免过度衰减" << std::endl;
|
||||
filtered_signal = simple_notch_filter(test_signal, sample_rate, 50.0, 10.0);
|
||||
}
|
||||
|
||||
// 保存滤波后的信号
|
||||
std::string filtered_filename = "filtered_signal_" + std::to_string(sample_rate) + "Hz.csv";
|
||||
save_signal_to_csv(filtered_signal, filtered_filename, sample_rate);
|
||||
|
||||
// 计算信号统计信息
|
||||
float original_rms = 0, filtered_rms = 0;
|
||||
for (float val : test_signal) original_rms += val * val;
|
||||
for (float val : filtered_signal) filtered_rms += val * val;
|
||||
original_rms = sqrt(original_rms / test_signal.size());
|
||||
filtered_rms = sqrt(filtered_rms / filtered_signal.size());
|
||||
|
||||
std::cout << "原始信号RMS: " << original_rms << std::endl;
|
||||
std::cout << "滤波后信号RMS: " << filtered_rms << std::endl;
|
||||
std::cout << "信号保留率: " << (filtered_rms / original_rms * 100) << "%" << std::endl;
|
||||
}
|
||||
|
||||
std::cout << "\n=== 测试完成 ===" << std::endl;
|
||||
std::cout << "请检查生成的CSV文件，验证滤波效果" << std::endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue