geopro/external/gpr3dviewer/RadarProcessor.cpp

/*
 * RadarProcessor.cpp
 * GPR探地雷达预处理算法实现源文件
 * 配套头文件 RadarProcessor.h
 * 依赖库：kiss_fftr(快速傅里叶变换)、Qt容器、标准数学库
 * 整体流水线顺序：零时校正 → 零漂消除 → 背景压制 → TVG深度增益 → 带通滤波 → 剖面平滑 → 道内AGC → 道间均衡 → 希尔伯特包络/相位
 * 数据存储格式：short 16位整型振幅，全程运算double浮点，截断回short防止溢出
 */
#include "RadarProcessor.h"
#include "GPRDataModel.h"
#include "PerformanceLogger.h"
#include <algorithm>
#include <cmath>
#include <limits>
#include <QVector>
#include <vector>
#include "kiss_fftr.h"
#include <omp.h>

// 匿名命名空间：仅内部工具函数，对外不可见
namespace {

/**
 * @brief 将浮点运算结果限制在short量程 [-32768,32767]
 * @param value 运算后浮点振幅值
 * @return 钳位取整后的16位短整型
 */
short clampToShort(double value) {
    if (value > 32767.0) return 32767;
    if (value < -32768.0) return -32768;
    return static_cast<short>(std::lround(value));
}

/**
 * @brief 计算单道波形整体RMS均方根振幅
 * @param trace 单道雷达波形数据
 * @return 该道RMS能量值
 */
double traceRms(const RadarTrace &trace) {
    if (trace.amplitudes.isEmpty()) return 0.0;
    double sumSquares = 0.0;
    for (short amp : trace.amplitudes) {
        double v = amp;
        sumSquares += v * v;
    }
    return std::sqrt(sumSquares / trace.amplitudes.size());
}

/**
 * @brief 全局整份数据所有采样点总RMS
 * @param model 完整GPR多道数据集
 * @return 全局平均能量
 */
double globalRms(const GPRDataModel &model) {
    double sumSq = 0.0;
    int cnt = 0;
    for (const auto& tr : model.traces) {
        for (short a : tr.amplitudes) {
            double v = a;
            sumSq += v * v;
            cnt++;
        }
    }
    if (cnt == 0) return 0.0;
    return std::sqrt(sumSq / cnt);
}

/**
 * @brief 安全窗口尺寸限制，防止窗口超出数据长度
 * @param req 请求窗口大小
 * @param limit 最大允许长度
 * @return 合法区间[1,limit]内窗口值
 */
int safeWindowSize(int req, int limit) {
    if (limit <= 0) return 0;
    return std::clamp(req, 1, limit);
}

/**
 * @brief 单道自动识别直达波零点起跳位置
 * @param trace 单道波形
 * @param frontWindow 仅在道前N个采样内搜索零点
 * @param sigmaMultiple 噪声标准差倍数起跳阈值
 * @return 零点对应的采样下标，-1代表识别失败
 * 逻辑：取道最前端小段计算基线噪声，使用正负极性对称的振幅+梯度条件判定起跳点
 */
int detectSingleTraceZeroPoint(const RadarTrace &trace, int frontWindow, double sigmaMultiple)
{
    const int sampleCount = trace.amplitudes.size();
    const int searchEnd = std::min(frontWindow, sampleCount);
    // 搜索区间太短无法判定基线噪声，直接返回失败
    if (searchEnd < 10) return -1;

    // 取搜索窗前1/4且不少于10点作为安静基线，避免固定10点对噪声估计过敏
    const int baselineLength = std::clamp(searchEnd / 4, 10, std::min(searchEnd, 40));
    double mean = 0.0;
    for (int i = 0; i < baselineLength; ++i) {
        mean += trace.amplitudes[i];
    }
    mean /= baselineLength;

    // 计算基线方差、标准差
    double variance = 0.0;
    for (int i = 0; i < baselineLength; ++i) {
        const double diff = trace.amplitudes[i] - mean;
        variance += diff * diff;
    }
    const double sigma = std::sqrt(variance / baselineLength);
    const double threshold = std::max(1.0, sigmaMultiple * std::max(sigma, 1.0));
    const double gradientThreshold = std::max(1.0, 0.5 * threshold);

    // 逐点扫描梯度+振幅，正负极性直达波均可识别
    for (int sample = 1; sample < searchEnd; ++sample) {
        const double current = trace.amplitudes[sample] - mean;
        const double previous = trace.amplitudes[sample - 1] - mean;
        const double gradient = current - previous;
        if (std::abs(gradient) > gradientThreshold && std::abs(current) > threshold) {
            return sample;
        }
    }
    return -1;
}

/**
 * @brief 升余弦过渡带频域响应函数（带通滤波器窗函数）
 * @param freq 当前频率点Hz
 * @param low 通带下限Hz
 * @param high 通带上限Hz
 * @param transition 高低频过渡带宽Hz
 * @return 频率增益系数[0,1]
 * 作用：平缓滚降，消除矩形窗吉布斯振荡，滤波波形畸变更小
 */
double bandpassResponse(double freq, double low, double high, double transition)
{
    // 升余弦渐入：低频侧从阻带到通带平滑抬升
    auto raisedCosineFadeIn = [](double t) {
        return 0.5 * (1.0 - std::cos(M_PI * t));
    };
    // 升余弦渐出：高频侧从通带到阻带平滑衰减
    auto raisedCosineFadeOut = [](double t) {
        return 0.5 * (1.0 + std::cos(M_PI * t));
    };

    // 完全阻带，增益0
    if (freq < low - transition || freq > high + transition) return 0.0;
    // 完全通带，增益1
    if (freq >= low + transition && freq <= high - transition) return 1.0;

    // 低频过渡区
    if (freq < low + transition) {
        double t = (freq - (low - transition)) / (2.0 * transition);
        t = std::clamp(t, 0.0, 1.0);
        return raisedCosineFadeIn(t);
    }
    // 高频过渡区
    else {
        double t = (freq - (high - transition)) / (2.0 * transition);
        t = std::clamp(t, 0.0, 1.0);
        return raisedCosineFadeOut(t);
    }
}

} // end 匿名命名空间

