This commit is contained in:
ZhangJinLong 2025-10-10 10:33:05 +08:00
parent 1ce6cbffb3
commit 767d182d01
6 changed files with 295 additions and 172 deletions

View File

@ -42,6 +42,10 @@ public:
float calculate_signal_quality(const std::vector<float>& signal); float calculate_signal_quality(const std::vector<float>& signal);
float estimate_heart_rate_from_frequency(const std::vector<float>& signal, float sample_rate); float estimate_heart_rate_from_frequency(const std::vector<float>& signal, float sample_rate);
// 频域分析相关函数
float calculate_snr_frequency_domain(const std::vector<float>& signal);
float calculate_snr_time_domain(const std::vector<float>& signal);
// 综合指标计算 // 综合指标计算
std::map<std::string, float> calculate_all_ecg_metrics(const SensorData& ecg_data, float sample_rate); std::map<std::string, float> calculate_all_ecg_metrics(const SensorData& ecg_data, float sample_rate);
std::map<std::string, float> calculate_all_ppg_metrics(const SensorData& ppg_data, float sample_rate); std::map<std::string, float> calculate_all_ppg_metrics(const SensorData& ppg_data, float sample_rate);

View File

@ -10,6 +10,16 @@ enum class filtertype
{ {
lowpass,highpass,notchpass,bandpass,bandstop 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 { struct ECGChannelFeatures {
std::vector<int> r_peaks; // R波位置 std::vector<int> r_peaks; // R波位置
@ -82,6 +92,13 @@ public:
// 设置滤波器参数 // 设置滤波器参数
void setFilterParams(const std::string& filter_type, const FilterConfig& params); void setFilterParams(const std::string& filter_type, const FilterConfig& params);
// 一阶滤波器相关函数
FirstOrderFilterParams calculateFirstOrderFilterParams(double sample_rate, double cutoff_freq, filtertype type);
std::vector<float> applyFirstOrderFilter(const std::vector<float>& input, const FirstOrderFilterParams& params);
// 一阶带通滤波器辅助函数
std::pair<FirstOrderFilterParams, FirstOrderFilterParams> calculateFirstOrderBandpassParams(double sample_rate, double low_cutoff, double high_cutoff);
std::vector<float> applyFirstOrderBandpassFilter(const std::vector<float>& input, const FirstOrderFilterParams& highpass_params, const FirstOrderFilterParams& lowpass_params);
// 修改为public访问权限 // 修改为public访问权限
public: public:

View File

@ -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; return result;
} }

View File

@ -979,38 +979,37 @@ std::vector<float> MetricsCalculator::detect_pulse_peaks(const std::vector<float
return pulse_peaks; return pulse_peaks;
} }
// 新增:信号质量评估 // 新增:信号质量评估(基于频域分析)
float MetricsCalculator::calculate_signal_quality(const std::vector<float>& signal) { float MetricsCalculator::calculate_signal_quality(const std::vector<float>& signal) {
if (signal.empty()) return 0.0f; if (signal.empty()) return 0.0f;
const size_t n = signal.size(); const size_t n = signal.size();
// 1. 计算信噪比 (SNR) // 1. 计算信噪比 (SNR) - 使用频域分析
float signal_power = 0.0f; float snr = calculate_snr_frequency_domain(signal);
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);
}
// 2. 计算信号连续性 // 2. 计算信号连续性
float continuity_score = 0.0f; float continuity_score = 0.0f;
int discontinuity_count = 0; 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++) { for (size_t i = 1; i < n; i++) {
float diff = std::abs(signal[i] - signal[i-1]); 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++; discontinuity_count++;
} }
} }
@ -1200,3 +1199,173 @@ std::map<std::string, float> MetricsCalculator::calculate_all_hrv_metrics(const
return metrics; return metrics;
} }
// ============================================================================
// 频域分析相关函数
// 简单的FFT实现使用Cooley-Tukey算法
std::vector<std::complex<float>> fft(const std::vector<float>& signal) {
const size_t n = signal.size();
// 确保信号长度是2的幂次
size_t padded_size = 1;
while (padded_size < n) {
padded_size <<= 1;
}
std::vector<std::complex<float>> padded_signal(padded_size, 0.0f);
for (size_t i = 0; i < n; i++) {
padded_signal[i] = std::complex<float>(signal[i], 0.0f);
}
// Cooley-Tukey FFT
std::vector<std::complex<float>> 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<float> wlen(std::cos(angle), std::sin(angle));
for (size_t i = 0; i < padded_size; i += len) {
std::complex<float> w(1.0f, 0.0f);
for (size_t j = 0; j < len / 2; j++) {
std::complex<float> u = result[i + j];
std::complex<float> 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<float> calculate_power_spectrum(const std::vector<std::complex<float>>& fft_result) {
std::vector<float> 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<float>& 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<size_t>(low_freq / freq_resolution);
size_t high_bin = static_cast<size_t>(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<float>& 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<float>& 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;
}

View File

@ -12,157 +12,6 @@ T clamp(T value, T min_val, T max_val) {
return value; return value;
} }
// 新增:简化的通道级滤波处理(测试版本)
std::vector<SensorData> SignalProcessor::process_channel_based_filtering_simple(
const std::vector<SensorData>& data_packets) {
std::cout << "=== 开始简化版通道级滤波处理 ===" << std::endl;
if (data_packets.empty()) {
std::cout << "输入数据包为空,返回空结果" << std::endl;
return {};
}
std::cout << "输入数据包数量: " << data_packets.size() << std::endl;
// 按数据类型分组
std::map<DataType, std::vector<SensorData>> 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<SensorData> processed_packets;
// 对每种数据类型分别处理
for (auto& [data_type, packets] : grouped_data) {
if (packets.empty()) continue;
std::cout << "处理数据类型: " << static_cast<int>(data_type) << ",数据包数量: " << packets.size() << std::endl;
// 获取第一个数据包作为模板
SensorData template_packet = packets[0];
// 检查通道数据类型
if (!std::holds_alternative<std::vector<std::vector<float>>>(packets[0].channel_data)) {
std::cout << "警告: 数据类型 " << static_cast<int>(data_type) << " 不是多通道格式,跳过" << std::endl;
continue;
}
// 获取通道信息
auto& first_channels = std::get<std::vector<std::vector<float>>>(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<std::vector<std::vector<float>>>();
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<std::vector<std::vector<float>>>(&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<int>(data_type) << " 处理完成" << std::endl;
}
std::cout << "=== 简化版通道级滤波处理完成 ===" << std::endl;
std::cout << "总共创建 " << processed_packets.size() << " 个合并数据包" << std::endl;
return processed_packets;
}
// 新增:对完整通道数据应用滤波器 // 新增:对完整通道数据应用滤波器
std::vector<std::vector<float>> SignalProcessor::apply_channel_filters( std::vector<std::vector<float>> SignalProcessor::apply_channel_filters(
@ -1688,3 +1537,87 @@ std::vector<float> SignalProcessor::apply_smoothing(
return smoothed_signal; 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<float> SignalProcessor::applyFirstOrderFilter(const std::vector<float>& input, const FirstOrderFilterParams& params) {
if (input.empty()) {
return input;
}
std::vector<float> 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<float>(y0);
// 更新延迟
x1 = x0;
y1 = y0;
}
return output;
}
// 辅助函数:计算一阶带通滤波器参数(高通 + 低通)
std::pair<FirstOrderFilterParams, FirstOrderFilterParams>
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<float> SignalProcessor::applyFirstOrderBandpassFilter(const std::vector<float>& input,
const FirstOrderFilterParams& highpass_params,
const FirstOrderFilterParams& lowpass_params) {
// 先应用高通滤波
std::vector<float> highpass_output = applyFirstOrderFilter(input, highpass_params);
// 再应用低通滤波
std::vector<float> bandpass_output = applyFirstOrderFilter(highpass_output, lowpass_params);
return bandpass_output;
}