diff --git a/.kotlin/sessions/kotlin-compiler-15274677907984621681.salive b/.kotlin/sessions/kotlin-compiler-15274677907984621681.salive deleted file mode 100644 index e69de29..0000000 diff --git a/app/src/main/cpp/include/cpp/indicator_cal.h b/app/src/main/cpp/include/cpp/indicator_cal.h index 6107f3a..0393ec1 100644 --- a/app/src/main/cpp/include/cpp/indicator_cal.h +++ b/app/src/main/cpp/include/cpp/indicator_cal.h @@ -42,6 +42,10 @@ public: float calculate_signal_quality(const std::vector& signal); float estimate_heart_rate_from_frequency(const std::vector& signal, float sample_rate); + // 频域分析相关函数 + float calculate_snr_frequency_domain(const std::vector& signal); + float calculate_snr_time_domain(const std::vector& signal); + // 综合指标计算 std::map calculate_all_ecg_metrics(const SensorData& ecg_data, float sample_rate); std::map calculate_all_ppg_metrics(const SensorData& ppg_data, float sample_rate); diff --git a/app/src/main/cpp/include/cpp/signal_processor.h b/app/src/main/cpp/include/cpp/signal_processor.h index 7580a5d..9fd1e94 100644 --- a/app/src/main/cpp/include/cpp/signal_processor.h +++ b/app/src/main/cpp/include/cpp/signal_processor.h @@ -10,6 +10,16 @@ enum class filtertype { lowpass,highpass,notchpass,bandpass,bandstop }; + +// 一阶滤波器参数结构 +struct FirstOrderFilterParams { + double b0, b1; // 分子系数 + double a1; // 分母系数 (a0=1) + + FirstOrderFilterParams() : b0(1.0), b1(0.0), a1(0.0) {} + FirstOrderFilterParams(double b0_, double b1_, double a1_) + : b0(b0_), b1(b1_), a1(a1_) {} +}; // 特征数据结构 struct ECGChannelFeatures { std::vector r_peaks; // R波位置 @@ -82,6 +92,13 @@ public: // 设置滤波器参数 void setFilterParams(const std::string& filter_type, const FilterConfig& params); + // 一阶滤波器相关函数 + FirstOrderFilterParams calculateFirstOrderFilterParams(double sample_rate, double cutoff_freq, filtertype type); + std::vector applyFirstOrderFilter(const std::vector& input, const FirstOrderFilterParams& params); + + // 一阶带通滤波器辅助函数 + std::pair calculateFirstOrderBandpassParams(double sample_rate, double low_cutoff, double high_cutoff); + std::vector applyFirstOrderBandpassFilter(const std::vector& input, const FirstOrderFilterParams& highpass_params, const FirstOrderFilterParams& lowpass_params); // 修改为public访问权限 public: diff --git a/app/src/main/cpp/src/data_praser.cpp b/app/src/main/cpp/src/data_praser.cpp index 4aa4baf..bfe5dd4 100644 --- a/app/src/main/cpp/src/data_praser.cpp +++ b/app/src/main/cpp/src/data_praser.cpp @@ -311,7 +311,7 @@ SensorData parse_snore(const uint8_t* data) { } // 赋值原始二进制数据 - result.raw_data.assign(data, data + 6 + 232); + result.raw_data.assign(data, data + 238); return result; } diff --git a/app/src/main/cpp/src/indicator_cal.cpp b/app/src/main/cpp/src/indicator_cal.cpp index a4de769..2d1ed0d 100644 --- a/app/src/main/cpp/src/indicator_cal.cpp +++ b/app/src/main/cpp/src/indicator_cal.cpp @@ -979,38 +979,37 @@ std::vector MetricsCalculator::detect_pulse_peaks(const std::vector& signal) { if (signal.empty()) return 0.0f; const size_t n = signal.size(); - // 1. 计算信噪比 (SNR) - float signal_power = 0.0f; - float noise_power = 0.0f; - - // 使用简化的功率谱密度计算(避免复杂的FFT) - // 假设低频部分是信号,高频部分是噪声 - size_t mid_point = n / 2; - for (size_t i = 0; i < mid_point; i++) { - signal_power += signal[i] * signal[i]; - } - for (size_t i = mid_point; i < n; i++) { - noise_power += signal[i] * signal[i]; - } - - float snr = 0.0f; - if (noise_power > 0.0001f) { - snr = 10.0f * std::log10(signal_power / noise_power); - } + // 1. 计算信噪比 (SNR) - 使用频域分析 + float snr = calculate_snr_frequency_domain(signal); // 2. 计算信号连续性 float continuity_score = 0.0f; int discontinuity_count = 0; + // 计算信号的标准差作为自适应阈值的基础 + float mean = 0.0f; + for (float val : signal) { + mean += val; + } + mean /= n; + + float variance = 0.0f; + for (float val : signal) { + float diff = val - mean; + variance += diff * diff; + } + variance /= n; + float std_dev = std::sqrt(variance); + for (size_t i = 1; i < n; i++) { float diff = std::abs(signal[i] - signal[i-1]); - if (diff > 10.0f * std::sqrt(noise_power / n)) { // 自适应阈值 + if (diff > 3.0f * std_dev) { // 使用3倍标准差作为阈值 discontinuity_count++; } } @@ -1199,4 +1198,174 @@ std::map MetricsCalculator::calculate_all_hrv_metrics(const metrics["total_power"] = total_power; return metrics; +} + +// ============================================================================ +// 频域分析相关函数 + +// 简单的FFT实现(使用Cooley-Tukey算法) +std::vector> fft(const std::vector& signal) { + const size_t n = signal.size(); + + // 确保信号长度是2的幂次 + size_t padded_size = 1; + while (padded_size < n) { + padded_size <<= 1; + } + + std::vector> padded_signal(padded_size, 0.0f); + for (size_t i = 0; i < n; i++) { + padded_signal[i] = std::complex(signal[i], 0.0f); + } + + // Cooley-Tukey FFT + std::vector> result = padded_signal; + + // 位反转 + for (size_t i = 1, j = 0; i < padded_size; i++) { + size_t bit = padded_size >> 1; + while (j & bit) { + j ^= bit; + bit >>= 1; + } + j ^= bit; + if (i < j) { + std::swap(result[i], result[j]); + } + } + + // FFT计算 + for (size_t len = 2; len <= padded_size; len <<= 1) { + float angle = -2.0f * M_PI / len; + std::complex wlen(std::cos(angle), std::sin(angle)); + + for (size_t i = 0; i < padded_size; i += len) { + std::complex w(1.0f, 0.0f); + for (size_t j = 0; j < len / 2; j++) { + std::complex u = result[i + j]; + std::complex v = result[i + j + len / 2] * w; + result[i + j] = u + v; + result[i + j + len / 2] = u - v; + w *= wlen; + } + } + } + + return result; +} + +// 计算功率谱密度 +std::vector calculate_power_spectrum(const std::vector>& fft_result) { + std::vector power_spectrum(fft_result.size() / 2); + + for (size_t i = 0; i < power_spectrum.size(); i++) { + float magnitude = std::abs(fft_result[i]); + power_spectrum[i] = magnitude * magnitude; + } + + return power_spectrum; +} + +// 计算指定频带的功率 +float calculate_power_in_frequency_band(const std::vector& power_spectrum, + float sample_rate, + float low_freq, + float high_freq) { + const size_t n = power_spectrum.size(); + const float freq_resolution = sample_rate / (2.0f * n); + + size_t low_bin = static_cast(low_freq / freq_resolution); + size_t high_bin = static_cast(high_freq / freq_resolution); + + // 确保索引在有效范围内 + low_bin = std::min(low_bin, n - 1); + high_bin = std::min(high_bin, n - 1); + + if (low_bin >= high_bin) return 0.0f; + + float power = 0.0f; + for (size_t i = low_bin; i <= high_bin; i++) { + power += power_spectrum[i]; + } + + return power; +} + +// 基于频域分析计算信噪比 +float MetricsCalculator::calculate_snr_frequency_domain(const std::vector& signal) { + if (signal.size() < 64) { + // 信号太短,使用简化的时域方法 + return calculate_snr_time_domain(signal); + } + + // 假设采样率为250Hz(可以根据实际情况调整) + const float sample_rate = 250.0f; + + // 进行FFT变换 + auto fft_result = fft(signal); + + // 计算功率谱 + auto power_spectrum = calculate_power_spectrum(fft_result); + + // 定义生物电信号的频带(Hz) + const float signal_low = 0.5f; // 生物电信号低频下限 + const float signal_high = 30.0f; // 生物电信号高频上限 + const float noise_low = 50.0f; // 噪声频带下限 + const float noise_high = 125.0f; // 噪声频带上限(奈奎斯特频率的一半) + + // 计算信号频带功率 + float signal_power = calculate_power_in_frequency_band(power_spectrum, sample_rate, + signal_low, signal_high); + + // 计算噪声频带功率 + float noise_power = calculate_power_in_frequency_band(power_spectrum, sample_rate, + noise_low, noise_high); + + // 如果噪声功率太小,使用总功率的10%作为噪声估计 + if (noise_power < signal_power * 0.01f) { + float total_power = 0.0f; + for (float p : power_spectrum) { + total_power += p; + } + noise_power = total_power * 0.1f; + } + + // 计算SNR (dB) + float snr = 0.0f; + if (noise_power > 0.0001f) { + snr = 10.0f * std::log10(signal_power / noise_power); + } + + return snr; +} + +// 简化的时域SNR计算(用于短信号) +float MetricsCalculator::calculate_snr_time_domain(const std::vector& signal) { + if (signal.empty()) return 0.0f; + + // 计算信号均值 + float mean = 0.0f; + for (float val : signal) { + mean += val; + } + mean /= signal.size(); + + // 计算信号功率(去直流分量) + float signal_power = 0.0f; + for (float val : signal) { + float centered = val - mean; + signal_power += centered * centered; + } + signal_power /= signal.size(); + + // 使用信号功率的10%作为噪声估计 + float noise_power = signal_power * 0.1f; + + // 计算SNR + float snr = 0.0f; + if (noise_power > 0.0001f) { + snr = 10.0f * std::log10(signal_power / noise_power); + } + + return snr; } \ No newline at end of file diff --git a/app/src/main/cpp/src/signal_processor.cpp b/app/src/main/cpp/src/signal_processor.cpp index 3362222..d322104 100644 --- a/app/src/main/cpp/src/signal_processor.cpp +++ b/app/src/main/cpp/src/signal_processor.cpp @@ -12,157 +12,6 @@ T clamp(T value, T min_val, T max_val) { return value; } -// 新增:简化的通道级滤波处理(测试版本) -std::vector SignalProcessor::process_channel_based_filtering_simple( - const std::vector& data_packets) { - - std::cout << "=== 开始简化版通道级滤波处理 ===" << std::endl; - - if (data_packets.empty()) { - std::cout << "输入数据包为空,返回空结果" << std::endl; - return {}; - } - - std::cout << "输入数据包数量: " << data_packets.size() << std::endl; - - // 按数据类型分组 - std::map> grouped_data; - for (const auto& packet : data_packets) { - grouped_data[packet.data_type].push_back(packet); - } - - std::cout << "按数据类型分组完成,共 " << grouped_data.size() << " 种类型" << std::endl; - - std::vector processed_packets; - - // 对每种数据类型分别处理 - for (auto& [data_type, packets] : grouped_data) { - if (packets.empty()) continue; - - std::cout << "处理数据类型: " << static_cast(data_type) << ",数据包数量: " << packets.size() << std::endl; - - // 获取第一个数据包作为模板 - SensorData template_packet = packets[0]; - - // 检查通道数据类型 - if (!std::holds_alternative>>(packets[0].channel_data)) { - std::cout << "警告: 数据类型 " << static_cast(data_type) << " 不是多通道格式,跳过" << std::endl; - continue; - } - - // 获取通道信息 - auto& first_channels = std::get>>(packets[0].channel_data); - size_t num_channels = first_channels.size(); - size_t samples_per_packet = first_channels.empty() ? 0 : first_channels[0].size(); - - std::cout << "通道数量: " << num_channels << ", 每包采样点数: " << samples_per_packet << std::endl; - - if (num_channels == 0 || samples_per_packet == 0) { - std::cout << "警告: 通道数据无效,跳过此类型" << std::endl; - continue; - } - - // 计算总采样点数 - size_t total_samples = samples_per_packet * packets.size(); - std::cout << "总采样点数: " << total_samples << std::endl; - - // 创建合并后的数据包 - SensorData merged_packet = template_packet; - merged_packet.packet_sn = 0; // 合并后的包序号为0 - - // 创建合并后的通道数据 - auto& merged_channels = merged_packet.channel_data.emplace>>(); - merged_channels.resize(num_channels); - - // 合并所有数据包的通道数据 - for (size_t ch = 0; ch < num_channels; ch++) { - merged_channels[ch].reserve(total_samples); - - for (const auto& packet : packets) { - if (auto* channels = std::get_if>>(&packet.channel_data)) { - if (ch < channels->size() && !(*channels)[ch].empty()) { - merged_channels[ch].insert(merged_channels[ch].end(), - (*channels)[ch].begin(), - (*channels)[ch].end()); - } - } - } - - std::cout << "通道 " << ch << " 合并完成,采样点数: " << merged_channels[ch].size() << std::endl; - } - - // 对合并后的数据进行基本信号处理 - std::cout << "开始基本信号处理..." << std::endl; - for (size_t ch = 0; ch < num_channels; ch++) { - if (merged_channels[ch].empty()) continue; - - try { - // 根据数据类型应用不同的基本处理 - switch (data_type) { - case DataType::ECG_2LEAD: - // ECG 2导联基本处理:去除直流分量 - merged_channels[ch] = remove_dc_offset(merged_channels[ch]); - break; - case DataType::ECG_12LEAD: - // ECG 12导联专业滤波处理 - std::cout << "通道 " << ch << " 开始12导联心电专业滤波..." << std::endl; - - // 1. 去除直流分量 - merged_channels[ch] = remove_dc_offset(merged_channels[ch]); - std::cout << " - 直流分量去除完成" << std::endl; - - // 2. 0.1Hz高通滤波(去除基线漂移,比0.5Hz更温和) - merged_channels[ch] = filter(merged_channels[ch], 250.0, 0, 0.1, filtertype::highpass); - std::cout << " - 0.1Hz高通滤波完成" << std::endl; - - // 3. 50Hz陷波滤波(去除工频干扰,带宽1.0Hz更精确) - merged_channels[ch] = filter(merged_channels[ch], 250.0, 49.5, 50.5, filtertype::notchpass); - std::cout << " - 50Hz陷波滤波完成" << std::endl; - - std::cout << "通道 " << ch << " 12导联心电滤波处理完成" << std::endl; - break; - case DataType::EEG: - // EEG基本处理:去除直流分量 - merged_channels[ch] = remove_dc_offset(merged_channels[ch]); - break; - case DataType::PPG: - // PPG基本处理:去除直流分量 - merged_channels[ch] = remove_dc_offset(merged_channels[ch]); - // 2. 0.5-8Hz带通滤波(更精确的PPG频带) - merged_channels[ch] = bandpass_filter(merged_channels[ch], 50, 0.5, 8.0); - // // 3. 50Hz陷波滤波(去除工频干扰) - //merged_channels[ch] = filter(merged_channels[ch], 50, 49.5, 50.5, filtertype::notchpass); - // // 4. 运动伪迹检测和去除(优化版本) - merged_channels[ch] = remove_motion_artifacts(merged_channels[ch], 50); - break; - case DataType::RESPIRATION: - // 呼吸信号基本处理:去除直流分量 - merged_channels[ch] = remove_dc_offset(merged_channels[ch]); - break; - default: - // 通用处理:去除直流分量 - merged_channels[ch] = remove_dc_offset(merged_channels[ch]); - break; - } - - std::cout << "通道 " << ch << " 基本处理完成" << std::endl; - } catch (const std::exception& e) { - std::cerr << "通道 " << ch << " 处理失败: " << e.what() << std::endl; - // 继续处理其他通道 - } - } - - // 添加到处理结果中 - processed_packets.push_back(merged_packet); - - std::cout << "数据类型 " << static_cast(data_type) << " 处理完成" << std::endl; - } - - std::cout << "=== 简化版通道级滤波处理完成 ===" << std::endl; - std::cout << "总共创建 " << processed_packets.size() << " 个合并数据包" << std::endl; - - return processed_packets; -} // 新增:对完整通道数据应用滤波器 std::vector> SignalProcessor::apply_channel_filters( @@ -1688,3 +1537,87 @@ std::vector SignalProcessor::apply_smoothing( return smoothed_signal; } +// ============================================================================ +// 一阶滤波器实现 + +// 计算一阶滤波器参数 +FirstOrderFilterParams SignalProcessor::calculateFirstOrderFilterParams(double sample_rate, double cutoff_freq, filtertype type) { + if (sample_rate <= 0 || cutoff_freq <= 0 || cutoff_freq >= sample_rate / 2) { + // 返回单位滤波器参数(不滤波) + return FirstOrderFilterParams(1.0, 0.0, 0.0); + } + //α = exp(-ωc * T) = exp(-2π * fc / fs) + // 计算归一化频率 + double omega = 2.0 * M_PI * cutoff_freq / sample_rate; + double alpha = std::exp(-omega); + switch (type) { + case filtertype::lowpass: { + return FirstOrderFilterParams(1.0 - alpha, 0.0, -alpha); + } + + case filtertype::highpass: { // y[n] = b0*x[n] + b1*x[n-1] - a1*y[n-1] + double gain = (1.0 + alpha) / 2.0; + return FirstOrderFilterParams(gain, -gain, -alpha); + } + + case filtertype::bandpass: { + double gain = (1.0 + alpha) / 2.0; + return FirstOrderFilterParams(gain, -gain, -alpha); + } + + default: + return FirstOrderFilterParams(1.0, 0.0, 0.0); + } +} + +// 应用一阶滤波器 +std::vector SignalProcessor::applyFirstOrderFilter(const std::vector& input, const FirstOrderFilterParams& params) { + if (input.empty()) { + return input; + } + + std::vector output(input.size()); + + // 初始化状态变量 + double x1 = 0.0; // 输入延迟 + double y1 = 0.0; // 输出延迟 + // 应用差分方程:y[n] = b0*x[n] + b1*x[n-1] - a1*y[n-1] + for (size_t i = 0; i < input.size(); i++) { + double x0 = input[i]; + double y0 = params.b0 * x0 + params.b1 * x1 - params.a1 * y1; + + output[i] = static_cast(y0); + + // 更新延迟 + x1 = x0; + y1 = y0; + } + + return output; +} + +// 辅助函数:计算一阶带通滤波器参数(高通 + 低通) +std::pair +SignalProcessor::calculateFirstOrderBandpassParams(double sample_rate, double low_cutoff, double high_cutoff) { + // 高通滤波器参数 + FirstOrderFilterParams highpass_params = calculateFirstOrderFilterParams(sample_rate, low_cutoff, filtertype::highpass); + + // 低通滤波器参数 + FirstOrderFilterParams lowpass_params = calculateFirstOrderFilterParams(sample_rate, high_cutoff, filtertype::lowpass); + + return std::make_pair(highpass_params, lowpass_params); +} + +// 应用一阶带通滤波器(级联高通和低通) +std::vector SignalProcessor::applyFirstOrderBandpassFilter(const std::vector& input, + const FirstOrderFilterParams& highpass_params, + const FirstOrderFilterParams& lowpass_params) { + // 先应用高通滤波 + std::vector highpass_output = applyFirstOrderFilter(input, highpass_params); + + // 再应用低通滤波 + std::vector bandpass_output = applyFirstOrderFilter(highpass_output, lowpass_params); + + return bandpass_output; +} +