/**
 * @brief 处理器构造函数，无资源初始化
 */
RadarProcessor::RadarProcessor()
{
}

//=====================================================================
// 步骤0：零时自动/手动裁切
//=====================================================================
/**
 * @brief 切除每道零点前无效采样，统一道起始位置
 * @param input 原始输入数据
 * @param param 零时校正配置参数
 * @return 裁切后数据集
 * 自动模式：统计所有道零点中位数作为统一裁切长度；手动模式直接使用配置cutSamples
 */
GPRDataModel RadarProcessor::cutTimeZero(const GPRDataModel& input, const TimeZeroCutParams& param)
{
    SCOPED_PERF_TIMER("Processing", "cutTimeZero");
    // 开关未开启直接返回原数据
    if (!param.enable) return input;
    GPRDataModel out = input;
    int cutN = std::max(0, param.cutSamples);

    // 自动零点识别模式
    if (param.autoDetect) {
        QVector<int> zeroPoints;
        zeroPoints.reserve(out.traces.size());
        const int nChannels = std::max(1, out.header.numberOfChannels);
        QVector<QVector<int>> channelZeroPoints(nChannels);

        // 逐道检测零点，失败值(-1)不参与统计，避免把全局中位数拉向0
        for (int ti = 0; ti < out.traces.size(); ++ti) {
            const int zeroPoint = detectSingleTraceZeroPoint(out.traces[ti], param.frontSearchWindow, param.noiseSigmaMultiple);
            if (zeroPoint <= 0) continue;
            zeroPoints.append(zeroPoint);
            channelZeroPoints[ti % nChannels].append(zeroPoint);
        }

        QVector<int> channelMedians;
        channelMedians.reserve(nChannels);
        for (QVector<int> &points : channelZeroPoints) {
            if (points.isEmpty()) continue;
            std::sort(points.begin(), points.end());
            channelMedians.append(points[points.size() / 2]);
        }

        if (!channelMedians.isEmpty()) {
            // 多通道数据先按通道求中位数，再对通道中位数取中位数，避免通道间耦合差异互相污染
            std::sort(channelMedians.begin(), channelMedians.end());
            cutN = channelMedians[channelMedians.size() / 2];
        } else if (!zeroPoints.isEmpty()) {
            std::sort(zeroPoints.begin(), zeroPoints.end());
            cutN = zeroPoints[zeroPoints.size() / 2];
        }
    }

    // 无需裁切
    if (cutN <= 0) return out;

    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : out.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }

    // 每道截断前面cutN个采样点
    RadarTrace *traces = out.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < out.traces.size(); ++ti)
    {
        RadarTrace &trace = traces[ti];
        int total = trace.amplitudes.size();
        if (cutN >= total) continue;
        trace.amplitudes = trace.amplitudes.mid(cutN);
    }
    // 更新文件头每道采样总数
    if (!out.traces.isEmpty())
    {
        out.header.samplesPerTrace = out.traces.first().amplitudes.size();
    }
    return out;
}

//=====================================================================
// 简易整道去直流（零漂DC模式底层实现）
//=====================================================================
/**
 * @brief 扣除单道整体直流均值基线
 * @param input 输入数据
 * @param param 开关配置
 * @return 去直流后数据
 */
GPRDataModel RadarProcessor::removeDCShift(const GPRDataModel &input, const DcShiftParams& param)
{
    SCOPED_PERF_TIMER("Processing", "removeDCShift");
    if (!param.enable) return input;
    GPRDataModel output = input;
    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }
    RadarTrace *traces = output.traces.data();
    #pragma omp parallel for
    for (int i = 0; i < output.traces.size(); ++i) {
        RadarTrace &trace = traces[i];
        if (trace.amplitudes.isEmpty()) continue;
        // 计算单道整体均值
        double sum = 0;
        for (auto amp : trace.amplitudes) sum += amp;
        double mean = sum / trace.amplitudes.size();
        // 逐采样减去均值并钳位short
        short *amps = trace.amplitudes.data();
        int sampleCount = trace.amplitudes.size();
        for (int s = 0; s < sampleCount; ++s) {
            amps[s] = clampToShort(amps[s] - mean);
        }
    }
    return output;
}

//=====================================================================
// 零漂消除：DC全局均值 / Sliding滑动窗口基线
//=====================================================================
/**
 * @brief 消除波形基线漂移
 * @param input 输入数据
 * @param param 模式与窗口参数
 * @return 校正后数据
 * DC模式直接复用removeDCShift；滑动窗口用前缀和快速计算局部均值基线
 */
