diff --git a/ENHANCED_MOTION_ARTIFACT_DETECTION_README.md b/ENHANCED_MOTION_ARTIFACT_DETECTION_README.md index 19530aa..dfa5e87 100644 --- a/ENHANCED_MOTION_ARTIFACT_DETECTION_README.md +++ b/ENHANCED_MOTION_ARTIFACT_DETECTION_README.md @@ -221,3 +221,4 @@ g++ -o test_motion_artifact test_motion_artifact.cpp - 质量改善：SNR提升8-15dB + diff --git a/main.cpp b/main.cpp index 22d1693..5603ff5 100644 --- a/main.cpp +++ b/main.cpp @@ -2,7 +2,7 @@ #include #include #include - +#include // 引入 Windows API 头文件 std::vector 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; diff --git a/src/signal_processor/signal_processor.cpp b/src/signal_processor/signal_processor.cpp index aa04cb3..5ed623f 100644 --- a/src/signal_processor/signal_processor.cpp +++ b/src/signal_processor/signal_processor.cpp @@ -1067,300 +1067,44 @@ float SignalProcessor::calculate_PPG_sqi(const std::vector& red_channel, std::vector SignalProcessor::remove_motion_artifacts(const std::vector& signal, double sample_rate) { if (signal.empty()) return signal; - std::cout << "开始增强版运动伪迹检测和去除，信号长度: " << signal.size() << " 样本" << std::endl; - std::vector result = signal; - const size_t min_window_size = static_cast(0.1 * sample_rate); // 最小窗口0.1秒 - const size_t max_window_size = static_cast(2.0 * sample_rate); // 最大窗口2秒 + const size_t window_size = static_cast(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; - // 1. 多尺度运动伪迹检测 - std::vector artifact_mask(signal.size(), false); - - // 使用多个窗口大小进行检测 - std::vector window_sizes = { - static_cast(0.1 * sample_rate), // 0.1秒 - 快速运动 - static_cast(0.5 * sample_rate), // 0.5秒 - 中等运动 - static_cast(1.0 * sample_rate), // 1.0秒 - 慢速运动 - static_cast(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 = 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]; + } - for (size_t i = ws; i < signal.size(); ++i) { - // 计算滑动窗口统计特征 - std::vector window_data; - for (size_t j = i - ws; j < i; ++j) { - window_data.push_back(signal[j]); + float mean = sum / window_size; + float variance = (sum_sq / window_size) - (mean * mean); + float std_dev = std::sqrt(std::max(0.0f, variance)); + + // 检测异常值（可能的运动伪迹） + float current_val = signal[i]; + float z_score = std::abs(current_val - mean) / (std_dev + 1e-6f); + + // 如果Z-score > 3，认为是运动伪迹 + if (z_score > 3.0f) { + // 使用中值滤波替换异常值 + std::vector window_values; + for (size_t j = std::max(static_cast(0), i - window_size/2); + j < std::min(signal.size(), i + window_size/2); ++j) { + window_values.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; - } - - // 条件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(1024), signal.size()); - if (fft_size >= 64) { - std::vector fft_window(fft_size); - std::vector 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 neighbors; - for (int offset = -2; offset <= 2; ++offset) { - int idx = static_cast(i) + offset; - if (idx >= 0 && idx < static_cast(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(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 smoothed = result; - const size_t smooth_window = static_cast(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; } diff --git a/test_motion_artifact.cpp b/test_motion_artifact.cpp deleted file mode 100644 index ba84e45..0000000 --- a/test_motion_artifact.cpp +++ /dev/null @@ -1,132 +0,0 @@ -#include -#include -#include -#include -#include - -// 生成包含运动伪迹的测试信号 -std::vector generate_test_signal(int sample_rate, float duration) { - int num_samples = static_cast(sample_rate * duration); - std::vector signal(num_samples); - - std::random_device rd; - std::mt19937 gen(rd()); - std::normal_distribution noise(0.0f, 0.05f); - - for (int i = 0; i < num_samples; ++i) { - float t = static_cast(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(t)) * 10, 2)); - } - - signal[i] = base + artifact + noise(gen); - } - - return signal; -} - -// 简化的增强版运动伪迹检测 -std::vector remove_motion_artifacts_simple(const std::vector& signal, double sample_rate) { - if (signal.empty()) return signal; - - std::vector result = signal; - const size_t window_size = static_cast(0.5 * sample_rate); - - if (signal.size() < window_size) return result; - - std::vector artifact_mask(signal.size(), false); - - // 多尺度检测 - for (size_t i = window_size; i < signal.size(); ++i) { - std::vector 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 neighbors; - for (int offset = -2; offset <= 2; ++offset) { - int idx = static_cast(i) + offset; - if (idx >= 0 && idx < static_cast(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& 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(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 original = generate_test_signal(sample_rate, duration); - std::cout << "生成测试信号，长度: " << original.size() << " 样本" << std::endl; - - // 应用运动伪迹检测和去除 - std::vector 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; -} - - diff --git a/test_motion_artifact_removal.cpp b/test_motion_artifact_removal.cpp deleted file mode 100644 index 85faae6..0000000 --- a/test_motion_artifact_removal.cpp +++ /dev/null @@ -1,339 +0,0 @@ -#include -#include -#include -#include -#include -#include - -// 模拟PPG信号生成函数，包含运动伪迹 -std::vector generate_test_ppg_with_artifacts(int sample_rate, float duration_seconds) { - int num_samples = static_cast(sample_rate * duration_seconds); - std::vector signal(num_samples); - - std::random_device rd; - std::mt19937 gen(rd()); - std::normal_distribution noise_dist(0.0f, 0.05f); - std::uniform_real_distribution artifact_dist(0.0f, 1.0f); - - for (int i = 0; i < num_samples; ++i) { - float t = static_cast(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(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(t)) * 5, 2)); - } - } - - // 组合所有成分 - signal[i] = ppg_component + respiration + baseline_drift + high_freq_noise + motion_artifact + noise_dist(gen); - } - - return signal; -} - -// 简化的增强版运动伪迹检测算法（用于测试） -std::vector enhanced_motion_artifact_removal(const std::vector& signal, double sample_rate) { - if (signal.empty()) return signal; - - std::vector result = signal; - const size_t min_window_size = static_cast(0.1 * sample_rate); - - if (signal.size() < min_window_size) return result; - - // 多尺度检测窗口 - std::vector window_sizes = { - static_cast(0.1 * sample_rate), - static_cast(0.5 * sample_rate), - static_cast(1.0 * sample_rate) - }; - - std::vector 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 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 neighbors; - for (int offset = -2; offset <= 2; ++offset) { - int idx = static_cast(i) + offset; - if (idx >= 0 && idx < static_cast(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& 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(artifact_count) / signal.size(); - - return metrics; -} - -// 保存信号到CSV文件 -void save_signal_to_csv(const std::vector& 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(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 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 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; -} - - diff --git a/test_notch_filter.cpp b/test_notch_filter.cpp deleted file mode 100644 index 3dcff93..0000000 --- a/test_notch_filter.cpp +++ /dev/null @@ -1,123 +0,0 @@ -#include -#include -#include -#include - -// 模拟PPG信号生成函数 -std::vector generate_test_ppg_signal(int sample_rate, float duration_seconds) { - int num_samples = static_cast(sample_rate * duration_seconds); - std::vector signal(num_samples); - - for (int i = 0; i < num_samples; ++i) { - float t = static_cast(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 simple_notch_filter(const std::vector& 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 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& 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(i) / sample_rate; - file << time << "," << signal[i] << "\n"; - } - file.close(); - std::cout << "信号已保存到: " << filename << std::endl; -} - -int main() { - std::cout << "=== 50Hz陷波滤波器测试程序 ===" << std::endl; - - // 测试不同的采样率 - std::vector test_sample_rates = {50, 100, 250, 500}; - - for (int sample_rate : test_sample_rates) { - std::cout << "\n--- 测试采样率: " << sample_rate << "Hz ---" << std::endl; - - // 生成测试信号 - std::vector 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 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; -}