GPRDataModel RadarProcessor::removeZeroDrift(const GPRDataModel &input, const ZeroDriftParams& param)
{
    SCOPED_PERF_TIMER("Processing", "removeZeroDrift");
    if (!param.enable) return input;
    // 直流均值模式，调用简易去直流接口
    if (param.mode == ZeroDriftMode::DC) {
        DcShiftParams dc;
        dc.enable = true;
        return removeDCShift(input, dc);
    }

    GPRDataModel output = input;
    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }
    RadarTrace *traces = output.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < output.traces.size(); ++ti) {
        RadarTrace &trace = traces[ti];
        const QVector<short> source = trace.amplitudes;
        const int sampleCount = source.size();
        if (sampleCount == 0) continue;

        const int window = safeWindowSize(param.slidingWindowSamples, sampleCount);
        const int halfWindow = window / 2;
        // 前缀和数组，O(1)区间求和加速滑动窗口均值
        QVector<double> prefix(sampleCount + 1, 0.0);
        for (int sample = 0; sample < sampleCount; ++sample) {
            prefix[sample + 1] = prefix[sample] + source[sample];
        }

        // 逐采样计算窗口内局部基线并扣除
        short *amps = trace.amplitudes.data();
        for (int sample = 0; sample < sampleCount; ++sample) {
            const int start = std::max(0, sample - halfWindow);
            const int end = std::min(sampleCount, start + window);
            const int adjStart = std::max(0, end - window);
            const int count = end - adjStart;
            const double localMean = count > 0 ? (prefix[end] - prefix[adjStart]) / count : 0.0;
            amps[sample] = clampToShort(source[sample] - localMean);
        }
    }
    return output;
}

//=====================================================================
// 全局统一倍率增益
//=====================================================================
/**
 * @brief 整份数据统一放大缩小倍数
 * @param input 输入数据
 * @param param 增益系数与开关
 * @return 缩放后数据
 */
GPRDataModel RadarProcessor::applyGain(const GPRDataModel &input, const GlobalGainParams& param)
{
    SCOPED_PERF_TIMER("Processing", "applyGain");
    if (!param.enable) return input;
    GPRDataModel output = input;
    double f = param.factor;
    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }
    RadarTrace *traces = output.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < output.traces.size(); ++ti) {
        RadarTrace &tr = traces[ti];
        short *amps = tr.amplitudes.data();
        int sampleCount = tr.amplitudes.size();
        for (int s = 0; s < sampleCount; ++s) {
            amps[s] = clampToShort(amps[s] * f);
        }
    }
    return output;
}

//=====================================================================
// TVG 球面扩散+介质吸收深度增益补偿
//=====================================================================
/**
 * @brief 随深度自动补偿电磁波衰减
 * @param input 原始数据
 * @param params 速度、吸收系数、最大增益限制等配置
 * @return 增益补偿后剖面
 * 公式：增益 = 球面扩散项 × 指数吸收项，上限maxGain防止噪声爆炸放大
 */
GPRDataModel RadarProcessor::sphericalTvg(const GPRDataModel &input, const SphericalTvgParams &params)
{
    SCOPED_PERF_TIMER("Processing", "sphericalTvg");
    if (!params.enable) return input;
    GPRDataModel output = input;
    if (output.traces.isEmpty() || output.header.samplesPerTrace <= 0 || output.header.timeWindowNs <= 0.0) return output;

    // 波速优先级：参数指定 > 文件头读取 > 兜底0.1 m/ns
    const double velocity = params.velocityMPerNs > 0.0 ? params.velocityMPerNs : (output.header.waveVelocity > 0.0 ? output.header.waveVelocity : 0.1);
    const double referenceDepth = std::max(params.referenceDepthM, std::numeric_limits<double>::epsilon());
    const double exponent = std::max(params.exponent, 0.0);
    const double maxGain = params.maxGain > 0.0 ? params.maxGain : 1.0;
    // 单采样时间间隔ns
    const double dtNs = output.header.timeWindowNs / output.header.samplesPerTrace;

    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }
    RadarTrace *traces = output.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < output.traces.size(); ++ti) {
        RadarTrace &trace = traces[ti];
        short *amps = trace.amplitudes.data();
        const int sampleCount = trace.amplitudes.size();
        for (int sample = 0; sample < sampleCount; ++sample) {
            const double timeNs = sample * dtNs;
            // 单程深度 = 往返时间 × 波速 / 2
            const double depthM = timeNs * velocity / 2.0;
            double gain = 1.0;
            // 球面扩散补偿 (r/r0)^n
            if (params.enableSpherical) {
                gain *= std::pow(std::max(depthM, referenceDepth) / referenceDepth, exponent);
            }
            // 介质指数吸收补偿 exp(β·t)
            if (params.enableAbsorption && params.absorptionBeta > 0.0) {
                gain *= std::exp(params.absorptionBeta * timeNs);
            }
            // 增益封顶
            gain = std::min(gain, maxGain);
            amps[sample] = clampToShort(amps[sample] * gain);
        }
    }
    return output;
}

//=====================================================================
// 二维剖面平滑：垂直(深度) + 水平(测线道)均值平滑
//=====================================================================
/**
 * @brief 时空二维滑动均值降噪
 * @param input 输入数据集
 * @param params 窗口半宽、横竖平滑开关
 * @return 平滑后剖面
 * 先深度方向逐道平滑，再道方向同通道横向平滑；横向使用原始未平滑数据做基准避免二次模糊
 */
GPRDataModel RadarProcessor::profileSmooth(const GPRDataModel &input, const ProfileSmoothParams &params)
{
    SCOPED_PERF_TIMER("Processing", "profileSmooth");
    if (!params.enable) return input;
    GPRDataModel output = input;
    if (output.traces.isEmpty()) return output;

    const int nChannels = output.header.numberOfChannels > 0 ? output.header.numberOfChannels : 1;
    const int tracesPerChannel = output.getTracesPerChannel();
    const int winHalf = std::max(0, params.smoothWindow);
    if (winHalf == 0) return output;

    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }

    // 第一步：垂直深度方向单道内平滑
    if (params.verticalSmooth) {
        RadarTrace *traces = output.traces.data();
        #pragma omp parallel for
        for (int ti = 0; ti < output.traces.size(); ++ti) {
            RadarTrace &trace = traces[ti];
            const QVector<short> source = trace.amplitudes;
            const int sampleCount = source.size();
            short *amps = trace.amplitudes.data();
            for (int sample = 0; sample < sampleCount; ++sample) {
                // 窗口左右边界钳位
                const int start = std::max(0, sample - winHalf);
                const int end = std::min(sampleCount - 1, sample + winHalf);
                double sum = 0.0;
                for (int idx = start; idx <= end; ++idx) {
                    sum += source[idx];
                }
                amps[sample] = clampToShort(sum / (end - start + 1));
            }
        }
    }

    // 第二步：水平测线方向同通道多道平滑，使用原始模型做参考源
    if (params.horizontalSmooth && tracesPerChannel > 1) {
        const GPRDataModel sourceModel = output;
        RadarTrace *traces = output.traces.data();
        #pragma omp parallel for
        for (int globalIndex = 0; globalIndex < output.traces.size(); ++globalIndex) {
            int channel = globalIndex % nChannels;
            int traceInChannel = globalIndex / nChannels;

            RadarTrace &targetTrace = traces[globalIndex];
            const int sampleCount = targetTrace.amplitudes.size();
            const int startTrace = std::max(0, traceInChannel - winHalf);
            const int endTrace = std::min(tracesPerChannel - 1, traceInChannel + winHalf);
            // 每个采样点横向多道平均
            for (int sample = 0; sample < sampleCount; ++sample) {
                double sum = 0.0;
                int count = 0;
                for (int neighbor = startTrace; neighbor <= endTrace; ++neighbor) {
                    const int neighborIndex = sourceModel.getTraceIndex(channel, neighbor);
                    if (neighborIndex < 0 || neighborIndex >= sourceModel.traces.size()) continue;
                    const QVector<short> &amps = sourceModel.traces[neighborIndex].amplitudes;
                    if (sample >= amps.size()) continue;
                    sum += amps[sample];
                    ++count;
                }
                if (count > 0) {
                    short *tgtAmps = targetTrace.amplitudes.data();
                    tgtAmps[sample] = clampToShort(sum / count);
                }
            }
        }
    }
    return output;
}

//=====================================================================
// 道内滑动窗口AGC振幅均衡
//=====================================================================
/**
 * @brief 单道内部深浅振幅自适应均衡
 * @param input 原始数据
 * @param params 窗口大小、目标RMS、最大增益限制
 * @return AGC均衡后数据
 * 前缀平方和快速计算局部RMS，增益=目标RMS/局部RMS，封顶maxGain防止噪声放大
 */
GPRDataModel RadarProcessor::traceInnerAgc(const GPRDataModel &input, const TraceInnerAgcParams &params)
{
    SCOPED_PERF_TIMER("Processing", "traceInnerAgc");
    if (!params.enable) return input;
    GPRDataModel output = input;
    if (output.traces.isEmpty()) return output;

    // 目标RMS：0则自动取整份数据全局RMS
    const double targetRms = params.targetRms > 0.0 ? params.targetRms : globalRms(input);
    if (targetRms <= 0.0) return output;
    const double epsilon = params.epsilon > 0.0 ? params.epsilon : 1.0;
    const double maxGain = params.maxGain > 0.0 ? params.maxGain : 1.0;

    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }
    const RadarTrace *inputTraces = input.traces.data();
    RadarTrace *outputTraces = output.traces.data();
    #pragma omp parallel for
    for (int traceIndex = 0; traceIndex < input.traces.size(); ++traceIndex) {
        const RadarTrace &sourceTrace = inputTraces[traceIndex];
        RadarTrace &targetTrace = outputTraces[traceIndex];
        const int sampleCount = sourceTrace.amplitudes.size();
        if (sampleCount == 0) continue;

        const int window = safeWindowSize(params.windowSamples, sampleCount);
        const int halfWindow = window / 2;
        // 平方前缀和数组，快速区间RMS计算
        QVector<double> prefixSquares(sampleCount + 1, 0.0);
        const short *srcAmps = sourceTrace.amplitudes.data();
        for (int s = 0; s < sampleCount; s++) {
            double v = srcAmps[s];
            prefixSquares[s+1] = prefixSquares[s] + v * v;
        }
        // 逐采样计算局部增益并缩放振幅
        short *tgtAmps = targetTrace.amplitudes.data();
        for (int s = 0; s < sampleCount; s++) {
            int start = std::max(0, s - halfWindow);
            int end = std::min(sampleCount, start + window);
            int adjStart = std::max(0, end - window);
            int cnt = end - adjStart;
            double sumSq = prefixSquares[end] - prefixSquares[adjStart];
            double localRms = cnt > 0 ? std::sqrt(sumSq / cnt) : 0.0;
            // 防除零极小值epsilon
            double gain = std::min(targetRms / std::max(localRms, epsilon), maxGain);
            tgtAmps[s] = clampToShort(srcAmps[s] * gain);
        }
    }
    return output;
}

//=====================================================================
// 道间AGC：整条测线道与道振幅均衡
//=====================================================================
/**
 * @brief 横向测线多道之间整体能量拉平
 * @param input 输入数据
 * @param params 全局/局部滑动窗口模式
 * @return 均衡后剖面
 * Global：全部道共用一个基准RMS；Local：相邻窗口道平均RMS做基准
 */
GPRDataModel RadarProcessor::traceToTraceEqualization(const GPRDataModel &input, const TraceToTraceAgcParams &params)
{
    SCOPED_PERF_TIMER("Processing", "traceToTraceEqualization");
    if (!params.enable) return input;
    GPRDataModel output = input;
    int traceCnt = input.traces.size();
    if (traceCnt == 0) return output;

    const int nChannels = std::max(1, input.header.numberOfChannels);
    const int tracesPerChannel = input.getTracesPerChannel();

    // 预计算每道整体RMS能量
    QVector<double> trRms(traceCnt, 0.0);
    for (int i = 0; i < traceCnt; i++) trRms[i] = traceRms(input.traces[i]);
    double targetRms = params.targetRms > 0 ? params.targetRms : globalRms(input);
    if (targetRms <= 0) return output;

    double eps = params.epsilon > 0 ? params.epsilon : 1.0;
    double maxG = params.maxGain > 0 ? params.maxGain : 1.0;
    int win = safeWindowSize(params.horizontalWindowTraces, std::max(1, tracesPerChannel));
    int halfWin = win / 2;

    QVector<QVector<double>> channelPrefixSq(nChannels);
    if (params.mode == TraceToTraceMode::Local && tracesPerChannel > 0) {
        for (int ch = 0; ch < nChannels; ++ch) {
            channelPrefixSq[ch].resize(tracesPerChannel + 1);
            channelPrefixSq[ch][0] = 0.0;
            for (int tr = 0; tr < tracesPerChannel; ++tr) {
                const int idx = input.getTraceIndex(ch, tr);
                const double rms = (idx >= 0 && idx < traceCnt) ? trRms[idx] : 0.0;
                channelPrefixSq[ch][tr + 1] = channelPrefixSq[ch][tr] + rms * rms;
            }
        }
    }

    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }
    RadarTrace *traces = output.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < traceCnt; ti++)
    {
        double base = trRms[ti];
        // 局部模式只沿同一通道的测线方向取滑动窗口RMS，避免多通道交错数据互相污染
        if (params.mode == TraceToTraceMode::Local && tracesPerChannel > 0)
        {
            const int channel = ti % nChannels;
            const int traceInChannel = ti / nChannels;
            if (channel >= 0 && channel < channelPrefixSq.size() && traceInChannel >= 0 && traceInChannel < tracesPerChannel) {
                int st = std::max(0, traceInChannel - halfWin);
                int ed = std::min(tracesPerChannel, st + win);
                int adjSt = std::max(0, ed - win);
                int c = ed - adjSt;
                double sq = channelPrefixSq[channel][ed] - channelPrefixSq[channel][adjSt];
                base = c > 0 ? std::sqrt(sq / c) : base;
            }
        }
        // 计算整道统一增益
        double g = std::min(targetRms / std::max(base, eps), maxG);
        // 整道所有采样同倍率缩放
        short *amps = traces[ti].amplitudes.data();
        int sampleCount = traces[ti].amplitudes.size();
        for (int s = 0; s < sampleCount; ++s)
        {
            amps[s] = clampToShort(amps[s] * g);
        }
    }
    return output;
}

//=====================================================================
// 希尔伯特变换：振幅包络 / 正交虚部
//=====================================================================
/**
 * @brief 时域直接希尔伯特变换（简易实现，适合中小采样长度）
 * @param input 输入波形
 * @param params 输出包络/正交分量开关
 * @return 变换后数据
 * 公式：虚部h[n] = Σ(2/(π·(n-k)))·x[k]；包络=√(实部²+虚部²)
 */
GPRDataModel RadarProcessor::hilbertTransform(const GPRDataModel &input, const HilbertParams &params)
{
    SCOPED_PERF_TIMER("Processing", "hilbertTransform");
    if (!params.enable) return input;
    GPRDataModel output = input;
    if (output.traces.isEmpty()) return output;

    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : output.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }
    RadarTrace *traces = output.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < output.traces.size(); ++ti) {
        RadarTrace &trace = traces[ti];
        const QVector<short> source = trace.amplitudes;
        const int sampleCount = source.size();
        if (sampleCount < 3) continue;

        QVector<double> hilbert(sampleCount, 0.0);
        // 时域卷积计算希尔伯特虚部
        for (int n = 0; n < sampleCount; ++n) {
            double sum = 0.0;
            for (int k = 0; k < sampleCount; ++k) {
                const int diff = n - k;
                // 差分0、偶数项系数为0，只累加奇数差分
                if (diff == 0 || diff % 2 == 0) continue;
                sum += (2.0 / (M_PI * diff)) * source[k];
            }
            hilbert[n] = sum;
        }

        // 输出模式分支
        short *amps = trace.amplitudes.data();
        for (int sample = 0; sample < sampleCount; ++sample) {
            if (params.computeEnvelope) {
                // 成像常用：振幅包络
                const double real = source[sample];
                amps[sample] = clampToShort(std::sqrt(real * real + hilbert[sample] * hilbert[sample]));
            } else {
                // 相位分析：仅输出正交虚部
                amps[sample] = clampToShort(hilbert[sample]);
            }
        }
    }
    return output;
}

//=====================================================================
// FFT频域带通滤波（kiss_fftr快速实数FFT）
//=====================================================================
/**
 * @brief 频域升余弦窗带通滤波
 * @param input 原始数据
 * @param params 自动/手动频带配置
 * @return 滤波后波形
 * 流程：时域→FFT频域乘增益响应→逆FFT回时域；补零2倍长度避免循环卷积混叠
 */
GPRDataModel RadarProcessor::bandpassFilter(const GPRDataModel &input, const BandpassParams &params)
{
    SCOPED_PERF_TIMER("Processing", "bandpassFilter");
    if (!params.enable) return input;
    GPRDataModel out = input;
    if (out.traces.isEmpty() || out.header.samplesPerTrace <= 0) return out;

    const double dtNs = out.header.timeWindowNs / out.header.samplesPerTrace;
    if (dtNs <= 0.0) return out;
    const double fs = 1e9 / dtNs;

    double lowFreq = params.lowFreqHz;
    double highFreq = params.highFreqHz;
    if (params.autoFreq && params.antennaFreqMHz > 0.0) {
        const double fcenter = params.antennaFreqMHz * 1e6;
        lowFreq = fcenter * 0.1;
        highFreq = fcenter * 2.5;
    }
    lowFreq = std::max(lowFreq, 0.0);
    highFreq = std::min(highFreq, fs * 0.5);
    if (lowFreq >= highFreq) return out;

    const double transition = std::max(50.0, (highFreq - lowFreq) * 0.1);
    const int sampleCount = out.header.samplesPerTrace;

    int nfft = 1;
    while (nfft < sampleCount * 2) nfft <<= 1;

    std::vector<kiss_fft_scalar> filterResponse(nfft / 2 + 1);
    for (int k = 0; k <= nfft / 2; ++k) {
        const double freq = k * fs / nfft;
        filterResponse[k] = static_cast<kiss_fft_scalar>(bandpassResponse(freq, lowFreq, highFreq, transition));
    }

    for (auto &trace : out.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }

    #pragma omp parallel
    {
        kiss_fftr_cfg fwd = kiss_fftr_alloc(nfft, 0, nullptr, nullptr);
        kiss_fftr_cfg inv = kiss_fftr_alloc(nfft, 1, nullptr, nullptr);
        std::vector<kiss_fft_scalar> timedata(nfft, 0.0f);
        std::vector<kiss_fft_cpx> freqdata(nfft / 2 + 1);

        RadarTrace *traces = out.traces.data();
        if (fwd && inv) {
            #pragma omp for
            for (int ti = 0; ti < out.traces.size(); ++ti) {
                RadarTrace &trace = traces[ti];
                if (trace.amplitudes.size() < 4) continue;

                short *amps = trace.amplitudes.data();
                for (int i = 0; i < sampleCount; ++i) {
                    timedata[i] = static_cast<kiss_fft_scalar>(amps[i]);
                }
                for (int i = sampleCount; i < nfft; ++i) {
                    timedata[i] = 0.0f;
                }

                kiss_fftr(fwd, timedata.data(), freqdata.data());

                for (int k = 0; k <= nfft / 2; ++k) {
                    freqdata[k].r *= filterResponse[k];
                    freqdata[k].i *= filterResponse[k];
                }

                kiss_fftri(inv, freqdata.data(), timedata.data());

                const double scale = 1.0 / nfft;
                for (int i = 0; i < sampleCount; ++i) {
                    amps[i] = clampToShort(timedata[i] * scale);
                }
            }
        }

        if (fwd) kiss_fftr_free(fwd);
        if (inv) kiss_fftr_free(inv);
    }
    return out;
}

//=====================================================================
// 背景杂波去除：均值法 / 奇异值过滤法
//=====================================================================
/**
 * @brief 多道平均背景相减，压制固定耦合、地层静态杂波
 * @param input 原始剖面
 * @param params 模式、平均窗口道数、奇异阈值
 * @return 去除背景后数据
 * MeanAverage：全部道参与背景平均；SingularityFilter：过滤高能异常道再平均，保护空洞/管线强反射
 */
GPRDataModel RadarProcessor::backgroundRemoval(const GPRDataModel &input, const BackgroundParams &params)
{
    SCOPED_PERF_TIMER("Processing", "backgroundRemoval");
    if (!params.enable) return input;
    GPRDataModel out = input;
    const int traceTotal = out.traces.size();
    if (traceTotal < 3) return out;

    const int nChannels = std::max(1, out.header.numberOfChannels);
    const int tracesPerChannel = out.getTracesPerChannel();
    if (tracesPerChannel < 3) return out;

    QVector<QVector<bool>> validTraceMask(nChannels);
    for (int ch = 0; ch < nChannels; ++ch) {
        validTraceMask[ch].fill(true, tracesPerChannel);
    }

    // 奇异滤波模式：每个通道内部筛掉能量过高的异常道，只用平稳道算背景
    if (params.mode == BackgroundMode::SingularityFilter) {
        for (int ch = 0; ch < nChannels; ++ch) {
            QVector<double> energies;
            energies.reserve(tracesPerChannel);
            for (int tr = 0; tr < tracesPerChannel; ++tr) {
                const int idx = input.getTraceIndex(ch, tr);
                if (idx >= 0 && idx < traceTotal) {
                    energies.append(traceRms(input.traces[idx]));
                }
            }
            if (energies.isEmpty()) continue;

            QVector<double> sortedEnergies = energies;
            std::sort(sortedEnergies.begin(), sortedEnergies.end());
            const double median = sortedEnergies[sortedEnergies.size() / 2];
            const double limit = median * std::max(params.singularityThreshold, 1.0);

            for (int tr = 0; tr < tracesPerChannel; ++tr) {
                const int idx = input.getTraceIndex(ch, tr);
                validTraceMask[ch][tr] = idx >= 0 && idx < traceTotal && traceRms(input.traces[idx]) <= limit;
            }
        }
    }

    // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach
    for (auto &trace : out.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }

    // 滑动半道窗口只沿同一通道的测线方向取背景，避免多通道交错数据互相污染
    const int halfWindow = std::max(1, params.averageTraceCount) / 2;
    RadarTrace *outTraces = out.traces.data();
    const RadarTrace *inTraces = input.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < traceTotal; ++ti)
    {
        const int channel = ti % nChannels;
        const int traceInChannel = ti / nChannels;
        if (channel < 0 || channel >= validTraceMask.size() || traceInChannel < 0 || traceInChannel >= tracesPerChannel) continue;

        const int left = std::max(0, traceInChannel - halfWindow);
        const int right = std::min(tracesPerChannel - 1, traceInChannel + halfWindow);
        const int sampleCount = outTraces[ti].amplitudes.size();
        QVector<double> background(sampleCount, 0.0);
        int backgroundCount = 0;

        for (int tr = left; tr <= right; ++tr) {
            if (!validTraceMask[channel][tr]) continue;
            const int idx = input.getTraceIndex(channel, tr);
            if (idx < 0 || idx >= traceTotal) continue;
            const RadarTrace &bgTrace = inTraces[idx];
            const short *bgAmps = bgTrace.amplitudes.data();
            int bgSize = bgTrace.amplitudes.size();
            for (int sample = 0; sample < sampleCount && sample < bgSize; ++sample) {
                background[sample] += bgAmps[sample];
            }
            ++backgroundCount;
        }
        if (backgroundCount == 0) continue;

        RadarTrace &currentTrace = outTraces[ti];
        const QVector<short> &source = inTraces[ti].amplitudes;
        short *curAmps = currentTrace.amplitudes.data();
        const short *srcAmps = source.data();
        int srcSize = source.size();
        for (int sample = 0; sample < sampleCount && sample < srcSize; ++sample)
        {
            const double averageBackground = background[sample] / backgroundCount;
            curAmps[sample] = clampToShort(srcAmps[sample] - averageBackground);
        }
    }
    return out;
}

//=====================================================================
// Kirchhoff偏移：x-t域双曲线求和
//=====================================================================
/**
 * @brief Kirchhoff偏移在x-t域实现
 * @param input 输入剖面数据
 * @param params 求和宽度（两侧道数）与波速参数
 * @return 偏移后数据
 * 算法流程：
 * 1. 对原始剖面逐道时间方向微分
 * 2. 对每个输出点(t0, x0)沿双曲线 t = sqrt(t0^2 + (dx*n/v)^2) 求和
 * 3. 加权因子采用1/t几何扩散补偿
 * 4. sumWidth < 128时预计算双曲线查表，避免重复sqrt运算
 */
GPRDataModel RadarProcessor::kirchhoffMigration(const GPRDataModel &input, const MigrationParams &params)
{
    SCOPED_PERF_TIMER("Processing", "kirchhoffMigration");
    if (!params.enable) return input;
    GPRDataModel out = input;
    int traceTotal = out.traces.size();
    if (traceTotal < 3) return out;
    int sampleCount = out.header.samplesPerTrace;
    if (sampleCount < 3) return out;

    const double dtNs = out.header.timeWindowNs / sampleCount;
    double dx = out.header.distanceInc;
    if (dx <= 0.0) dx = 0.01;
    double v = params.velocityMPerNs > 0.0 ? params.velocityMPerNs : out.header.waveVelocity;
    if (v <= 0.0) v = 0.1;

    const int sumWidth = std::max(1, params.sumWidth);

    // 1. 逐道时间方向微分，结果存为double浮点数组
    QVector<QVector<double>> diffData(traceTotal);
    for (int ti = 0; ti < traceTotal; ++ti) {
        diffData[ti].resize(sampleCount);
        const QVector<short> &src = out.traces[ti].amplitudes;
        if (sampleCount > 0) diffData[ti][0] = 0.0;
        for (int s = 1; s < sampleCount; ++s) {
            diffData[ti][s] = static_cast<double>(src[s]) - static_cast<double>(src[s - 1]);
        }
    }

    // 2. 当求和宽度小于128时，预计算双曲线采样位置查表（只进行一次）
    QVector<QVector<double>> hyperbolaTable;
    if (sumWidth < 128) {
        hyperbolaTable.resize(sampleCount);
        for (int t0 = 0; t0 < sampleCount; ++t0) {
            hyperbolaTable[t0].resize(sumWidth * 2 + 1);
            const double t0Ns = t0 * dtNs;
            for (int k = -sumWidth; k <= sumWidth; ++k) {
                const double x = k * dx;
                const double tNs = std::sqrt(t0Ns * t0Ns + (x / v) * (x / v));
                hyperbolaTable[t0][k + sumWidth] = tNs / dtNs;
            }
        }
    }

    // 预先detach所有trace的amplitudes
    for (auto &trace : out.traces) {
        if (!trace.amplitudes.isEmpty())
            trace.amplitudes.data();
    }

    // 3. 对每个输出点沿双曲线加权求和
    RadarTrace *traces = out.traces.data();
    #pragma omp parallel for
    for (int ti = 0; ti < traceTotal; ++ti) {
        QVector<short> result(sampleCount, 0);
        for (int s = 0; s < sampleCount; ++s) {
            double sum = 0.0;
            int validCount = 0;
            for (int k = -sumWidth; k <= sumWidth; ++k) {
                const int srcTrace = ti + k;
                if (srcTrace < 0 || srcTrace >= traceTotal) continue;

                double srcSampleD;
                if (sumWidth < 128 && !hyperbolaTable.isEmpty()) {
                    srcSampleD = hyperbolaTable[s][k + sumWidth];
                } else {
                    const double t0Ns = s * dtNs;
                    const double x = k * dx;
                    const double tNs = std::sqrt(t0Ns * t0Ns + (x / v) * (x / v));
                    srcSampleD = tNs / dtNs;
                }

                int srcSample = static_cast<int>(srcSampleD);
                double frac = srcSampleD - srcSample;
                if (srcSample < 0) {
                    srcSample = 0;
                    frac = 0.0;
                }
                if (srcSample >= sampleCount - 1) {
                    srcSample = sampleCount - 1;
                    frac = 0.0;
                }

                double amp = diffData[srcTrace][srcSample];
                if (frac > 0.0 && srcSample < sampleCount - 1) {
                    amp += frac * (diffData[srcTrace][srcSample + 1] - diffData[srcTrace][srcSample]);
                }

                // 1/t 几何扩散加权（t>0时）
                const double tNs = s * dtNs;
                if (tNs > 1e-12) {
                    amp /= tNs;
                }

                sum += amp;
                ++validCount;
            }
            if (validCount > 0) {
                result[s] = clampToShort(sum);
            }
        }
        traces[ti].amplitudes = std::move(result);
    }

    return out;
}

//=====================================================================
// 流水线总调度入口：按步骤队列顺序串行执行所有开启算法
//=====================================================================
/**
 * @brief 批量串行执行整套预处理流水线
 * @param rawData 未处理原始GPR数据
 * @param pipeline 有序步骤配置容器
 * @return 逐步骤处理完成后的成品数据
 */
GPRDataModel RadarProcessor::runPipeline(const GPRDataModel& rawData, const ProcPipeline& pipeline)
{
    return runPipeline(rawData, pipeline, nullptr);
}

GPRDataModel RadarProcessor::runPipeline(const GPRDataModel& rawData,
                                         const ProcPipeline& pipeline,
                                         const PipelineRuntimeContext *context)
{
    SCOPED_PERF_TIMER("Processing", "runPipeline(Total)");
    auto canceled = [context]() {
        return context && context->isCanceled && context->isCanceled();
    };
    auto stepName = [](ProcStepType type) -> QString {
        switch (type) {
        case ProcStepType::StepTimeZeroCut: return QStringLiteral("Set Time Zero");
        case ProcStepType::StepZeroDrift: return QStringLiteral("Zero Drift");
        case ProcStepType::StepRemoveDC: return QStringLiteral("Remove DC");
        case ProcStepType::StepBackgroundRemove: return QStringLiteral("Background Removal");
        case ProcStepType::StepSphericalTVG: return QStringLiteral("Spherical TVG");
        case ProcStepType::StepBandpassFilter: return QStringLiteral("Bandpass Filter");
        case ProcStepType::StepProfileSmooth: return QStringLiteral("Profile Smooth");
        case ProcStepType::StepInnerAGC: return QStringLiteral("Trace Inner AGC");
        case ProcStepType::StepTraceToTraceAGC: return QStringLiteral("Trace-to-Trace AGC");
        case ProcStepType::StepHilbertTransform: return QStringLiteral("Hilbert Transform");
        case ProcStepType::StepGlobalGain: return QStringLiteral("Global Gain");
        case ProcStepType::StepMigration: return QStringLiteral("Migration");
        default: return QStringLiteral("Unknown Step");
        }
    };

    GPRDataModel current = rawData;
    // 处理流水线只操作 traces，volumeData 不参与任何算法步骤
    // 提前清除可显著减少各步骤内部 GPRDataModel out = input 时的深拷贝开销
    current.volumeData.clear();
    // 严格按照steps数组存储顺序依次执行
    const int stepCount = pipeline.steps.size();
    for (int i = 0; i < pipeline.steps.size(); ++i)
    {
        if (canceled())
            break;

        const auto& unit = pipeline.steps[i];
        qDebug() << "[RadarProcessor] Step" << (i + 1) << "/" << stepCount << "starting:" << stepName(unit.type)
                 << "| traces=" << current.traces.size() << "samples=" << current.header.samplesPerTrace;

        switch (unit.type)
        {
        case ProcStepType::StepTimeZeroCut:
            current = cutTimeZero(current, unit.zeroCut);
            break;
        case ProcStepType::StepZeroDrift:
            current = removeZeroDrift(current, unit.zeroDrift);
            break;
        case ProcStepType::StepSphericalTVG:
            current = sphericalTvg(current, unit.tvg);
            break;
        case ProcStepType::StepBandpassFilter:
            current = bandpassFilter(current, unit.band);
            break;
        case ProcStepType::StepProfileSmooth:
            current = profileSmooth(current, unit.smooth);
            break;
        case ProcStepType::StepBackgroundRemove:
            current = backgroundRemoval(current, unit.bg);
            break;
        case ProcStepType::StepRemoveDC:
            current = removeDCShift(current, unit.dc);
            break;
        case ProcStepType::StepInnerAGC:
            current = traceInnerAgc(current, unit.innerAgc);
            break;
        case ProcStepType::StepTraceToTraceAGC:
            current = traceToTraceEqualization(current, unit.traceAgc);
            break;
        case ProcStepType::StepHilbertTransform:
            current = hilbertTransform(current, unit.hilbert);
            break;
        case ProcStepType::StepGlobalGain:
            current = applyGain(current, unit.gain);
            break;
        case ProcStepType::StepMigration:
            current = kirchhoffMigration(current, unit.migration);
            break;
        default:
            // 未知步骤类型跳过
            break;
        }

        qDebug() << "[RadarProcessor] Step" << (i + 1) << "finished:" << stepName(unit.type);

        if (context && context->reportProgress)
            context->reportProgress(i + 1, stepCount, stepName(unit.type));

        if (canceled())
            break;
    }
    return current;
}