feat(gpr3dv): 拷入 3DGPRViewer 数据生成链(geopro_gpr3dv 静态库)+冒烟

把参考实现的"多通道测线→GPRDataModel立方体→RadarProcessor处理"链
原样 vendored 进 external/gpr3dviewer/(算法零改动),生产管线A地基。

- 拷入: GPRDataModel.h/SurveyGeometry.h/RadarTypes.h/IprhParser/
  ImpulseMultiChannelConverter/Rd3Parser/RadarProcessor/PerformanceLogger
  + third_party/kissfft/*(逐字复制,未动算法体)。
- CMake: geopro_gpr3dv 静态库,链 Qt6::Core+Gui(QVector3D)+OpenMP+kissfft;
  enable_language(C)使kissfft .c入编;接进根构建。
- .gitignore: /external/* + 例外 !/external/gpr3dviewer/ 使 vendored 入库,
  qwt-src/vtk-install 仍忽略。
- 冒烟 tools/gpr3dv_smoke: 走原版 API loadImpulseMultiChannel→buildVolumeData
  →runPipeline(默认流水线)。线001冒烟: 通道数=14(读自数据)、立方体
  14x45305x821、处理前后平均绝对幅值 393.58→360.34(处理生效)。
- 全量构建通过,425/425 测试通过。
This commit is contained in:
gaozheng 2026-06-24 20:19:24 +08:00
parent d539fc1b73
commit 0efd84544c
24 changed files with 4833 additions and 1 deletions

6
.gitignore vendored
View File

@ -9,7 +9,11 @@ CMakeUserPresets.json
/vcpkg/
# ---- external source-built deps (方案②-修订: VTK 源码编到 install 前缀) ----
/external/
# 用 /external/* 而非 /external/ 忽略目录内容,使下方 vendored 子目录可被例外重新纳入
# (git 无法重新纳入被整体忽略目录内的文件)。
/external/*
# 例外vendored 3DGPRViewer 数据生成链(原样拷贝的算法源码,版权自有)需入库。
!/external/gpr3dviewer/
# ---- Visual Studio / IDE ----
.vs/

View File

@ -76,10 +76,17 @@ if(EXISTS "${CMAKE_SOURCE_DIR}/external/qwt-src/src")
include("${CMAKE_SOURCE_DIR}/cmake/qwt.cmake")
endif()
# vendored 3DGPRViewer geopro_gpr3dv
# .iprh/.iprb GPRDataModel RadarProcessor 线 A
add_subdirectory(external/gpr3dviewer)
add_subdirectory(src)
# POC-B headless CLIgpr_poc io_gpr/core/store/render
add_subdirectory(tools/gpr_poc)
# gpr3dv CLI vendored APIloadImpulseMultiChannel buildVolumeData runPipeline
add_subdirectory(tools/gpr3dv_smoke)
enable_testing()
add_subdirectory(tests)

49
external/gpr3dviewer/CMakeLists.txt vendored Normal file
View File

@ -0,0 +1,49 @@
# =====================================================================
# geopro_gpr3dv vendored 3DGPRViewer
#
# _Axx.iprh/.iprb
# IprhParser::loadImpulseMultiChannel GPRDataModeltraces
# GPRDataModel::buildVolumeData volumeData[ch][trace][sample]
# RadarProcessor::runPipeline GPRDataModel
#
# lanbingtech vendored
# Qt6::CoreQVector/QString/QFile/QVector3D/QJson+ OpenMP + kissfft
# UI/Sql/Network GPRWidget/mainwindow
# =====================================================================
# enable CXXkissfft C kiss_fft.c/kiss_fftr.c C
# .c CMake kiss_fftr/kiss_fftr_alloc
enable_language(C)
find_package(OpenMP)
add_library(geopro_gpr3dv STATIC
IprhParser.cpp
ImpulseMultiChannelConverter.cpp
Rd3Parser.cpp
RadarProcessor.cpp
PerformanceLogger.cpp
third_party/kissfft/kiss_fft.c
third_party/kissfft/kiss_fftr.c
)
# GPRDataModel.h + kissfftRadarProcessor.cpp #include "kiss_fftr.h"
target_include_directories(geopro_gpr3dv PUBLIC
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_SOURCE_DIR}/third_party/kissfft
)
# Qt6::Gui QVector3DGPRDataModel.h Core QVector/QString/QFile/QJson
target_link_libraries(geopro_gpr3dv PUBLIC Qt6::Core Qt6::Gui)
if(OpenMP_CXX_FOUND)
target_link_libraries(geopro_gpr3dv PUBLIC OpenMP::OpenMP_CXX)
endif()
target_compile_features(geopro_gpr3dv PRIVATE cxx_std_17)
# vendored Qt QObject/moc
set_target_properties(geopro_gpr3dv PROPERTIES AUTOMOC OFF AUTOUIC OFF AUTORCC OFF)
if(MSVC)
target_compile_options(geopro_gpr3dv PRIVATE /W0)
endif()

297
external/gpr3dviewer/GPRDataModel.h vendored Normal file
View File

@ -0,0 +1,297 @@
#ifndef GPRDATAMODEL_H
#define GPRDATAMODEL_H
#include <QVector>
#include <QString>
#include <QMap>
#include <QVector3D>
#include <QMetaType>
#include <utility>
#include "SurveyGeometry.h"
struct RadarTrace {
QVector<short> amplitudes; // 【重要】改为 short (16-bit),之前是 float
double startTime = 0.0;
double timeStep = 0.0;
// 三维位置信息
QVector3D position; // 该道的三维空间位置
int channelNumber = 0; // 通道号
RadarTrace() = default;
RadarTrace(const RadarTrace& other)
: amplitudes(QVector<short>(other.amplitudes.cbegin(), other.amplitudes.cend()))
, startTime(other.startTime)
, timeStep(other.timeStep)
, position(other.position)
, channelNumber(other.channelNumber)
{}
RadarTrace& operator=(const RadarTrace& other) {
if (this != &other) {
amplitudes = QVector<short>(other.amplitudes.cbegin(), other.amplitudes.cend());
startTime = other.startTime;
timeStep = other.timeStep;
position = other.position;
channelNumber = other.channelNumber;
}
return *this;
}
RadarTrace(RadarTrace&& other) noexcept
: amplitudes(std::move(other.amplitudes))
, startTime(other.startTime)
, timeStep(other.timeStep)
, position(other.position)
, channelNumber(other.channelNumber)
{}
RadarTrace& operator=(RadarTrace&& other) noexcept {
if (this != &other) {
amplitudes = std::move(other.amplitudes);
startTime = other.startTime;
timeStep = other.timeStep;
position = other.position;
channelNumber = other.channelNumber;
}
return *this;
}
};
class GPRDataModel {
public:
struct Header {
int numTraces = 0; // LAST TRACE
int samplesPerTrace = 0; // SAMPLES
double timeWindowNs = 0.0; // TIMEWINDOW (ns)
double distanceInc = 0.0; // DISTANCE INTERVAL (m)
double antennaFreq = 0.0; // FREQUENCY (MHz)
QString antennaType; // ANTENNAS
QString date; // DATE
QString timeStr; // TIME
// 三维相关参数
int numberOfChannels = 1; // NUMBER_OF_CH
QVector<float> chXOffsets; // CH_X_OFFSETS
QVector<float> chYOffsets; // CH_Y_OFFSETS
QString units; // UNITS
double startPosition = 0.0; // START POSITION
double stopPosition = 0.0; // STOP POSITION
double waveVelocity = 0.1; // 波速,用于深度计算
double timeIntervalNs = 0.0; // TIME INTERVAL (ns),用于计算时窗
// 原始映射,用于调试或扩展
QMap<QString, QString> rawParams;
};
Header header;
QVector<RadarTrace> traces;
// 三维数据相关属性
int tracesPerChannel = 0; // 每个通道的道数
int channels = 0; // 实际通道数
double totalDistance = 0.0; // 总距离
// 三维数据体 - 存储为[channel][trace][sample]格式
QVector<QVector<QVector<short>>> volumeData; // [channel][trace][sample]
bool isEmpty() const { return traces.isEmpty(); }
void clear() {
traces.clear();
header = Header{};
tracesPerChannel = 0;
channels = 0;
totalDistance = 0.0;
volumeData.clear();
}
// 获取指定通道的道数
int getTracesPerChannel() const {
if (header.numTraces <= 0) {
return 0;
}
if (header.numberOfChannels > 0) {
return header.numTraces / header.numberOfChannels;
}
return header.numTraces;
}
int getTraceIndex(int channel, int traceInChannel) const {
const int nChannels = header.numberOfChannels > 0 ? header.numberOfChannels : 1;
if (channel < 0 || channel >= nChannels ||
traceInChannel < 0 || traceInChannel >= getTracesPerChannel()) {
return -1;
}
return traceInChannel * nChannels + channel;
}
// 获取指定通道和道号的数据
const RadarTrace* getTrace(int channel, int traceInChannel) const {
const int index = getTraceIndex(channel, traceInChannel);
if (index >= 0 && index < traces.size()) {
return &traces[index];
}
return nullptr;
}
// 计算深度(米)
double calculateDepth(int sampleIndex) const {
if (header.samplesPerTrace > 0) {
// 使用默认波速0.1m/ns除以2是因为往返时间
double timePerSample = header.timeWindowNs / header.samplesPerTrace;
double timeNs = sampleIndex * timePerSample;
return timeNs * header.waveVelocity / 2.0; // 深度 = 时间 * 波速 / 2
}
return 0.0;
}
// 计算距离(米)
double calculateDistance(int traceIndex) const {
if (header.distanceInc > 0) {
return traceIndex * header.distanceInc;
}
return 0.0;
}
// 构建三维数据体
void buildVolumeData() {
if (traces.isEmpty()) return;
int nChannels = header.numberOfChannels > 0 ? header.numberOfChannels : 1;
int nTracesPerChannel = getTracesPerChannel();
int nSamples = header.samplesPerTrace;
volumeData.resize(nChannels);
for (int c = 0; c < nChannels; ++c) {
volumeData[c].resize(nTracesPerChannel);
for (int t = 0; t < nTracesPerChannel; ++t) {
volumeData[c][t].resize(nSamples);
}
}
// 填充数据
for (int i = 0; i < traces.size(); ++i) {
int channel = i % nChannels;
int traceInChannel = i / nChannels;
if (channel < volumeData.size() && traceInChannel < volumeData[channel].size()) {
for (int s = 0; s < traces[i].amplitudes.size() && s < nSamples; ++s) {
volumeData[channel][traceInChannel][s] = traces[i].amplitudes[s];
}
}
}
}
// 获取三维数据点
short getVolumeValue(int channel, int trace, int sample) const {
if (channel < 0 || channel >= volumeData.size() ||
trace < 0 || trace >= volumeData[channel].size() ||
sample < 0 || sample >= volumeData[channel][trace].size()) {
return 0;
}
return volumeData[channel][trace][sample];
}
// 获取全局最小最大值
void getGlobalMinMax(short& minVal, short& maxVal) const {
if (volumeData.isEmpty()) {
minVal = 0;
maxVal = 0;
return;
}
minVal = 32767;
maxVal = -32768;
for (const auto& channelData : volumeData) {
for (const auto& traceData : channelData) {
for (short val : traceData) {
if (val < minVal) minVal = val;
if (val > maxVal) maxVal = val;
}
}
}
}
// 获取通道数
int getChannelCount() const {
if (!volumeData.isEmpty()) return volumeData.size();
return qMax(1, header.numberOfChannels);
}
// 获取每通道的迹数
int getTraceCountPerChannel() const {
if (!volumeData.isEmpty() && !volumeData[0].isEmpty()) return volumeData[0].size();
int nCh = qMax(1, header.numberOfChannels);
return header.numTraces / nCh;
}
// 获取每迹的样本数
int getSampleCount() const {
if (!volumeData.isEmpty() && !volumeData[0].isEmpty()) return volumeData[0][0].size();
return header.samplesPerTrace;
}
GPRDataModel() = default;
GPRDataModel(const GPRDataModel&) = default;
GPRDataModel& operator=(const GPRDataModel&) = default;
GPRDataModel(GPRDataModel&&) = default;
GPRDataModel& operator=(GPRDataModel&&) = default;
};
struct TrajectoryEditPoint {
QVector3D rawPosition;
QVector3D editedPosition;
bool enabled = true;
bool manuallyEdited = false;
bool distanceOutlier = false;
bool speedOutlier = false;
bool angleOutlier = false;
};
struct TrajectoryEditSettings {
bool distanceFilterEnabled = true;
double maxSegmentDistanceM = 1.0;
bool speedFilterEnabled = false;
double maxSpeedMps = 5.0;
double traceIntervalSec = 0.05;
bool angleFilterEnabled = true;
double maxTurnAngleDeg = 60.0;
int interpolationMode = 0; // 0 linear, 1 spline
bool preserveManualEdits = true;
};
// 单条测线结构含GPS
struct SurveyLine {
QString name; // 测线编号,如 "_000", "_001"
QString displayName; // 显示名称,如 "碧新路_000"
QString radFilePath; // .rad 文件完整路径
GPRDataModel data; // 雷达数据
QVector<QVector3D> gpsPositions; // 每个trace的GPS坐标 (X,Y,Z)
QVector<TrajectoryEditPoint> trajectoryEditPoints; // 可编辑中心线轨迹点
TrajectoryEditSettings trajectoryEditSettings;
// 三维轨迹相关(新增)
SurveyGeometry geometry; // 天线几何参数
QVector<QVector<QVector3D>> channelTrajectories; // [channel][trace] 绝对坐标
bool hasGps() const { return !gpsPositions.isEmpty(); }
bool hasValidTrajectories() const { return !channelTrajectories.isEmpty() && !channelTrajectories[0].isEmpty(); }
bool isEmpty() const { return data.isEmpty(); }
void clear() {
name.clear();
displayName.clear();
radFilePath.clear();
data.clear();
gpsPositions.clear();
trajectoryEditPoints.clear();
trajectoryEditSettings = TrajectoryEditSettings{};
geometry = SurveyGeometry{};
channelTrajectories.clear();
}
SurveyLine() = default;
SurveyLine(const SurveyLine&) = default;
SurveyLine& operator=(const SurveyLine&) = default;
SurveyLine(SurveyLine&&) = default;
SurveyLine& operator=(SurveyLine&&) = default;
};
Q_DECLARE_METATYPE(SurveyLine)
#endif // GPRDATAMODEL_H

View File

@ -0,0 +1,406 @@
#include "ImpulseMultiChannelConverter.h"
#include "IprhParser.h"
#include <QDebug>
#include <QDir>
#include <QFile>
#include <QFileInfo>
#include <QLocale>
#include <QMap>
#include <QRegularExpression>
#include <QSharedPointer>
#include <QStringList>
#include <QTextStream>
#include <algorithm>
#include <limits>
bool ImpulseMultiChannelConverter::isMultiChannelImpulseHeader(const QString &headerPath,
QString *dirPath,
QString *baseName)
{
QFileInfo fi(headerPath);
QRegularExpression reMultiChannel(QStringLiteral("^(.*)_A(\\d+)$"),
QRegularExpression::CaseInsensitiveOption);
QRegularExpressionMatch match = reMultiChannel.match(fi.completeBaseName());
if (!match.hasMatch())
return false;
const QString surveyBase = match.captured(1);
QDir dir(fi.absolutePath());
const QStringList channels = dir.entryList(QStringList() << surveyBase + QStringLiteral("_A*.iprh"),
QDir::Files);
if (channels.size() <= 1)
return false;
if (dirPath)
*dirPath = fi.absolutePath();
if (baseName)
*baseName = surveyBase;
return true;
}
bool ImpulseMultiChannelConverter::buildPlan(const QString &dirPath,
const QString &baseName,
ConversionPlan &plan,
QString *errorMessage)
{
plan = ConversionPlan{};
plan.dirPath = dirPath;
plan.baseName = baseName;
QDir dir(dirPath);
dir.setFilter(QDir::Files | QDir::NoDotAndDotDot);
QRegularExpression reChannel(QStringLiteral("^%1_A(\\d+)\\.iprh$").arg(QRegularExpression::escape(baseName)),
QRegularExpression::CaseInsensitiveOption);
QMap<int, QString> channelHeaders;
for (const QFileInfo &fi : dir.entryInfoList()) {
QRegularExpressionMatch match = reChannel.match(fi.fileName());
if (!match.hasMatch())
continue;
const int chNum = match.captured(1).toInt();
const QString iprhPath = fi.absoluteFilePath();
QString iprbPath = iprhPath;
iprbPath.replace(QStringLiteral(".iprh"), QStringLiteral(".iprb"), Qt::CaseInsensitive);
if (QFile::exists(iprbPath)) {
channelHeaders.insert(chNum, iprhPath);
} else {
qDebug() << "Missing .iprb for" << iprhPath;
}
}
if (channelHeaders.isEmpty()) {
if (errorMessage)
*errorMessage = QStringLiteral("未找到 Impulse 多通道文件:%1").arg(baseName);
return false;
}
GPRDataModel::Header masterHeader;
if (!IprhParser::parseHeaderOnly(channelHeaders.first(), masterHeader)) {
if (errorMessage)
*errorMessage = QStringLiteral("无法解析主通道头文件:%1").arg(channelHeaders.first());
return false;
}
QVector<float> xOffsets;
QVector<float> yOffsets;
xOffsets.reserve(channelHeaders.size());
yOffsets.reserve(channelHeaders.size());
qint64 minTraceCount = std::numeric_limits<qint64>::max();
double antennaSeparation = 0.0;
for (auto it = channelHeaders.begin(); it != channelHeaders.end(); ++it) {
ChannelInfo info;
info.channelNumber = it.key();
info.iprhPath = it.value();
info.iprbPath = info.iprhPath;
info.iprbPath.replace(QStringLiteral(".iprh"), QStringLiteral(".iprb"), Qt::CaseInsensitive);
if (!IprhParser::parseHeaderOnly(info.iprhPath, info.header)) {
if (errorMessage)
*errorMessage = QStringLiteral("无法解析通道头文件:%1").arg(info.iprhPath);
return false;
}
if (info.header.samplesPerTrace != masterHeader.samplesPerTrace) {
if (errorMessage)
*errorMessage = QStringLiteral("多通道 SAMPLES 不一致:%1").arg(info.iprhPath);
return false;
}
if (!qFuzzyCompare(info.header.timeIntervalNs, masterHeader.timeIntervalNs) &&
info.header.timeIntervalNs > 0 && masterHeader.timeIntervalNs > 0) {
qDebug() << "Warning: Inconsistent TIME INTERVAL across channels" << info.iprhPath;
}
if (!qFuzzyCompare(info.header.distanceInc, masterHeader.distanceInc) &&
info.header.distanceInc > 0 && masterHeader.distanceInc > 0) {
qDebug() << "Warning: Inconsistent DISTANCE INTERVAL across channels" << info.iprhPath;
}
if (info.header.rawParams.contains(QStringLiteral("CH_X_OFFSET"))) {
info.xOffset = info.header.rawParams.value(QStringLiteral("CH_X_OFFSET")).toFloat();
} else if (info.header.rawParams.contains(QStringLiteral("CH_OFFSET_X"))) {
info.xOffset = info.header.rawParams.value(QStringLiteral("CH_OFFSET_X")).toFloat();
} else if (!info.header.chXOffsets.isEmpty()) {
info.xOffset = info.header.chXOffsets.first();
}
if (info.header.rawParams.contains(QStringLiteral("CH_Y_OFFSET"))) {
info.yOffset = info.header.rawParams.value(QStringLiteral("CH_Y_OFFSET")).toFloat();
} else if (info.header.rawParams.contains(QStringLiteral("CH_OFFSET_Y"))) {
info.yOffset = info.header.rawParams.value(QStringLiteral("CH_OFFSET_Y")).toFloat();
} else if (!info.header.chYOffsets.isEmpty()) {
info.yOffset = info.header.chYOffsets.first();
}
xOffsets.append(info.xOffset);
yOffsets.append(info.yOffset);
if (info.header.rawParams.contains(QStringLiteral("ANTENNA SEPARATION"))) {
antennaSeparation = info.header.rawParams.value(QStringLiteral("ANTENNA SEPARATION")).toDouble();
}
const qint64 traceBytes = static_cast<qint64>(info.header.samplesPerTrace) * sizeof(qint16);
if (traceBytes <= 0) {
if (errorMessage)
*errorMessage = QStringLiteral("通道采样数无效:%1").arg(info.iprhPath);
return false;
}
info.traceCount = QFileInfo(info.iprbPath).size() / traceBytes;
if (info.traceCount <= 0) {
if (errorMessage)
*errorMessage = QStringLiteral("通道无有效数据:%1").arg(info.iprbPath);
return false;
}
minTraceCount = std::min(minTraceCount, info.traceCount);
plan.channels.append(info);
}
if (plan.channels.isEmpty() || minTraceCount <= 0 || minTraceCount == std::numeric_limits<qint64>::max()) {
if (errorMessage)
*errorMessage = QStringLiteral("Impulse 多通道没有有效 trace%1").arg(baseName);
return false;
}
for (const ChannelInfo &info : plan.channels) {
if (info.traceCount != minTraceCount) {
qDebug() << "Warning: Inconsistent trace count across channels. Using minimum:"
<< minTraceCount << "channel has:" << info.traceCount << info.iprhPath;
}
}
bool allZeroOffsets = true;
for (float value : xOffsets) {
if (!qFuzzyIsNull(value)) {
allZeroOffsets = false;
break;
}
}
if (allZeroOffsets && antennaSeparation > 1e-6 && plan.channels.size() > 1) {
const double totalWidth = (plan.channels.size() - 1) * antennaSeparation;
const double startX = -totalWidth / 2.0;
for (int i = 0; i < xOffsets.size(); ++i) {
xOffsets[i] = static_cast<float>(startX + i * antennaSeparation);
plan.channels[i].xOffset = xOffsets[i];
}
qDebug() << "Computed symmetric X offsets from antenna separation:" << antennaSeparation;
}
plan.channelCount = plan.channels.size();
plan.tracesPerChannel = static_cast<int>(minTraceCount);
plan.samplesPerTrace = masterHeader.samplesPerTrace;
plan.traceByteSize = plan.samplesPerTrace * static_cast<qint64>(sizeof(qint16));
plan.totalOutputTraces = static_cast<qint64>(plan.tracesPerChannel) * plan.channelCount;
plan.outputHeader = masterHeader;
plan.outputHeader.numberOfChannels = plan.channelCount;
plan.outputHeader.numTraces = static_cast<int>(plan.totalOutputTraces);
plan.outputHeader.chXOffsets = xOffsets;
plan.outputHeader.chYOffsets = yOffsets;
if (plan.outputHeader.timeWindowNs <= 0.0 && plan.outputHeader.timeIntervalNs > 0.0) {
plan.outputHeader.timeWindowNs = plan.outputHeader.samplesPerTrace * plan.outputHeader.timeIntervalNs;
}
const QString basePath = dir.absoluteFilePath(baseName + QStringLiteral("_mala_converted"));
plan.outputRadPath = basePath + QStringLiteral(".rad");
plan.outputRd3Path = basePath + QStringLiteral(".rd3");
return true;
}
bool ImpulseMultiChannelConverter::convertStreaming(const ConversionPlan &plan,
const Options &options,
QString *radFilePath,
QString *errorMessage,
CancelFn cancel,
ProgressFn progress)
{
if (plan.channelCount <= 0 || plan.tracesPerChannel <= 0 || plan.traceByteSize <= 0) {
if (errorMessage)
*errorMessage = QStringLiteral("Impulse 转换计划无效");
return false;
}
if (options.reuseExistingIfValid && convertedFilesAreValid(plan)) {
if (radFilePath)
*radFilePath = plan.outputRadPath;
return true;
}
if (!options.overwriteExisting && (QFile::exists(plan.outputRadPath) || QFile::exists(plan.outputRd3Path))) {
if (errorMessage)
*errorMessage = QStringLiteral("转换输出文件已存在:%1").arg(plan.outputRadPath);
return false;
}
const QString tmpRd3Path = plan.outputRd3Path + QStringLiteral(".tmp");
const QString tmpRadPath = plan.outputRadPath + QStringLiteral(".tmp");
QFile::remove(tmpRd3Path);
QFile::remove(tmpRadPath);
QVector<QSharedPointer<QFile>> inputFiles;
inputFiles.reserve(plan.channels.size());
for (const ChannelInfo &channel : plan.channels) {
auto file = QSharedPointer<QFile>::create(channel.iprbPath);
if (!file->open(QIODevice::ReadOnly)) {
if (errorMessage)
*errorMessage = QStringLiteral("无法读取 Impulse 通道数据:%1").arg(channel.iprbPath);
return false;
}
inputFiles.append(file);
}
QFile rd3File(tmpRd3Path);
if (!rd3File.open(QIODevice::WriteOnly | QIODevice::Truncate)) {
if (errorMessage)
*errorMessage = QStringLiteral("无法写入转换后的 RD3 文件:%1").arg(tmpRd3Path);
return false;
}
const qint64 bytesPerTracePosition = plan.traceByteSize * plan.channelCount;
const qint64 chunkBudget = qMax<qint64>(bytesPerTracePosition, options.maxChunkBytes);
const int chunkTraces = qMax<qint64>(1, chunkBudget / bytesPerTracePosition);
QVector<QByteArray> channelBuffers(plan.channelCount);
qint64 tracesDone = 0;
while (tracesDone < plan.tracesPerChannel) {
if (cancel && cancel()) {
rd3File.close();
QFile::remove(tmpRd3Path);
QFile::remove(tmpRadPath);
if (errorMessage)
*errorMessage = QStringLiteral("Impulse 多通道转换已取消");
return false;
}
const int currentChunk = qMin<qint64>(chunkTraces, plan.tracesPerChannel - tracesDone);
const qint64 expectedChannelBytes = currentChunk * plan.traceByteSize;
for (int ch = 0; ch < plan.channelCount; ++ch) {
channelBuffers[ch] = inputFiles[ch]->read(expectedChannelBytes);
if (channelBuffers[ch].size() != expectedChannelBytes) {
rd3File.close();
QFile::remove(tmpRd3Path);
if (errorMessage) {
*errorMessage = QStringLiteral("读取通道数据不完整:%1").arg(plan.channels[ch].iprbPath);
}
return false;
}
}
for (int localTrace = 0; localTrace < currentChunk; ++localTrace) {
const qint64 offset = localTrace * plan.traceByteSize;
for (int ch = 0; ch < plan.channelCount; ++ch) {
const qint64 written = rd3File.write(channelBuffers[ch].constData() + offset, plan.traceByteSize);
if (written != plan.traceByteSize) {
rd3File.close();
QFile::remove(tmpRd3Path);
if (errorMessage)
*errorMessage = QStringLiteral("写入转换后的 RD3 文件失败:%1").arg(tmpRd3Path);
return false;
}
}
}
tracesDone += currentChunk;
if (progress) {
progress(Progress{tracesDone, plan.tracesPerChannel,
QStringLiteral("正在转换 Impulse 多通道 %1/%2")
.arg(tracesDone)
.arg(plan.tracesPerChannel)});
}
}
rd3File.close();
if (rd3File.error() != QFile::NoError) {
QFile::remove(tmpRd3Path);
if (errorMessage)
*errorMessage = QStringLiteral("保存转换后的 RD3 文件失败:%1").arg(rd3File.errorString());
return false;
}
ConversionPlan tmpPlan = plan;
tmpPlan.outputRadPath = tmpRadPath;
if (!writeRadHeader(tmpRadPath, tmpPlan, errorMessage)) {
QFile::remove(tmpRd3Path);
QFile::remove(tmpRadPath);
return false;
}
QFile::remove(plan.outputRd3Path);
QFile::remove(plan.outputRadPath);
if (!QFile::rename(tmpRd3Path, plan.outputRd3Path)) {
QFile::remove(tmpRd3Path);
QFile::remove(tmpRadPath);
if (errorMessage)
*errorMessage = QStringLiteral("无法替换转换后的 RD3 文件:%1").arg(plan.outputRd3Path);
return false;
}
if (!QFile::rename(tmpRadPath, plan.outputRadPath)) {
QFile::remove(plan.outputRd3Path);
QFile::remove(tmpRadPath);
if (errorMessage)
*errorMessage = QStringLiteral("无法替换转换后的 RAD 文件:%1").arg(plan.outputRadPath);
return false;
}
if (radFilePath)
*radFilePath = plan.outputRadPath;
qDebug() << "Impulse multi-channel streaming conversion OK:" << plan.outputRadPath << plan.outputRd3Path;
return true;
}
bool ImpulseMultiChannelConverter::writeRadHeader(const QString &radFilePath,
const ConversionPlan &plan,
QString *errorMessage)
{
QFile radFile(radFilePath);
if (!radFile.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate)) {
if (errorMessage)
*errorMessage = QStringLiteral("无法写入转换后的 RAD 文件:%1").arg(radFilePath);
return false;
}
QTextStream rad(&radFile);
const auto &header = plan.outputHeader;
rad << "# Converted from Impulse multi-channel survey: " << plan.baseName << '\n';
rad << "SAMPLES: " << header.samplesPerTrace << '\n';
rad << "LAST TRACE: " << plan.totalOutputTraces << '\n';
rad << "TIMEWINDOW: " << QLocale::c().toString(header.timeWindowNs, 'g', 12) << '\n';
rad << "TIME INTERVAL: " << QLocale::c().toString(header.timeIntervalNs, 'g', 12) << '\n';
rad << "DISTANCE INTERVAL: " << QLocale::c().toString(header.distanceInc, 'g', 12) << '\n';
rad << "ANTENNAS: " << QLocale::c().toString(header.antennaFreq, 'g', 12) << '\n';
if (!header.antennaType.isEmpty()) rad << "ANTENNA: " << header.antennaType << '\n';
if (!header.date.isEmpty()) rad << "DATE: " << header.date << '\n';
if (!header.timeStr.isEmpty()) rad << "TIME: " << header.timeStr << '\n';
rad << "NUMBER_OF_CH: " << qMax(1, header.numberOfChannels) << '\n';
rad << "CH_X_OFFSETS: " << formatFloatList(header.chXOffsets) << '\n';
rad << "CH_Y_OFFSETS: " << formatFloatList(header.chYOffsets) << '\n';
if (!header.units.isEmpty()) rad << "UNITS: " << header.units << '\n';
rad << "START POSITION: " << QLocale::c().toString(header.startPosition, 'g', 12) << '\n';
const double stopPosition = header.stopPosition > header.startPosition
? header.stopPosition
: header.startPosition + plan.tracesPerChannel * header.distanceInc;
rad << "STOP POSITION: " << QLocale::c().toString(stopPosition, 'g', 12) << '\n';
radFile.close();
return true;
}
bool ImpulseMultiChannelConverter::convertedFilesAreValid(const ConversionPlan &plan)
{
QFileInfo radInfo(plan.outputRadPath);
QFileInfo rd3Info(plan.outputRd3Path);
if (!radInfo.exists() || !rd3Info.exists())
return false;
const qint64 expectedBytes = plan.totalOutputTraces * plan.traceByteSize;
return rd3Info.size() == expectedBytes;
}
QString ImpulseMultiChannelConverter::formatFloatList(const QVector<float> &values)
{
QStringList parts;
parts.reserve(values.size());
for (float value : values) {
parts.append(QLocale::c().toString(value, 'g', 10));
}
return parts.join(' ');
}

View File

@ -0,0 +1,75 @@
#ifndef IMPULSEMULTICHANNELCONVERTER_H
#define IMPULSEMULTICHANNELCONVERTER_H
#include "GPRDataModel.h"
#include <QVector>
#include <QString>
#include <functional>
class ImpulseMultiChannelConverter {
public:
struct ChannelInfo {
int channelNumber = 0;
QString iprhPath;
QString iprbPath;
GPRDataModel::Header header;
qint64 traceCount = 0;
float xOffset = 0.0f;
float yOffset = 0.0f;
};
struct ConversionPlan {
QString dirPath;
QString baseName;
QVector<ChannelInfo> channels;
GPRDataModel::Header outputHeader;
int channelCount = 0;
int tracesPerChannel = 0;
qint64 totalOutputTraces = 0;
qint64 samplesPerTrace = 0;
qint64 traceByteSize = 0;
QString outputRadPath;
QString outputRd3Path;
};
struct Options {
qint64 maxChunkBytes = 32 * 1024 * 1024;
bool overwriteExisting = true;
bool reuseExistingIfValid = false;
};
struct Progress {
qint64 tracesDone = 0;
qint64 tracesTotal = 0;
QString message;
};
using CancelFn = std::function<bool()>;
using ProgressFn = std::function<void(const Progress&)>;
static bool isMultiChannelImpulseHeader(const QString &headerPath,
QString *dirPath = nullptr,
QString *baseName = nullptr);
static bool buildPlan(const QString &dirPath,
const QString &baseName,
ConversionPlan &plan,
QString *errorMessage = nullptr);
static bool convertStreaming(const ConversionPlan &plan,
const Options &options,
QString *radFilePath = nullptr,
QString *errorMessage = nullptr,
CancelFn cancel = {},
ProgressFn progress = {});
private:
static bool writeRadHeader(const QString &radFilePath,
const ConversionPlan &plan,
QString *errorMessage);
static bool convertedFilesAreValid(const ConversionPlan &plan);
static QString formatFloatList(const QVector<float> &values);
};
#endif // IMPULSEMULTICHANNELCONVERTER_H

625
external/gpr3dviewer/IprhParser.cpp vendored Normal file
View File

@ -0,0 +1,625 @@
#include "IprhParser.h"
#include "PerformanceLogger.h"
#include "ImpulseMultiChannelConverter.h"
#include <QFile>
#include <QTextStream>
#include <QDataStream>
#include <QDebug>
#include <QDir>
#include <QFileInfo>
#include <QRegularExpression>
#include <QStringList>
#include <QVector3D>
#include <QLocale>
#include <algorithm>
#include <cstring>
/**
* @brief Impulse IPRH iprh头文件 + iprb纯二进制波形
* @param iprhFilePath .iprh
* @param model GPR
* @return true false
* @note .iprh /线.iprb short16
*/
bool IprhParser::loadFromIprh(const QString &iprhFilePath, GPRDataModel &model) {
SCOPED_PERF_TIMER("Parser", "IprhParser::loadFromIprh");
model.clear();
model.header = GPRDataModel::Header{};
// 第一步:解析 iprh 文本头
if (!parseIprhHeader(iprhFilePath, model.header)) {
qDebug() << "Error: Failed to parse .iprh header file:" << iprhFilePath;
return false;
}
// 自动匹配同目录同名二进制数据文件 .iprb
QFileInfo iprhInfo(iprhFilePath);
QString iprbPath = iprhInfo.absolutePath() + "/" + iprhInfo.completeBaseName() + ".iprb";
if (!QFile::exists(iprbPath)) {
qDebug() << "Error: Matching .iprb binary file not found at:" << iprbPath;
return false;
}
// 如果头文件没有 LAST TRACE从二进制文件大小推算
if (model.header.numTraces <= 0) {
QFile binaryFile(iprbPath);
if (binaryFile.open(QIODevice::ReadOnly)) {
qint64 fileSize = binaryFile.size();
qint64 traceBytes = static_cast<qint64>(model.header.samplesPerTrace) * sizeof(short);
if (traceBytes > 0) {
model.header.numTraces = static_cast<int>(fileSize / traceBytes);
qDebug() << "Inferred numTraces from binary size:" << model.header.numTraces;
}
binaryFile.close();
}
}
if (model.header.timeWindowNs <= 0.0 && model.header.timeIntervalNs > 0.0) {
model.header.timeWindowNs = model.header.samplesPerTrace * model.header.timeIntervalNs;
qDebug() << "Inferred timeWindowNs from samples * timeInterval:" << model.header.timeWindowNs;
}
if (model.header.numTraces <= 0) {
qDebug() << "Error: Unable to determine numTraces from header or binary file";
return false;
}
// 第二步:读取纯二进制波形
return loadIprbBinary(iprbPath, model);
}
/**
* @brief IPRH
* @param iprhFilePath iprh
* @param header
* @return true
*/
bool IprhParser::parseHeaderOnly(const QString &iprhFilePath, GPRDataModel::Header &header)
{
return parseIprhHeader(iprhFilePath, header);
}
bool IprhParser::parseIprhHeader(const QString &iprhFilePath, GPRDataModel::Header &header) {
SCOPED_PERF_TIMER("Parser", "IprhParser::parseIprhHeader");
QFile file(iprhFilePath);
if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) {
qDebug() << "Error: Cannot open iprh file for read";
return false;
}
QTextStream in(&file);
while (!in.atEnd()) {
QString line = in.readLine().trimmed();
if (line.isEmpty()) continue;
int sepPos = line.indexOf(':');
if (sepPos == -1) continue;
QString key = line.left(sepPos).trimmed();
QString value = line.mid(sepPos + 1).trimmed();
header.rawParams[key] = value;
if (key == "SAMPLES") {
header.samplesPerTrace = extractInt(value);
} else if (key == "LAST TRACE") {
header.numTraces = extractInt(value);
} else if (key == "TIMEWINDOW") {
header.timeWindowNs = extractDouble(value);
} else if (key == "DISTANCE INTERVAL") {
header.distanceInc = extractDouble(value);
} else if (key == "TIME INTERVAL") {
header.timeIntervalNs = extractDouble(value);
} else if (key == "ANTENNAS") {
header.antennaFreq = extractDouble(value);
} else if (key == "ANTENNA") {
header.antennaType = value;
} else if (key == "DATE") {
header.date = value;
} else if (key == "START TIME" || key == "TIME") {
header.timeStr = value;
} else if (key == "CHANNELS" || key == "NUMBER_OF_CH") {
header.numberOfChannels = extractInt(value);
} else if (key == "CH_X_OFFSET" || key == "CH_OFFSET_X") {
// 单通道偏移量Impulse 单通道文件)
if (!value.isEmpty()) {
header.chXOffsets.append(value.toFloat());
}
} else if (key == "CH_Y_OFFSET" || key == "CH_OFFSET_Y") {
if (!value.isEmpty()) {
header.chYOffsets.append(value.toFloat());
}
} else if (key == "CH_X_OFFSETS") {
// 兼容 Mala Mira 风格多值空格分隔
QStringList offsets = value.split(' ', Qt::SkipEmptyParts);
for (const QString &offset : offsets) {
header.chXOffsets.append(offset.toFloat());
}
} else if (key == "CH_Y_OFFSETS") {
QStringList offsets = value.split(' ', Qt::SkipEmptyParts);
for (const QString &offset : offsets) {
header.chYOffsets.append(offset.toFloat());
}
} else if (key == "UNITS") {
header.units = value;
} else if (key == "START POSITION") {
header.startPosition = extractDouble(value);
} else if (key == "STOP POSITION") {
header.stopPosition = extractDouble(value);
}
}
file.close();
if (header.samplesPerTrace <= 0) {
qDebug() << "Error: Invalid SAMPLES value in iprh file";
return false;
}
header.waveVelocity = 0.1;
qDebug() << "==== IPRH Header Parse Complete ===="
<< "\nSamples per trace:" << header.samplesPerTrace
<< "\nTotal traces:" << header.numTraces
<< "\nTime window(ns):" << header.timeWindowNs
<< "\nChannel count:" << header.numberOfChannels
<< "\nX offset array size:" << header.chXOffsets.size()
<< "\nY offset array size:" << header.chYOffsets.size();
return true;
}
/**
* @brief IPRB short16
* @param iprbFilePath iprb
* @param model traces
* @return true
*/
bool IprhParser::loadIprbBinary(const QString &iprbFilePath, GPRDataModel &model) {
SCOPED_PERF_TIMER("Parser", "IprhParser::loadIprbBinary");
QFile file(iprbFilePath);
if (!file.open(QIODevice::ReadOnly)) {
qDebug() << "Error: Open iprb binary failed:" << iprbFilePath;
return false;
}
const int samplesPerTrace = model.header.samplesPerTrace;
const int totalTraceCount = model.header.numTraces;
const qint64 dataStartOffset = 0;
file.seek(dataStartOffset);
const qint64 singleTraceByteSize = samplesPerTrace * sizeof(short);
const qint64 fullExpectedBytes = totalTraceCount * singleTraceByteSize;
const qint64 realFileBytes = file.size();
if (realFileBytes < fullExpectedBytes) {
qDebug() << "Warning: IPRB file size smaller than theoretical data size!"
<< "Expected:" << fullExpectedBytes << "Actual:" << realFileBytes;
}
model.traces.reserve(totalTraceCount);
QByteArray traceBuffer;
traceBuffer.resize(singleTraceByteSize);
const int channelCnt = model.header.numberOfChannels > 0 ? model.header.numberOfChannels : 1;
model.channels = channelCnt;
model.tracesPerChannel = totalTraceCount / channelCnt;
if (model.header.distanceInc > 1e-6) {
model.totalDistance = static_cast<float>(model.tracesPerChannel * model.header.distanceInc);
} else {
model.totalDistance = static_cast<float>(model.header.stopPosition - model.header.startPosition);
}
for (int traceGlobalIdx = 0; traceGlobalIdx < totalTraceCount; ++traceGlobalIdx) {
if (file.atEnd()) {
qDebug() << "Warning: File ended early at global trace index" << traceGlobalIdx;
break;
}
qint64 readBytes = file.read(traceBuffer.data(), traceBuffer.size());
if (readBytes != traceBuffer.size()) {
qDebug() << "Warning: Trace" << traceGlobalIdx << "incomplete byte read";
break;
}
RadarTrace oneTrace;
oneTrace.amplitudes.resize(samplesPerTrace);
const short* rawShortBuf = reinterpret_cast<const short*>(traceBuffer.constData());
for (int s = 0; s < samplesPerTrace; s++) {
oneTrace.amplitudes[s] = rawShortBuf[s];
}
int chNo = traceGlobalIdx % channelCnt;
int traceInChIdx = traceGlobalIdx / channelCnt;
oneTrace.channelNumber = chNo;
float xOff = 0.0f;
float yOff = 0.0f;
if (chNo < model.header.chXOffsets.size()) xOff = model.header.chXOffsets[chNo];
if (chNo < model.header.chYOffsets.size()) yOff = model.header.chYOffsets[chNo];
float lineDist = static_cast<float>(model.header.startPosition + traceInChIdx * model.header.distanceInc);
oneTrace.position = QVector3D(xOff, lineDist - yOff, 0.0f);
model.traces.append(std::move(oneTrace));
}
file.close();
if (model.traces.isEmpty()) {
qDebug() << "Error: No valid traces loaded from iprb";
return false;
}
qDebug() << "IPRB Binary Load OK, total valid traces:" << model.traces.size();
return true;
}
/**
* @brief Impulse .iprh + .iprb
* @param dirPath
* @param baseName 线 "明星路_001"
* @param model GPR
* @return true
*
* baseName_A01.iprh / .iprb ... baseName_A14.iprh / .iprb
* traces Mala Mira trace0=ch0_pos0, trace1=ch1_pos0, ...
*/
bool IprhParser::convertImpulseMultiChannelToMala(const QString &dirPath,
const QString &baseName,
QString *radFilePath,
QString *errorMessage)
{
ImpulseMultiChannelConverter::ConversionPlan plan;
if (!ImpulseMultiChannelConverter::buildPlan(dirPath, baseName, plan, errorMessage)) {
if (errorMessage && errorMessage->isEmpty()) {
*errorMessage = QStringLiteral("多通道 Impulse 数据合并失败:%1").arg(baseName);
}
return false;
}
ImpulseMultiChannelConverter::Options options;
options.overwriteExisting = true;
options.reuseExistingIfValid = false;
if (!ImpulseMultiChannelConverter::convertStreaming(plan, options, radFilePath, errorMessage)) {
return false;
}
qDebug() << "Impulse multi-channel converted to Mala Mira files:" << plan.outputRadPath << plan.outputRd3Path;
return true;
}
/**
* @brief Impulse .iprh + .iprb
* @param dirPath
* @param baseName 线 "明星路_001"
* @param model GPR
* @return true
*
* baseName_A01.iprh / .iprb ... baseName_A14.iprh / .iprb
* traces Mala Mira trace0=ch0_pos0, trace1=ch1_pos0, ...
*/
bool IprhParser::loadImpulseMultiChannel(const QString &dirPath,
const QString &baseName,
GPRDataModel &model)
{
SCOPED_PERF_TIMER("Parser", "IprhParser::loadImpulseMultiChannel");
model.clear();
model.header = GPRDataModel::Header{};
// 1. 发现所有通道文件并排序
QDir dir(dirPath);
dir.setFilter(QDir::Files | QDir::NoDotAndDotDot);
QRegularExpression reChannel(QStringLiteral("^%1_A(\\d+)\\.iprh$").arg(QRegularExpression::escape(baseName)));
QMap<int, QString> channelHeaders; // channelNum -> iprh path
for (const QFileInfo &fi : dir.entryInfoList()) {
QRegularExpressionMatch match = reChannel.match(fi.fileName());
if (match.hasMatch()) {
int chNum = match.captured(1).toInt();
QString iprhPath = fi.absoluteFilePath();
QString iprbPath = iprhPath;
iprbPath.replace(".iprh", ".iprb");
if (QFile::exists(iprbPath)) {
channelHeaders.insert(chNum, iprhPath);
} else {
qDebug() << "Missing .iprb for" << iprhPath;
}
}
}
if (channelHeaders.isEmpty()) {
qDebug() << "No multi-channel .iprh/.iprb found for baseName:" << baseName;
return false;
}
const int channelCount = channelHeaders.size();
qDebug() << "Impulse multi-channel: found" << channelCount << "channels for" << baseName;
// 2. 解析第一个通道的头文件作为 master
auto it = channelHeaders.begin();
GPRDataModel::Header masterHeader;
if (!parseIprhHeader(it.value(), masterHeader)) {
qDebug() << "Failed to parse master header:" << it.value();
return false;
}
QVector<float> xOffsets;
QVector<float> yOffsets;
xOffsets.reserve(channelCount);
yOffsets.reserve(channelCount);
// 3. 验证各通道一致性并收集偏移量
QVector<int> channelTracesPerChannel;
channelTracesPerChannel.reserve(channelCount);
double antennaSeparation = 0.0;
for (auto cit = channelHeaders.begin(); cit != channelHeaders.end(); ++cit) {
int chNum = cit.key();
QString iprhPath = cit.value();
GPRDataModel::Header chHeader;
if (!parseIprhHeader(iprhPath, chHeader)) {
qDebug() << "Failed to parse channel header:" << iprhPath;
return false;
}
// 验证关键参数一致性
if (chHeader.samplesPerTrace != masterHeader.samplesPerTrace) {
qDebug() << "Inconsistent SAMPLES across channels";
return false;
}
if (qFuzzyCompare(chHeader.timeIntervalNs, masterHeader.timeIntervalNs) == false &&
chHeader.timeIntervalNs > 0 && masterHeader.timeIntervalNs > 0) {
qDebug() << "Warning: Inconsistent TIME INTERVAL across channels";
}
if (qFuzzyCompare(chHeader.distanceInc, masterHeader.distanceInc) == false &&
chHeader.distanceInc > 0 && masterHeader.distanceInc > 0) {
qDebug() << "Warning: Inconsistent DISTANCE INTERVAL across channels";
}
// 收集单通道偏移量
float xOff = 0.0f, yOff = 0.0f;
if (chHeader.rawParams.contains("CH_X_OFFSET")) {
xOff = chHeader.rawParams.value("CH_X_OFFSET").toFloat();
} else if (chHeader.rawParams.contains("CH_OFFSET_X")) {
xOff = chHeader.rawParams.value("CH_OFFSET_X").toFloat();
} else if (!chHeader.chXOffsets.isEmpty()) {
xOff = chHeader.chXOffsets.first();
}
if (chHeader.rawParams.contains("CH_Y_OFFSET")) {
yOff = chHeader.rawParams.value("CH_Y_OFFSET").toFloat();
} else if (chHeader.rawParams.contains("CH_OFFSET_Y")) {
yOff = chHeader.rawParams.value("CH_OFFSET_Y").toFloat();
} else if (!chHeader.chYOffsets.isEmpty()) {
yOff = chHeader.chYOffsets.first();
}
xOffsets.append(xOff);
yOffsets.append(yOff);
if (chHeader.rawParams.contains("ANTENNA SEPARATION")) {
antennaSeparation = chHeader.rawParams.value("ANTENNA SEPARATION").toDouble();
}
// 从 .iprb 大小计算该通道的道数
QString iprbPath = iprhPath;
iprbPath.replace(".iprh", ".iprb");
QFile iprbFile(iprbPath);
int tracesInThisChannel = 0;
if (iprbFile.open(QIODevice::ReadOnly)) {
qint64 fileSize = iprbFile.size();
qint64 traceBytes = static_cast<qint64>(chHeader.samplesPerTrace) * sizeof(short);
if (traceBytes > 0) {
tracesInThisChannel = static_cast<int>(fileSize / traceBytes);
}
iprbFile.close();
}
channelTracesPerChannel.append(tracesInThisChannel);
}
// 4. 以最小道数为准截断读取,容忍各通道微小差异
int tracesPerChannel = channelTracesPerChannel.isEmpty() ? 0 : *std::min_element(channelTracesPerChannel.begin(), channelTracesPerChannel.end());
for (int tc : channelTracesPerChannel) {
if (tc != tracesPerChannel) {
qDebug() << "Warning: Inconsistent trace count across channels. Using minimum:" << tracesPerChannel << "channel has:" << tc;
}
}
if (tracesPerChannel <= 0) {
qDebug() << "Error: No valid traces in any channel";
return false;
}
// 5. 如果偏移量全部为零,基于 ANTENNA SEPARATION 计算对称分布
bool allZeroOffsets = true;
for (float v : xOffsets) { if (qFuzzyIsNull(v) == false) { allZeroOffsets = false; break; } }
if (allZeroOffsets && antennaSeparation > 1e-6 && channelCount > 1) {
double totalWidth = (channelCount - 1) * antennaSeparation;
double startX = -totalWidth / 2.0;
for (int i = 0; i < channelCount; ++i) {
xOffsets[i] = static_cast<float>(startX + i * antennaSeparation);
}
qDebug() << "Computed symmetric X offsets from antenna separation:" << antennaSeparation;
}
// 6. 组装合并后的 Header
model.header = masterHeader;
model.header.numberOfChannels = channelCount;
model.header.numTraces = tracesPerChannel * channelCount;
model.header.chXOffsets = xOffsets;
model.header.chYOffsets = yOffsets;
if (model.header.timeWindowNs <= 0.0 && model.header.timeIntervalNs > 0.0) {
model.header.timeWindowNs = model.header.samplesPerTrace * model.header.timeIntervalNs;
}
// 7. 预分配 traces
model.traces.reserve(model.header.numTraces);
// 8. 读取每个通道的 .iprb 到临时数组
QVector<QVector<RadarTrace>> channelTraceArrays;
channelTraceArrays.resize(channelCount);
int chIdx = 0;
for (auto cit = channelHeaders.begin(); cit != channelHeaders.end(); ++cit, ++chIdx) {
QString iprhPath = cit.value();
QString iprbPath = iprhPath;
iprbPath.replace(".iprh", ".iprb");
QFile file(iprbPath);
if (!file.open(QIODevice::ReadOnly)) {
qDebug() << "Error: Cannot open" << iprbPath;
return false;
}
const int samplesPerTrace = model.header.samplesPerTrace;
const qint64 singleTraceByteSize = samplesPerTrace * sizeof(short);
QByteArray traceBuffer;
traceBuffer.resize(singleTraceByteSize);
channelTraceArrays[chIdx].reserve(tracesPerChannel);
for (int t = 0; t < tracesPerChannel; ++t) {
if (file.atEnd()) break;
qint64 readBytes = file.read(traceBuffer.data(), traceBuffer.size());
if (readBytes != traceBuffer.size()) break;
RadarTrace oneTrace;
oneTrace.amplitudes.resize(samplesPerTrace);
const short* rawShortBuf = reinterpret_cast<const short*>(traceBuffer.constData());
for (int s = 0; s < samplesPerTrace; s++) {
oneTrace.amplitudes[s] = rawShortBuf[s];
}
channelTraceArrays[chIdx].append(std::move(oneTrace));
}
file.close();
if (channelTraceArrays[chIdx].size() != tracesPerChannel) {
qDebug() << "Warning: Channel" << cit.key() << "has fewer traces than expected";
}
}
// 9. 按 Mala Mira 格式交错合并
model.channels = channelCount;
model.tracesPerChannel = tracesPerChannel;
for (int pos = 0; pos < tracesPerChannel; ++pos) {
for (int ch = 0; ch < channelCount; ++ch) {
if (pos < channelTraceArrays[ch].size()) {
RadarTrace trace = std::move(channelTraceArrays[ch][pos]);
trace.channelNumber = ch;
float xOff = (ch < xOffsets.size()) ? xOffsets[ch] : 0.0f;
float yOff = (ch < yOffsets.size()) ? yOffsets[ch] : 0.0f;
float lineDist = static_cast<float>(model.header.startPosition + pos * model.header.distanceInc);
trace.position = QVector3D(xOff, lineDist - yOff, 0.0f);
model.traces.append(std::move(trace));
}
}
}
if (model.header.distanceInc > 1e-6) {
model.totalDistance = static_cast<float>(model.tracesPerChannel * model.header.distanceInc);
} else {
model.totalDistance = static_cast<float>(model.header.stopPosition - model.header.startPosition);
}
qDebug() << "Impulse multi-channel load OK. Total traces:" << model.traces.size()
<< "Channels:" << channelCount << "Traces/Channel:" << tracesPerChannel;
return !model.traces.isEmpty();
}
bool IprhParser::writeMalaFiles(const QString &radFilePath,
const QString &rd3FilePath,
const GPRDataModel &model,
const QString &sourceBaseName,
QString *errorMessage)
{
QFile rd3File(rd3FilePath);
if (!rd3File.open(QIODevice::WriteOnly | QIODevice::Truncate)) {
if (errorMessage) {
*errorMessage = QStringLiteral("无法写入转换后的 RD3 文件:%1").arg(rd3FilePath);
}
return false;
}
QDataStream out(&rd3File);
out.setByteOrder(QDataStream::LittleEndian);
const int samplesPerTrace = model.header.samplesPerTrace;
for (const RadarTrace &trace : model.traces) {
for (int s = 0; s < samplesPerTrace; ++s) {
const qint16 sample = static_cast<qint16>(s < trace.amplitudes.size() ? trace.amplitudes[s] : 0);
out << sample;
}
}
rd3File.close();
QFile radFile(radFilePath);
if (!radFile.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate)) {
if (errorMessage) {
*errorMessage = QStringLiteral("无法写入转换后的 RAD 文件:%1").arg(radFilePath);
}
QFile::remove(rd3FilePath);
return false;
}
QTextStream rad(&radFile);
rad << "# Converted from Impulse multi-channel survey: " << sourceBaseName << '\n';
rad << "SAMPLES: " << model.header.samplesPerTrace << '\n';
rad << "LAST TRACE: " << model.traces.size() << '\n';
rad << "TIMEWINDOW: " << QLocale::c().toString(model.header.timeWindowNs, 'g', 12) << '\n';
rad << "TIME INTERVAL: " << QLocale::c().toString(model.header.timeIntervalNs, 'g', 12) << '\n';
rad << "DISTANCE INTERVAL: " << QLocale::c().toString(model.header.distanceInc, 'g', 12) << '\n';
rad << "ANTENNAS: " << QLocale::c().toString(model.header.antennaFreq, 'g', 12) << '\n';
if (!model.header.antennaType.isEmpty()) rad << "ANTENNA: " << model.header.antennaType << '\n';
if (!model.header.date.isEmpty()) rad << "DATE: " << model.header.date << '\n';
if (!model.header.timeStr.isEmpty()) rad << "TIME: " << model.header.timeStr << '\n';
rad << "NUMBER_OF_CH: " << qMax(1, model.header.numberOfChannels) << '\n';
rad << "CH_X_OFFSETS: " << formatFloatList(model.header.chXOffsets) << '\n';
rad << "CH_Y_OFFSETS: " << formatFloatList(model.header.chYOffsets) << '\n';
if (!model.header.units.isEmpty()) rad << "UNITS: " << model.header.units << '\n';
rad << "START POSITION: " << QLocale::c().toString(model.header.startPosition, 'g', 12) << '\n';
const double stopPosition = model.header.stopPosition > model.header.startPosition
? model.header.stopPosition
: model.header.startPosition + model.tracesPerChannel * model.header.distanceInc;
rad << "STOP POSITION: " << QLocale::c().toString(stopPosition, 'g', 12) << '\n';
radFile.close();
return true;
}
QString IprhParser::formatFloatList(const QVector<float> &values)
{
QStringList parts;
parts.reserve(values.size());
for (float value : values) {
parts.append(QLocale::c().toString(value, 'g', 10));
}
return parts.join(' ');
}
QString IprhParser::uniqueConvertedBasePath(const QString &dirPath, const QString &baseName)
{
return QDir(dirPath).absoluteFilePath(baseName + QStringLiteral("_mala_converted"));
}
double IprhParser::extractDouble(const QString &value) {
bool ok = false;
double res = value.toDouble(&ok);
return ok ? res : 0.0;
}
int IprhParser::extractInt(const QString &value) {
bool ok = false;
int res = value.toInt(&ok);
return ok ? res : 0;
}

41
external/gpr3dviewer/IprhParser.h vendored Normal file
View File

@ -0,0 +1,41 @@
#ifndef IPRHPARSER_H
#define IPRHPARSER_H
#include "GPRDataModel.h"
#include <QString>
#include <QStringList>
class IprhParser {
public:
// 主入口:传入 .iprh 文件路径,自动寻找同名的 .iprb 文件
static bool loadFromIprh(const QString &iprhFilePath, GPRDataModel &model);
// 多通道 Impulse 数据合并入口:一个文件夹下有多个 _A01.iprh/.iprb ... _A14.iprh/.iprb
static bool loadImpulseMultiChannel(const QString &dirPath,
const QString &baseName,
GPRDataModel &model);
// 将多通道 Impulse 数据转换为本地 Mala Mira 风格 .rad/.rd3再走现有 Mala Mira 解析流程
static bool convertImpulseMultiChannelToMala(const QString &dirPath,
const QString &baseName,
QString *radFilePath,
QString *errorMessage = nullptr);
static bool parseHeaderOnly(const QString &iprhFilePath, GPRDataModel::Header &header);
private:
static bool parseIprhHeader(const QString &iprhFilePath, GPRDataModel::Header &header);
static bool loadIprbBinary(const QString &iprbFilePath, GPRDataModel &model);
static bool writeMalaFiles(const QString &radFilePath,
const QString &rd3FilePath,
const GPRDataModel &model,
const QString &sourceBaseName,
QString *errorMessage);
static QString formatFloatList(const QVector<float> &values);
static QString uniqueConvertedBasePath(const QString &dirPath, const QString &baseName);
static double extractDouble(const QString &value);
static int extractInt(const QString &value);
};
#endif // IPRHPARSER_H

View File

@ -0,0 +1,122 @@
#include "PerformanceLogger.h"
PerformanceLogger &PerformanceLogger::instance()
{
static PerformanceLogger logger;
return logger;
}
void PerformanceLogger::beginSession(const QString &sessionName)
{
QMutexLocker locker(&m_mutex);
m_sessionName = sessionName;
m_sessionActive = true;
}
void PerformanceLogger::endSession()
{
QMutexLocker locker(&m_mutex);
m_sessionActive = false;
}
void PerformanceLogger::log(const QString &category, const QString &name, qint64 elapsedMs)
{
QMutexLocker locker(&m_mutex);
Record rec;
rec.category = category;
rec.name = name;
rec.elapsedMs = elapsedMs;
rec.timestamp = QDateTime::currentDateTime();
m_records.append(rec);
}
QVector<PerformanceLogger::Record> PerformanceLogger::records() const
{
QMutexLocker locker(&m_mutex);
return m_records;
}
void PerformanceLogger::clear()
{
QMutexLocker locker(&m_mutex);
m_records.clear();
}
QString PerformanceLogger::reportString() const
{
QMutexLocker locker(&m_mutex);
QString report;
QTextStream ts(&report);
ts << "==================================================\n";
ts << "Performance Report";
if (!m_sessionName.isEmpty())
ts << " - " << m_sessionName;
ts << "\n";
ts << "Generated: " << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss") << "\n";
ts << "==================================================\n";
if (m_records.isEmpty()) {
ts << "No records.\n";
return report;
}
// Group by category
QString currentCategory;
qint64 categoryTotal = 0;
for (const auto &rec : m_records) {
if (rec.category != currentCategory) {
if (!currentCategory.isEmpty()) {
ts << " [Category Total: " << categoryTotal << " ms]\n\n";
}
currentCategory = rec.category;
categoryTotal = 0;
ts << "[" << currentCategory << "]\n";
}
ts << " " << rec.name << ": " << rec.elapsedMs << " ms\n";
categoryTotal += rec.elapsedMs;
}
if (!currentCategory.isEmpty()) {
ts << " [Category Total: " << categoryTotal << " ms]\n";
}
ts << "\n==================================================\n";
ts << "Grand Total: ";
qint64 grandTotal = 0;
for (const auto &rec : m_records)
grandTotal += rec.elapsedMs;
ts << grandTotal << " ms\n";
ts << "==================================================\n";
return report;
}
void PerformanceLogger::printReport() const
{
const QString report = reportString();
qDebug().noquote() << report;
}
void PerformanceLogger::saveToFile(const QString &filePath) const
{
QFile file(filePath);
if (!file.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Append)) {
qDebug() << "PerformanceLogger: failed to open" << filePath;
return;
}
QTextStream ts(&file);
ts << reportString();
file.close();
}
ScopedPerfTimer::ScopedPerfTimer(const QString &category, const QString &name)
: m_category(category), m_name(name)
{
m_timer.start();
}
ScopedPerfTimer::~ScopedPerfTimer()
{
const qint64 elapsed = m_timer.elapsed();
PerformanceLogger::instance().log(m_category, m_name, elapsed);
}

View File

@ -0,0 +1,60 @@
#ifndef PERFORMANCELOGGER_H
#define PERFORMANCELOGGER_H
#include <QString>
#include <QElapsedTimer>
#include <QVector>
#include <QFile>
#include <QTextStream>
#include <QDateTime>
#include <QDebug>
#include <QMutex>
#include <memory>
class PerformanceLogger
{
public:
struct Record
{
QString category;
QString name;
qint64 elapsedMs = 0;
QDateTime timestamp;
};
static PerformanceLogger &instance();
void beginSession(const QString &sessionName = QString());
void endSession();
void log(const QString &category, const QString &name, qint64 elapsedMs);
void printReport() const;
void saveToFile(const QString &filePath) const;
QString reportString() const;
QVector<Record> records() const;
void clear();
private:
PerformanceLogger() = default;
mutable QMutex m_mutex;
QVector<Record> m_records;
QString m_sessionName;
bool m_sessionActive = false;
};
class ScopedPerfTimer
{
public:
ScopedPerfTimer(const QString &category, const QString &name);
~ScopedPerfTimer();
private:
QString m_category;
QString m_name;
QElapsedTimer m_timer;
};
#define SCOPED_PERF_TIMER(category, name) ScopedPerfTimer _perfTimer(category, name)
#endif // PERFORMANCELOGGER_H

1148
external/gpr3dviewer/RadarProcessor.cpp vendored Normal file

File diff suppressed because it is too large Load Diff

346
external/gpr3dviewer/RadarProcessor.h vendored Normal file
View File

@ -0,0 +1,346 @@
/*
* RadarProcessor.h
* GPR数据预处理算法核心类头文件
*
* + 线
* MainWindow UI界面参数绑定GPRDataModel三维道集数据模型
* RAD/RD3格式多通道线
*/
#ifndef RADARPROCESSOR_H
#define RADARPROCESSOR_H
#include "GPRDataModel.h"
#include <QVector>
#include <QString>
#include <functional>
class RadarProcessor
{
public:
/**
* @brief
* 线
*/
enum class ProcStepType
{
StepTimeZeroCut, // 0零时校正切除直达波零点
StepZeroDrift, // 1道基线零漂消除
StepRemoveDC, // 2简易整道直流偏移扣除
StepBackgroundRemove, // 3剖面背景均值压制
StepSphericalTVG, // 4球面扩散+介质吸收TVG深度增益
StepBandpassFilter, // 5FIR带通滤波剔除高低频噪声
StepProfileSmooth, // 6时空剖面平滑降噪
StepInnerAGC, // 7单道内自适应振幅均衡AGC
StepTraceToTraceAGC, // 8道与道之间整体振幅均衡
StepHilbertTransform, // 9希尔伯特变换包络/正交分量)
StepGlobalGain, // 10全局统一倍率增益放大
StepMigration // 11Kirchhoff偏移
};
//=========================================================================
// 1. 零时切除参数结构体 StepTimeZeroCut
// 作用:定位直达波起跳零点,裁掉零点前无效采样点,统一道起始位置
//=========================================================================
struct TimeZeroCutParams
{
bool autoDetect = true; // true=自动识别零点false=手动固定裁掉前N个采样
int cutSamples = 30; // 手动模式裁切采样点数
int frontSearchWindow = 180; // 零点搜索窗口只在道前180个采样内找起跳点
double noiseSigmaMultiple = 3.0; // 起跳判定阈值:噪声标准差倍数,大于判定为有效直达波
bool enable = false; // 流水线开关:是否启用本步骤
};
/**
* @brief
*/
enum class ZeroDriftMode
{
DC, // 模式1整道均值直流偏移一次性扣除速度快
Sliding // 模式2滑动窗口逐段扣除基线适合长时漂移严重数据
};
struct ZeroDriftParams
{
ZeroDriftMode mode = ZeroDriftMode::Sliding;
int slidingWindowSamples = 100; // 滑动窗口采样点数仅Sliding模式生效
bool enable = false;
};
//=========================================================================
// TVG球面深度增益 StepSphericalTVG
// 物理原理:电磁波随传播距离球面扩散衰减+介质吸收衰减
//=========================================================================
struct SphericalTvgParams {
bool enableSpherical = true; // 开启球面扩散补偿
bool enableAbsorption = true; // 开启介质吸收衰减补偿
double velocityMPerNs = 0.12; // 雷达波速(m/ns)从GPRDataHeader自动读取
double referenceDepthM = 0.01; // 参考基准深度(增益归一化基准)
double exponent = 1.5; // 扩散指数空气1.0、土体/混凝土常用1.0~2.0
double absorptionBeta = 1.0; // 吸收系数 β (dB/m),干土小、湿土大
double maxGain = 30.0; // 最大增益上限,防止深层噪声过度放大
bool enable = false;
};
//=========================================================================
// 带通滤波 StepBandpassFilter
// 滤除天线主频外低频漂移、高频仪器白噪声
//=========================================================================
struct BandpassParams
{
bool autoFreq = true; // true自动根据天线主频计算通带false手动填高低频
double lowFreqHz = 1000; // 通带下限(Hz)
double highFreqHz = 100; // 通带上限(Hz)
double antennaFreqMHz = 200.0; // 天线中心频率(MHz),自动模式读取头文件
bool enable = false;
};
//=========================================================================
// 剖面二维平滑 StepProfileSmooth
// 垂直=深度方向平滑;水平=测线道方向平滑
//=========================================================================
struct ProfileSmoothParams
{
int smoothWindow = 2; // 半窗口大小,实际总窗口=2*window+1
bool verticalSmooth = true; // 深度方向平滑开关
bool horizontalSmooth = true; // 测线道方向平滑开关
bool enable = false;
};
//=========================================================================
// 背景去除 StepBackgroundRemove
// 压制整剖面静态杂波、地面耦合固定反射、钢筋/管线持续干扰
//=========================================================================
enum class BackgroundMode
{
MeanAverage, // 均值法:多道平均作为背景相减,通用稳定
SingularityFilter // 奇异值滤波:保留强异常,压制平缓背景(空洞/管线优选)
};
struct BackgroundParams
{
BackgroundMode mode = BackgroundMode::MeanAverage;
int averageTraceCount = 301; // 参与平均背景的同通道道数量
double singularityThreshold = 1.8; // 奇异判定阈值,越大越不容易抹掉小异常
bool enable = false;
};
//=========================================================================
// 道内AGC StepInnerAGC
// 同一道内深浅振幅均衡消除天然深度衰减TVG补充均衡
//=========================================================================
struct TraceInnerAgcParams {
int windowSamples = 120; // 道内滑动RMS统计窗口
double targetRms = 0.0; // 目标RMS值0=自适应均衡
double maxGain = 6.0; // 单窗口最大增益限制
double epsilon = 1.0; // 防除零极小值
bool enable = false;
};
//=========================================================================
// 道间AGC StepTraceToTraceAGC
// 整条测线所有道之间振幅拉平,消除收发耦合波动、行走速度不均
//=========================================================================
enum class TraceToTraceMode
{
Global, // 全局:全部道统一均衡基准
Local // 局部:滑动窗口相邻道均衡(起伏大路面优选)
};
struct TraceToTraceAgcParams {
TraceToTraceMode mode = TraceToTraceMode::Local;
int horizontalWindowTraces = 31; // 局部模式滑动道窗口
double targetRms = 0.0;
double maxGain = 4.0;
double epsilon = 1.0;
bool enable = false;
};
//=========================================================================
// 希尔伯特变换 StepHilbertTransform
// 输出振幅包络(成像常用)或正交虚部(相位分析)
//=========================================================================
struct HilbertParams
{
bool computeEnvelope = true; // true=振幅包络false=正交相位分量
bool enable = false;
};
// 全局倍率增益
struct GlobalGainParams
{
double factor = 1.0; // 放大倍数0.1~10区间常用
bool enable = false;
};
// 简易去直流(轻量零漂)
struct DcShiftParams
{
bool enable = false;
};
//=========================================================================
// Kirchhoff偏移参数
//=========================================================================
struct MigrationParams
{
int sumWidth = 64; // 两侧各求和的迹线数量
double velocityMPerNs = 0.12; // 雷达波速(m/ns)
bool enable = false;
};
//=========================================================================
// 流水线单步单元:存储步骤类型 + 全套独立参数
//=========================================================================
struct ProcStepUnit
{
ProcStepType type;
// 全部算法参数容器运行时仅type对应的结构体生效
TimeZeroCutParams zeroCut;
SphericalTvgParams tvg;
BandpassParams band;
ProfileSmoothParams smooth;
BackgroundParams bg;
ZeroDriftParams zeroDrift;
TraceInnerAgcParams innerAgc;
TraceToTraceAgcParams traceAgc;
HilbertParams hilbert;
GlobalGainParams gain;
DcShiftParams dc;
MigrationParams migration;
};
//=========================================================================
// 完整预处理流水线:有序步骤队列 + 默认标准流程
//=========================================================================
struct ProcPipeline
{
QVector<ProcStepUnit> steps;
/**
* @brief 线
* TVG增益AGCAGC
*/
void setDefaultFlow()
{
steps.clear();
// 1. 零时校正
ProcStepUnit s0;
s0.type = ProcStepType::StepTimeZeroCut;
s0.zeroCut.autoDetect = true;
s0.zeroCut.cutSamples = 30;
s0.zeroCut.frontSearchWindow = 180;
s0.zeroCut.noiseSigmaMultiple = 3.0;
s0.zeroCut.enable = true;
steps.append(s0);
// 2. 零漂去除
ProcStepUnit s1;
s1.type = ProcStepType::StepZeroDrift;
s1.zeroDrift.mode = ZeroDriftMode::Sliding;
s1.zeroDrift.slidingWindowSamples = 100;
s1.zeroDrift.enable = true;
steps.append(s1);
// 3. 背景压制
ProcStepUnit s2;
s2.type = ProcStepType::StepBackgroundRemove;
s2.bg.mode = BackgroundMode::MeanAverage;
s2.bg.averageTraceCount = 301;
s2.bg.enable = true;
steps.append(s2);
// 4. TVG深度增益
ProcStepUnit s3;
s3.type = ProcStepType::StepSphericalTVG;
s3.tvg.enableSpherical = true;
s3.tvg.enableAbsorption = true;
s3.tvg.exponent = 1.0;
s3.tvg.maxGain = 20.0;
s3.tvg.absorptionBeta = 0.002;
s3.tvg.enable = true;
steps.append(s3);
// 5. 带通滤波
ProcStepUnit s4;
s4.type = ProcStepType::StepBandpassFilter;
s4.band.autoFreq = true;
s4.band.lowFreqHz = 100;
s4.band.highFreqHz = 1500;
s4.band.enable = true;
steps.append(s4);
// 6. 剖面平滑
ProcStepUnit s5;
s5.type = ProcStepType::StepProfileSmooth;
s5.smooth.smoothWindow = 2;
s5.smooth.verticalSmooth = true;
s5.smooth.horizontalSmooth = true;
s5.smooth.enable = true;
steps.append(s5);
// 7. 道内AGC
ProcStepUnit s6;
s6.type = ProcStepType::StepInnerAGC;
s6.innerAgc.windowSamples = 120;
s6.innerAgc.maxGain = 6.0;
s6.innerAgc.enable = true;
steps.append(s6);
// 8. 道间均衡AGC
ProcStepUnit s7;
s7.type = ProcStepType::StepTraceToTraceAGC;
s7.traceAgc.mode = TraceToTraceMode::Local;
s7.traceAgc.horizontalWindowTraces = 31;
s7.traceAgc.maxGain = 4.0;
s7.traceAgc.enable = true;
steps.append(s7);
}
};
public:
struct PipelineRuntimeContext {
std::function<bool()> isCanceled;
std::function<void(int stepIndex, int stepCount, const QString &stepName)> reportProgress;
};
RadarProcessor();
//==================== 单算法静态处理函数 ====================
/// 简易去除整道直流偏移
static GPRDataModel removeDCShift(const GPRDataModel &input, const DcShiftParams& param);
/// 零漂基线校正
static GPRDataModel removeZeroDrift(const GPRDataModel &input, const ZeroDriftParams& param);
/// 全局统一倍率增益
static GPRDataModel applyGain(const GPRDataModel &input, const GlobalGainParams& param);
/// TVG球面+吸收深度增益补偿
static GPRDataModel sphericalTvg(const GPRDataModel &input, const SphericalTvgParams &params);
/// 二维剖面均值平滑
static GPRDataModel profileSmooth(const GPRDataModel &input, const ProfileSmoothParams &params);
/// 道内滑动AGC均衡
static GPRDataModel traceInnerAgc(const GPRDataModel &input, const TraceInnerAgcParams &params);
/// 道与道之间振幅均衡
static GPRDataModel traceToTraceEqualization(const GPRDataModel &input, const TraceToTraceAgcParams &params);
/// 希尔伯特变换包络/相位
static GPRDataModel hilbertTransform(const GPRDataModel &input, const HilbertParams &params);
/// 零点自动/手动裁切算法实现
static GPRDataModel cutTimeZero(const GPRDataModel& input, const TimeZeroCutParams& param);
/// FIR带通滤波实现
static GPRDataModel bandpassFilter(const GPRDataModel &input, const BandpassParams &params);
/// 背景杂波去除实现(均值/奇异滤波双模式)
static GPRDataModel backgroundRemoval(const GPRDataModel &input, const BackgroundParams &params);
/// Kirchhoff偏移x-t域双曲线求和
static GPRDataModel kirchhoffMigration(const GPRDataModel &input, const MigrationParams &params);
/**
* @brief 线
* @param rawData GPR三维道集数据
* @param pipeline 线
* @return
*/
static GPRDataModel runPipeline(const GPRDataModel& rawData, const ProcPipeline& pipeline);
static GPRDataModel runPipeline(const GPRDataModel& rawData,
const ProcPipeline& pipeline,
const PipelineRuntimeContext *context);
};
#endif // RADARPROCESSOR_H

9
external/gpr3dviewer/RadarTypes.h vendored Normal file
View File

@ -0,0 +1,9 @@
#ifndef RADARTYPES_H
#define RADARTYPES_H
enum class RadarType {
Mala_Mira, // .rad + .rd3
Impulse // .iprh + .iprb
};
#endif // RADARTYPES_H

300
external/gpr3dviewer/Rd3Parser.cpp vendored Normal file
View File

@ -0,0 +1,300 @@
#include "Rd3Parser.h"
#include "PerformanceLogger.h"
#include <QFile>
#include <QTextStream>
#include <QDataStream>
#include <QDebug>
#include <QDir>
#include <QFileInfo>
#include <QRegularExpression>
#include <QStringList>
#include <QVector3D>
/**
* @brief Mala Mira系列雷达RD3数据rad头文件 + rd3纯二进制波形
* @param radFilePath .rad文本配置头文件路径
* @param model GPR全局数据模型
* @return true false
* @note .rad存储全部仪器/线.rd3无文件头short16振幅采样
*/
bool Rd3Parser::loadFromRad(const QString &radFilePath, GPRDataModel &model) {
SCOPED_PERF_TIMER("Parser", "Rd3Parser::loadFromRad");
// 清空旧数据与头部信息
model.clear();
model.header = GPRDataModel::Header{};
// 第一步解析rad文本头填充全部采集参数
if (!parseRadHeader(radFilePath, model.header)) {
qDebug() << "Error: Failed to parse .rad header file:" << radFilePath;
return false;
}
// 自动匹配同目录同名二进制数据文件 .rd3
QFileInfo radInfo(radFilePath);
QString binaryPath = radInfo.absolutePath() + "/" + radInfo.completeBaseName() + ".rd3";
if (!QFile::exists(binaryPath)) {
qDebug() << "Error: Matching binary file not found at:" << binaryPath;
return false;
}
// 推算后仍无有效道数则失败
if (model.header.numTraces <= 0) {
qDebug() << "Error: Unable to determine numTraces from header";
return false;
}
// 第二步:读取无头部的纯二进制波形
return loadRd3Binary(binaryPath, model);
}
/**
* @brief RAD文本头文件Header结构体
* @param radFilePath rad文本文件路径
* @param header
* @return true
*/
bool Rd3Parser::parseRadHeader(const QString &radFilePath, GPRDataModel::Header &header) {
SCOPED_PERF_TIMER("Parser", "Rd3Parser::parseRadHeader");
QFile file(radFilePath);
// 只读文本模式打开
if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) {
qDebug() << "Error: Cannot open rad file for read";
return false;
}
QTextStream in(&file);
// 逐行读取键值对 KEY: VALUE 或 KEY=VALUE
while (!in.atEnd()) {
QString line = in.readLine().trimmed();
if (line.isEmpty()) continue;
int sepPos = line.indexOf(':');
if (sepPos == -1) sepPos = line.indexOf('=');
if (sepPos == -1) continue;
QString key = line.left(sepPos).trimmed();
QString value = line.mid(sepPos + 1).trimmed();
// 原始字符串参数全存入map备用
header.rawParams[key] = value;
// 识别核心业务参数并类型转换赋值
if (key == "SAMPLES") {
header.samplesPerTrace = extractInt(value);
} else if (key == "LAST TRACE") {
header.numTraces = extractInt(value);
} else if (key == "TIMEWINDOW") {
header.timeWindowNs = extractDouble(value);
} else if (key == "DISTANCE INTERVAL") {
header.distanceInc = extractDouble(value);
} else if (key == "ANTENNAS") {
header.antennaFreq = extractDouble(value);
} else if (key == "ANTENNAS") {
header.antennaType = value;
} else if (key == "DATE") {
header.date = value;
} else if (key == "TIME") {
header.timeStr = value;
} else if (key == "NUMBER_OF_CH") {
header.numberOfChannels = extractInt(value);
} else if (key == "CH_X_OFFSETS") {
// 多通道天线水平X偏移量
QStringList offsets = value.split(' ', Qt::SkipEmptyParts);
for (const QString &offset : offsets) {
header.chXOffsets.append(offset.toFloat());
}
} else if (key == "CH_Y_OFFSETS") {
// 多通道天线行进Y偏移量
QStringList offsets = value.split(' ', Qt::SkipEmptyParts);
for (const QString &offset : offsets) {
header.chYOffsets.append(offset.toFloat());
}
} else if (key == "UNITS") {
header.units = value;
} else if (key == "START POSITION") {
header.startPosition = extractDouble(value);
} else if (key == "STOP POSITION") {
header.stopPosition = extractDouble(value);
} else if (key == "TIME INTERVAL") {
header.timeIntervalNs = extractDouble(value);
}
}
file.close();
// 基础合法性校验单道采样数必须大于0总道数允许从二进制文件推算
if (header.samplesPerTrace <= 0) {
qDebug() << "Error: Invalid SAMPLES value in rad file";
return false;
}
// 设置地层电磁波默认速度 0.1 m/ns空气近似值后续界面可手动修改
header.waveVelocity = 0.1;
qDebug() << "==== RAD Header Parse Complete ===="
<< "\nSamples per trace:" << header.samplesPerTrace
<< "\nTotal traces:" << header.numTraces
<< "\nTime window(ns):" << header.timeWindowNs
<< "\nChannel count:" << header.numberOfChannels
<< "\nX offset array size:" << header.chXOffsets.size()
<< "\nY offset array size:" << header.chYOffsets.size();
return true;
}
/**
* @brief RD3二进制数据short16振幅
* @param rd3FilePath rd3波形文件路径
* @param model traces波形数组
* @return true
* @detail
* 00,01...0N | 10...1N | ...
* short 16
*/
bool Rd3Parser::loadRd3Binary(const QString &rd3FilePath, GPRDataModel &model) {
SCOPED_PERF_TIMER("Parser", "Rd3Parser::loadRd3Binary");
QFile file(rd3FilePath);
if (!file.open(QIODevice::ReadOnly)) {
qDebug() << "Error: Open rd3 binary failed:" << rd3FilePath;
return false;
}
const int samplesPerTrace = model.header.samplesPerTrace;
const int totalTraceCount = model.header.numTraces;
// ========== 关键修正rd3无头部偏移直接从文件0位置开始读波形 ==========
const qint64 dataStartOffset = 0;
file.seek(dataStartOffset);
// 理论完整字节大小校验
const qint64 singleTraceByteSize = samplesPerTrace * sizeof(short);
const qint64 fullExpectedBytes = totalTraceCount * singleTraceByteSize;
const qint64 realFileBytes = file.size();
if (realFileBytes < fullExpectedBytes) {
qDebug() << "Warning: RD3 file size smaller than theoretical data size!"
<< "Expected:" << fullExpectedBytes << "Actual:" << realFileBytes;
}
// 预分配内存容器,减少动态扩容开销
model.traces.reserve(totalTraceCount);
// 单道读取缓冲区
QByteArray traceBuffer;
traceBuffer.resize(singleTraceByteSize);
// 通道数量兜底
const int channelCnt = model.header.numberOfChannels > 0 ? model.header.numberOfChannels : 1;
model.channels = channelCnt;
// 每通道独立道数量
model.tracesPerChannel = totalTraceCount / channelCnt;
// 计算整条测线总行进距离
if (model.header.distanceInc > 1e-6) {
model.totalDistance = static_cast<float>(model.tracesPerChannel * model.header.distanceInc);
} else {
model.totalDistance = static_cast<float>(model.header.stopPosition - model.header.startPosition);
}
// 循环逐道读取二进制波形
for (int traceGlobalIdx = 0; traceGlobalIdx < totalTraceCount; ++traceGlobalIdx) {
if (file.atEnd()) {
qDebug() << "Warning: File ended early at global trace index" << traceGlobalIdx;
break;
}
// 读取一整道所有采样字节
qint64 readBytes = file.read(traceBuffer.data(), traceBuffer.size());
if (readBytes != traceBuffer.size()) {
qDebug() << "Warning: Trace" << traceGlobalIdx << "incomplete byte read";
break;
}
RadarTrace oneTrace;
oneTrace.amplitudes.resize(samplesPerTrace);
// 二进制指针强转shortWindows平台原生小端序雷达标准字节序
const short* rawShortBuf = reinterpret_cast<const short*>(traceBuffer.constData());
// 填充单道振幅值
for (int s = 0; s < samplesPerTrace; s++) {
oneTrace.amplitudes[s] = rawShortBuf[s];
}
// 分配所属通道、通道内道序号
int chNo = traceGlobalIdx % channelCnt;
int traceInChIdx = traceGlobalIdx / channelCnt;
oneTrace.channelNumber = chNo;
// 计算空间三维坐标 X(天线横向偏移) Y(行进里程) Z=0地面平面
float xOff = 0.0f;
float yOff = 0.0f;
if (chNo < model.header.chXOffsets.size()) xOff = model.header.chXOffsets[chNo];
if (chNo < model.header.chYOffsets.size()) yOff = model.header.chYOffsets[chNo];
float lineDist = static_cast<float>(model.header.startPosition + traceInChIdx * model.header.distanceInc);
oneTrace.position = QVector3D(xOff, lineDist - yOff, 0.0f);
model.traces.append(std::move(oneTrace));
}
file.close();
if (model.traces.isEmpty()) {
qDebug() << "Error: No valid traces loaded from rd3";
return false;
}
qDebug() << "RD3 Binary Load OK, total valid traces:" << model.traces.size();
return true;
}
/**
* @brief double数值0
* @param value 200 MHz shielded
* @return
*/
double Rd3Parser::extractDouble(const QString& value) {
bool ok = false;
double res = 0.0;
// 1. 遍历字符串,截取开头连续的数字/小数点部分
int numLength = 0;
while (numLength < value.length()) {
QChar ch = value.at(numLength);
// 只保留 数字 和 小数点(支持小数,如 200.5 MHz
if (ch.isDigit() || ch == '.') {
numLength++;
}
else {
// 遇到非数字/小数点,停止截取
break;
}
}
// 2. 截取有效数字字符串并转换
if (numLength > 0) {
QString numStr = value.left(numLength);
res = numStr.toDouble(&ok);
}
// 3. 转换成功返回数值失败返回0.0
return ok ? res : 0.0;
}
/**
* @brief int整型0
* @param value rad文件字符串数值
* @return
*/
int Rd3Parser::extractInt(const QString &value) {
bool ok = false;
int res = value.toInt(&ok);
return ok ? res : 0;
}

21
external/gpr3dviewer/Rd3Parser.h vendored Normal file
View File

@ -0,0 +1,21 @@
#ifndef RD3PARSER_H
#define RD3PARSER_H
#include "GPRDataModel.h"
#include <QString>
class Rd3Parser {
public:
// 主入口:传入 .rad 文件路径,它会自动寻找同名的 .rd3 文件
static bool loadFromRad(const QString &radFilePath, GPRDataModel &model);
private:
static bool parseRadHeader(const QString &radFilePath, GPRDataModel::Header &header);
static bool loadRd3Binary(const QString &rd3FilePath, GPRDataModel &model);
// 辅助:从 .rad 的字符串中提取数值
static double extractDouble(const QString &value);
static int extractInt(const QString &value);
};
#endif // RD3PARSER_H

127
external/gpr3dviewer/SurveyGeometry.h vendored Normal file
View File

@ -0,0 +1,127 @@
#ifndef SURVEYGEOMETRY_H
#define SURVEYGEOMETRY_H
#include <QJsonArray>
#include <QJsonObject>
#include <QVector>
struct SurveyGeometry {
double rtkOffsetX = 0.0; // RTK相对天线中心的设备x轴偏移m头文件无对应字段保留给用户输入
double rtkOffsetY = 0.0; // RTK相对天线中心的设备y轴偏移m来自CH_Y_OFFSETS
double ch1XRel = 0.0; // 第1通道相对天线中心的x坐标m由CH_X_OFFSETS首尾值推导
QVector<double> channelXRel; // 各通道相对天线中心的x坐标m来自CH_X_OFFSETS归一到天线中心
int channelCount = 1; // 通道总数来自NUMBER_OF_CH
int cgcsZone = 0; // CGCS2000 3度带带号0表示自动检测
double centralMeridianDeg = 0.0; // 中央经线(度),自动推导
bool operator==(const SurveyGeometry &other) const {
return rtkOffsetX == other.rtkOffsetX &&
rtkOffsetY == other.rtkOffsetY &&
ch1XRel == other.ch1XRel &&
channelXRel == other.channelXRel &&
channelCount == other.channelCount &&
cgcsZone == other.cgcsZone &&
centralMeridianDeg == other.centralMeridianDeg;
}
void applyHeaderOffsets(int numberOfChannels,
const QVector<float> &chXOffsets,
const QVector<float> &chYOffsets)
{
const int offsetChannelCount = qMax(chXOffsets.size(), chYOffsets.size());
channelCount = qMax(1, numberOfChannels > 0 ? numberOfChannels : offsetChannelCount);
rtkOffsetY = chYOffsets.isEmpty() ? 0.0 : static_cast<double>(chYOffsets.first());
channelXRel.clear();
channelXRel.reserve(channelCount);
if (chXOffsets.size() >= 2) {
const double centerX = (static_cast<double>(chXOffsets.first()) + chXOffsets.last()) / 2.0;
for (int i = 0; i < chXOffsets.size() && i < channelCount; ++i) {
channelXRel.append(static_cast<double>(chXOffsets[i]) - centerX);
}
ch1XRel = channelXRel.isEmpty() ? 0.0 : channelXRel.first();
} else {
ch1XRel = 0.0;
}
if (channelXRel.size() < channelCount) {
const int existing = channelXRel.size();
channelXRel.resize(channelCount);
if (existing == 0) {
for (int c = 0; c < channelCount; ++c) {
channelXRel[c] = 0.0;
}
} else if (channelCount > 1) {
const double lastXRel = -ch1XRel;
for (int c = existing; c < channelCount; ++c) {
channelXRel[c] = ch1XRel + (lastXRel - ch1XRel) * c / (channelCount - 1.0);
}
}
}
}
static SurveyGeometry fromHeaderOffsets(int numberOfChannels,
const QVector<float> &chXOffsets,
const QVector<float> &chYOffsets)
{
SurveyGeometry g;
g.applyHeaderOffsets(numberOfChannels, chXOffsets, chYOffsets);
return g;
}
QJsonObject toJson() const {
QJsonObject obj;
obj["rtkOffsetX"] = rtkOffsetX;
obj["rtkOffsetY"] = rtkOffsetY;
obj["ch1XRel"] = ch1XRel;
QJsonArray channelXRelArray;
for (double x : channelXRel) {
channelXRelArray.append(x);
}
obj["channelXRel"] = channelXRelArray;
obj["channelCount"] = channelCount;
obj["cgcsZone"] = cgcsZone;
obj["centralMeridianDeg"] = centralMeridianDeg;
return obj;
}
static SurveyGeometry fromJson(const QJsonObject &obj) {
SurveyGeometry g;
g.rtkOffsetX = obj.value("rtkOffsetX").toDouble(0.0);
g.rtkOffsetY = obj.value("rtkOffsetY").toDouble(0.0);
g.ch1XRel = obj.value("ch1XRel").toDouble(0.0);
g.channelCount = qMax(1, obj.value("channelCount").toInt(1));
g.cgcsZone = obj.value("cgcsZone").toInt(0);
g.centralMeridianDeg = obj.value("centralMeridianDeg").toDouble(0.0);
const QJsonArray channelXRelArray = obj.value("channelXRel").toArray();
for (const QJsonValue &value : channelXRelArray) {
g.channelXRel.append(value.toDouble(0.0));
}
if (g.channelXRel.isEmpty() && obj.contains("ch16XRel")) {
const double legacyLastXRel = obj.value("ch16XRel").toDouble(-g.ch1XRel);
g.channelXRel.resize(g.channelCount);
for (int c = 0; c < g.channelCount; ++c) {
g.channelXRel[c] = (g.channelCount > 1)
? g.ch1XRel + (legacyLastXRel - g.ch1XRel) * c / (g.channelCount - 1.0)
: g.ch1XRel;
}
} else if (g.channelXRel.isEmpty()) {
g.channelXRel.resize(g.channelCount);
const double lastXRel = -g.ch1XRel;
for (int c = 0; c < g.channelCount; ++c) {
g.channelXRel[c] = (g.channelCount > 1)
? g.ch1XRel + (lastXRel - g.ch1XRel) * c / (g.channelCount - 1.0)
: g.ch1XRel;
}
} else if (g.channelXRel.size() != g.channelCount) {
g.channelCount = qMax(1, g.channelXRel.size());
}
return g;
}
};
#endif // SURVEYGEOMETRY_H

View File

@ -0,0 +1,189 @@
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
/* kiss_fft.h
defines kiss_fft_scalar as either short or a float type
and defines
typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */
#ifndef _kiss_fft_guts_h
#define _kiss_fft_guts_h
#include "kiss_fft.h"
#include "kiss_fft_log.h"
#include <limits.h>
#define MAXFACTORS 32
/* e.g. an fft of length 128 has 4 factors
as far as kissfft is concerned
4*4*4*2
*/
struct kiss_fft_state{
int nfft;
int inverse;
int factors[2*MAXFACTORS];
kiss_fft_cpx twiddles[1];
};
/*
Explanation of macros dealing with complex math:
C_MUL(m,a,b) : m = a*b
C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise
C_SUB( res, a,b) : res = a - b
C_SUBFROM( res , a) : res -= a
C_ADDTO( res , a) : res += a
* */
#ifdef FIXED_POINT
#include <stdint.h>
#if (FIXED_POINT==32)
# define FRACBITS 31
# define SAMPPROD int64_t
#define SAMP_MAX INT32_MAX
#define SAMP_MIN INT32_MIN
#else
# define FRACBITS 15
# define SAMPPROD int32_t
#define SAMP_MAX INT16_MAX
#define SAMP_MIN INT16_MIN
#endif
#if defined(CHECK_OVERFLOW)
# define CHECK_OVERFLOW_OP(a,op,b) \
if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
KISS_FFT_WARNING("overflow (%d " #op" %d) = %ld", (a),(b),(SAMPPROD)(a) op (SAMPPROD)(b)); }
#endif
# define smul(a,b) ( (SAMPPROD)(a)*(b) )
# define sround( x ) (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS )
# define S_MUL(a,b) sround( smul(a,b) )
# define C_MUL(m,a,b) \
do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
(m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
# define DIVSCALAR(x,k) \
(x) = sround( smul( x, SAMP_MAX/k ) )
# define C_FIXDIV(c,div) \
do { DIVSCALAR( (c).r , div); \
DIVSCALAR( (c).i , div); }while (0)
# define C_MULBYSCALAR( c, s ) \
do{ (c).r = sround( smul( (c).r , s ) ) ;\
(c).i = sround( smul( (c).i , s ) ) ; }while(0)
#else /* not FIXED_POINT*/
# define S_MUL(a,b) ( (a)*(b) )
#define C_MUL(m,a,b) \
do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
(m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
# define C_FIXDIV(c,div) /* NOOP */
# define C_MULBYSCALAR( c, s ) \
do{ (c).r *= (s);\
(c).i *= (s); }while(0)
#endif
#ifndef CHECK_OVERFLOW_OP
# define CHECK_OVERFLOW_OP(a,op,b) /* noop */
#endif
#define C_ADD( res, a,b)\
do { \
CHECK_OVERFLOW_OP((a).r,+,(b).r)\
CHECK_OVERFLOW_OP((a).i,+,(b).i)\
(res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \
}while(0)
#define C_SUB( res, a,b)\
do { \
CHECK_OVERFLOW_OP((a).r,-,(b).r)\
CHECK_OVERFLOW_OP((a).i,-,(b).i)\
(res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \
}while(0)
#define C_ADDTO( res , a)\
do { \
CHECK_OVERFLOW_OP((res).r,+,(a).r)\
CHECK_OVERFLOW_OP((res).i,+,(a).i)\
(res).r += (a).r; (res).i += (a).i;\
}while(0)
#define C_SUBFROM( res , a)\
do {\
CHECK_OVERFLOW_OP((res).r,-,(a).r)\
CHECK_OVERFLOW_OP((res).i,-,(a).i)\
(res).r -= (a).r; (res).i -= (a).i; \
}while(0)
#ifdef FIXED_POINT
# define KISS_FFT_COS(phase) floor(.5+SAMP_MAX * cos (phase))
# define KISS_FFT_SIN(phase) floor(.5+SAMP_MAX * sin (phase))
# define HALF_OF(x) ((x)>>1)
#elif defined(USE_SIMD)
#if defined(HAVE_LASX)
#define KISS_FFT_COS(phase) ({ \
float __cos_val = cosf(phase); \
(__m256)(__lasx_xvldrepl_w(&__cos_val, 0)); \
})
#define KISS_FFT_SIN(phase) ({ \
float __sin_val = sinf(phase); \
(__m256)(__lasx_xvldrepl_w(&__sin_val, 0)); \
})
#define HALF_OF(x) ((x) * (__m256)(__lasx_xvreplgr2vr_w(0x3F000000))) // 0.5f
#elif defined(HAVE_LSX)
#define KISS_FFT_COS(phase) ({ \
float __cos_val = cosf(phase); \
(__m128)(__lsx_vldrepl_w(&__cos_val, 0)); \
})
#define KISS_FFT_SIN(phase) ({ \
float __sin_val = sinf(phase); \
(__m128)(__lsx_vldrepl_w(&__sin_val, 0)); \
})
#define HALF_OF(x) ((x) * (__m128)(__lsx_vreplgr2vr_w(0x3F000000))) // 0.5f
#else
# define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) )
# define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) )
# define HALF_OF(x) ((x)*_mm_set1_ps(.5))
#endif
#else
# define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase)
# define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase)
# define HALF_OF(x) ((x)*((kiss_fft_scalar).5))
#endif
#define kf_cexp(x,phase) \
do{ \
(x)->r = KISS_FFT_COS(phase);\
(x)->i = KISS_FFT_SIN(phase);\
}while(0)
/* a debugging function */
#define pcpx(c)\
KISS_FFT_DEBUG("%g + %gi\n",(double)((c)->r),(double)((c)->i))
#ifdef KISS_FFT_USE_ALLOCA
// define this to allow use of alloca instead of malloc for temporary buffers
// Temporary buffers are used in two case:
// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5
// 2. "in-place" FFTs. Notice the quotes, since kissfft does not really do an in-place transform.
#include <alloca.h>
#define KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes)
#define KISS_FFT_TMP_FREE(ptr)
#else
#define KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes)
#define KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr)
#endif
#endif /* _kiss_fft_guts_h */

View File

@ -0,0 +1,424 @@
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#include <stdint.h>
#include "_kiss_fft_guts.h"
/* The guts header contains all the multiplication and addition macros that are defined for
fixed or floating point complex numbers. It also delares the kf_ internal functions.
*/
static void kf_bfly2(
kiss_fft_cpx * Fout,
const size_t fstride,
const kiss_fft_cfg st,
int m
)
{
kiss_fft_cpx * Fout2;
kiss_fft_cpx * tw1 = st->twiddles;
kiss_fft_cpx t;
Fout2 = Fout + m;
do{
C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2);
C_MUL (t, *Fout2 , *tw1);
tw1 += fstride;
C_SUB( *Fout2 , *Fout , t );
C_ADDTO( *Fout , t );
++Fout2;
++Fout;
}while (--m);
}
static void kf_bfly4(
kiss_fft_cpx * Fout,
const size_t fstride,
const kiss_fft_cfg st,
const size_t m
)
{
kiss_fft_cpx *tw1,*tw2,*tw3;
kiss_fft_cpx scratch[6];
size_t k=m;
const size_t m2=2*m;
const size_t m3=3*m;
tw3 = tw2 = tw1 = st->twiddles;
do {
C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4);
C_MUL(scratch[0],Fout[m] , *tw1 );
C_MUL(scratch[1],Fout[m2] , *tw2 );
C_MUL(scratch[2],Fout[m3] , *tw3 );
C_SUB( scratch[5] , *Fout, scratch[1] );
C_ADDTO(*Fout, scratch[1]);
C_ADD( scratch[3] , scratch[0] , scratch[2] );
C_SUB( scratch[4] , scratch[0] , scratch[2] );
C_SUB( Fout[m2], *Fout, scratch[3] );
tw1 += fstride;
tw2 += fstride*2;
tw3 += fstride*3;
C_ADDTO( *Fout , scratch[3] );
if(st->inverse) {
Fout[m].r = scratch[5].r - scratch[4].i;
Fout[m].i = scratch[5].i + scratch[4].r;
Fout[m3].r = scratch[5].r + scratch[4].i;
Fout[m3].i = scratch[5].i - scratch[4].r;
}else{
Fout[m].r = scratch[5].r + scratch[4].i;
Fout[m].i = scratch[5].i - scratch[4].r;
Fout[m3].r = scratch[5].r - scratch[4].i;
Fout[m3].i = scratch[5].i + scratch[4].r;
}
++Fout;
}while(--k);
}
static void kf_bfly3(
kiss_fft_cpx * Fout,
const size_t fstride,
const kiss_fft_cfg st,
size_t m
)
{
size_t k=m;
const size_t m2 = 2*m;
kiss_fft_cpx *tw1,*tw2;
kiss_fft_cpx scratch[5];
kiss_fft_cpx epi3;
epi3 = st->twiddles[fstride*m];
tw1=tw2=st->twiddles;
do{
C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
C_MUL(scratch[1],Fout[m] , *tw1);
C_MUL(scratch[2],Fout[m2] , *tw2);
C_ADD(scratch[3],scratch[1],scratch[2]);
C_SUB(scratch[0],scratch[1],scratch[2]);
tw1 += fstride;
tw2 += fstride*2;
Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
C_MULBYSCALAR( scratch[0] , epi3.i );
C_ADDTO(*Fout,scratch[3]);
Fout[m2].r = Fout[m].r + scratch[0].i;
Fout[m2].i = Fout[m].i - scratch[0].r;
Fout[m].r -= scratch[0].i;
Fout[m].i += scratch[0].r;
++Fout;
}while(--k);
}
static void kf_bfly5(
kiss_fft_cpx * Fout,
const size_t fstride,
const kiss_fft_cfg st,
int m
)
{
kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
int u;
kiss_fft_cpx scratch[13];
kiss_fft_cpx * twiddles = st->twiddles;
kiss_fft_cpx *tw;
kiss_fft_cpx ya,yb;
ya = twiddles[fstride*m];
yb = twiddles[fstride*2*m];
Fout0=Fout;
Fout1=Fout0+m;
Fout2=Fout0+2*m;
Fout3=Fout0+3*m;
Fout4=Fout0+4*m;
tw=st->twiddles;
for ( u=0; u<m; ++u ) {
C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
scratch[0] = *Fout0;
C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
C_ADD( scratch[7],scratch[1],scratch[4]);
C_SUB( scratch[10],scratch[1],scratch[4]);
C_ADD( scratch[8],scratch[2],scratch[3]);
C_SUB( scratch[9],scratch[2],scratch[3]);
Fout0->r += scratch[7].r + scratch[8].r;
Fout0->i += scratch[7].i + scratch[8].i;
scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
C_SUB(*Fout1,scratch[5],scratch[6]);
C_ADD(*Fout4,scratch[5],scratch[6]);
scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
C_ADD(*Fout2,scratch[11],scratch[12]);
C_SUB(*Fout3,scratch[11],scratch[12]);
++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
}
}
/* perform the butterfly for one stage of a mixed radix FFT */
static void kf_bfly_generic(
kiss_fft_cpx * Fout,
const size_t fstride,
const kiss_fft_cfg st,
int m,
int p
)
{
int u,k,q1,q;
kiss_fft_cpx * twiddles = st->twiddles;
kiss_fft_cpx t;
int Norig = st->nfft;
kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p);
if (scratch == NULL){
KISS_FFT_ERROR("Memory allocation failed.");
return;
}
for ( u=0; u<m; ++u ) {
k=u;
for ( q1=0 ; q1<p ; ++q1 ) {
scratch[q1] = Fout[ k ];
C_FIXDIV(scratch[q1],p);
k += m;
}
k=u;
for ( q1=0 ; q1<p ; ++q1 ) {
int twidx=0;
Fout[ k ] = scratch[0];
for (q=1;q<p;++q ) {
twidx += fstride * k;
if (twidx>=Norig) twidx-=Norig;
C_MUL(t,scratch[q] , twiddles[twidx] );
C_ADDTO( Fout[ k ] ,t);
}
k += m;
}
}
KISS_FFT_TMP_FREE(scratch);
}
static
void kf_work(
kiss_fft_cpx * Fout,
const kiss_fft_cpx * f,
const size_t fstride,
int in_stride,
int * factors,
const kiss_fft_cfg st
)
{
kiss_fft_cpx * Fout_beg=Fout;
const int p=*factors++; /* the radix */
const int m=*factors++; /* stage's fft length/p */
const kiss_fft_cpx * Fout_end = Fout + p*m;
#ifdef _OPENMP
// use openmp extensions at the
// top-level (not recursive)
if (fstride==1 && p<=5 && m!=1)
{
int k;
// execute the p different work units in different threads
# pragma omp parallel for
for (k=0;k<p;++k)
kf_work( Fout +k*m, f+ fstride*in_stride*k,fstride*p,in_stride,factors,st);
// all threads have joined by this point
switch (p) {
case 2: kf_bfly2(Fout,fstride,st,m); break;
case 3: kf_bfly3(Fout,fstride,st,m); break;
case 4: kf_bfly4(Fout,fstride,st,m); break;
case 5: kf_bfly5(Fout,fstride,st,m); break;
default: kf_bfly_generic(Fout,fstride,st,m,p); break;
}
return;
}
#endif
if (m==1) {
do{
*Fout = *f;
f += fstride*in_stride;
}while(++Fout != Fout_end );
}else{
do{
// recursive call:
// DFT of size m*p performed by doing
// p instances of smaller DFTs of size m,
// each one takes a decimated version of the input
kf_work( Fout , f, fstride*p, in_stride, factors,st);
f += fstride*in_stride;
}while( (Fout += m) != Fout_end );
}
Fout=Fout_beg;
// recombine the p smaller DFTs
switch (p) {
case 2: kf_bfly2(Fout,fstride,st,m); break;
case 3: kf_bfly3(Fout,fstride,st,m); break;
case 4: kf_bfly4(Fout,fstride,st,m); break;
case 5: kf_bfly5(Fout,fstride,st,m); break;
default: kf_bfly_generic(Fout,fstride,st,m,p); break;
}
}
/* facbuf is populated by p1,m1,p2,m2, ...
where
p[i] * m[i] = m[i-1]
m0 = n */
static
void kf_factor(int n,int * facbuf)
{
int p=4;
double floor_sqrt;
floor_sqrt = floor( sqrt((double)n) );
/*factor out powers of 4, powers of 2, then any remaining primes */
do {
while (n % p) {
switch (p) {
case 4: p = 2; break;
case 2: p = 3; break;
default: p += 2; break;
}
if (p > floor_sqrt)
p = n; /* no more factors, skip to end */
}
n /= p;
*facbuf++ = p;
*facbuf++ = n;
} while (n > 1);
}
/*
*
* User-callable function to allocate all necessary storage space for the fft.
*
* The return value is a contiguous block of memory, allocated with malloc. As such,
* It can be freed with free(), rather than a kiss_fft-specific function.
* */
kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
{
KISS_FFT_ALIGN_CHECK(mem)
kiss_fft_cfg st=NULL;
// check for overflow condition {memneeded > SIZE_MAX}.
if (nfft >= (SIZE_MAX - 2*sizeof(struct kiss_fft_state))/sizeof(kiss_fft_cpx))
return NULL;
size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fft_state)
+ sizeof(kiss_fft_cpx)*(nfft-1)); /* twiddle factors*/
if ( lenmem==NULL ) {
st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
}else{
if (mem != NULL && *lenmem >= memneeded)
st = (kiss_fft_cfg)mem;
*lenmem = memneeded;
}
if (st) {
int i;
st->nfft=nfft;
st->inverse = inverse_fft;
for (i=0;i<nfft;++i) {
const double pi=3.141592653589793238462643383279502884197169399375105820974944;
double phase = -2*pi*i / nfft;
if (st->inverse)
phase *= -1;
kf_cexp(st->twiddles+i, phase );
}
kf_factor(nfft,st->factors);
}
return st;
}
void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
{
if (fin == fout) {
//NOTE: this is not really an in-place FFT algorithm.
//It just performs an out-of-place FFT into a temp buffer
if (fout == NULL){
KISS_FFT_ERROR("fout buffer NULL.");
return;
}
kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft);
if (tmpbuf == NULL){
KISS_FFT_ERROR("Memory allocation error.");
return;
}
kf_work(tmpbuf,fin,1,in_stride, st->factors,st);
memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft);
KISS_FFT_TMP_FREE(tmpbuf);
}else{
kf_work( fout, fin, 1,in_stride, st->factors,st );
}
}
void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
{
kiss_fft_stride(cfg,fin,fout,1);
}
void kiss_fft_cleanup(void)
{
// nothing needed any more
}
int kiss_fft_next_fast_size(int n)
{
while(1) {
int m=n;
while ( (m%2) == 0 ) m/=2;
while ( (m%3) == 0 ) m/=3;
while ( (m%5) == 0 ) m/=5;
if (m<=1)
break; /* n is completely factorable by twos, threes, and fives */
n++;
}
return n;
}

View File

@ -0,0 +1,184 @@
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#ifndef KISS_FFT_H
#define KISS_FFT_H
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
// Define KISS_FFT_SHARED macro to properly export symbols
#ifdef KISS_FFT_SHARED
# ifdef _WIN32
# ifdef KISS_FFT_BUILD
# define KISS_FFT_API __declspec(dllexport)
# else
# define KISS_FFT_API __declspec(dllimport)
# endif
# else
# define KISS_FFT_API __attribute__ ((visibility ("default")))
# endif
#else
# define KISS_FFT_API
#endif
#ifdef __cplusplus
extern "C" {
#endif
/*
ATTENTION!
If you would like a :
-- a utility that will handle the caching of fft objects
-- real-only (no imaginary time component ) FFT
-- a multi-dimensional FFT
-- a command-line utility to perform ffts
-- a command-line utility to perform fast-convolution filtering
Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c
in the tools/ directory.
*/
/* User may override KISS_FFT_MALLOC and/or KISS_FFT_FREE. */
#ifdef USE_SIMD
#ifdef HAVE_LASX
# include <lasxintrin.h>
# define kiss_fft_scalar __m256
# ifndef KISS_FFT_MALLOC
# define KISS_FFT_MALLOC(nbytes) aligned_alloc(32, KISS_FFT_ALIGN_SIZE_UP(nbytes))
# define KISS_FFT_ALIGN_CHECK(ptr)
# define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 31UL) & ~0x1FUL)
# endif
# ifndef KISS_FFT_FREE
# define KISS_FFT_FREE free
# endif
#elif defined(HAVE_LSX)
# include <lsxintrin.h>
# define kiss_fft_scalar __m128
# ifndef KISS_FFT_MALLOC
# define KISS_FFT_MALLOC(nbytes) aligned_alloc(16, KISS_FFT_ALIGN_SIZE_UP(nbytes))
# define KISS_FFT_ALIGN_CHECK(ptr)
# define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL)
# endif
# ifndef KISS_FFT_FREE
# define KISS_FFT_FREE free
# endif
#else
# include <xmmintrin.h>
# define kiss_fft_scalar __m128
# ifndef KISS_FFT_MALLOC
# define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16)
# define KISS_FFT_ALIGN_CHECK(ptr)
# define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL)
# endif
# ifndef KISS_FFT_FREE
# define KISS_FFT_FREE _mm_free
# endif
#endif
#else
# define KISS_FFT_ALIGN_CHECK(ptr)
# define KISS_FFT_ALIGN_SIZE_UP(size) (size)
# ifndef KISS_FFT_MALLOC
# define KISS_FFT_MALLOC malloc
# endif
# ifndef KISS_FFT_FREE
# define KISS_FFT_FREE free
# endif
#endif
#ifdef FIXED_POINT
#include <stdint.h>
# if (FIXED_POINT == 32)
# define kiss_fft_scalar int32_t
# else
# define kiss_fft_scalar int16_t
# endif
#else
# ifndef kiss_fft_scalar
/* default is float */
# define kiss_fft_scalar float
# endif
#endif
typedef struct {
kiss_fft_scalar r;
kiss_fft_scalar i;
}kiss_fft_cpx;
typedef struct kiss_fft_state* kiss_fft_cfg;
/*
* kiss_fft_alloc
*
* Initialize a FFT (or IFFT) algorithm's cfg/state buffer.
*
* typical usage: kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL);
*
* The return value from fft_alloc is a cfg buffer used internally
* by the fft routine or NULL.
*
* If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc.
* The returned value should be free()d when done to avoid memory leaks.
*
* The state can be placed in a user supplied buffer 'mem':
* If lenmem is not NULL and mem is not NULL and *lenmem is large enough,
* then the function places the cfg in mem and the size used in *lenmem
* and returns mem.
*
* If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough),
* then the function returns NULL and places the minimum cfg
* buffer size in *lenmem.
* */
kiss_fft_cfg KISS_FFT_API kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem);
/*
* kiss_fft(cfg,in_out_buf)
*
* Perform an FFT on a complex input buffer.
* for a forward FFT,
* fin should be f[0] , f[1] , ... ,f[nfft-1]
* fout will be F[0] , F[1] , ... ,F[nfft-1]
* Note that each element is complex and can be accessed like
f[k].r and f[k].i
* */
void KISS_FFT_API kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);
/*
A more generic version of the above function. It reads its input from every Nth sample.
* */
void KISS_FFT_API kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride);
/* If kiss_fft_alloc allocated a buffer, it is one contiguous
buffer and can be simply free()d when no longer needed*/
#define kiss_fft_free KISS_FFT_FREE
/*
Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up
your compiler output to call this before you exit.
*/
void KISS_FFT_API kiss_fft_cleanup(void);
/*
* Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5)
*/
int KISS_FFT_API kiss_fft_next_fast_size(int n);
/* for real ffts, we need an even size */
#define kiss_fftr_next_fast_size_real(n) \
(kiss_fft_next_fast_size( ((n)+1)>>1)<<1)
#ifdef __cplusplus
}
#endif
#endif

View File

@ -0,0 +1,36 @@
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#ifndef kiss_fft_log_h
#define kiss_fft_log_h
#define ERROR 1
#define WARNING 2
#define INFO 3
#define DEBUG 4
#define STRINGIFY(x) #x
#define TOSTRING(x) STRINGIFY(x)
#if defined(NDEBUG)
# define KISS_FFT_LOG_MSG(severity, ...) ((void)0)
#else
# define KISS_FFT_LOG_MSG(severity, ...) \
fprintf(stderr, "[" #severity "] " __FILE__ ":" TOSTRING(__LINE__) " "); \
fprintf(stderr, __VA_ARGS__); \
fprintf(stderr, "\n")
#endif
#define KISS_FFT_ERROR(...) KISS_FFT_LOG_MSG(ERROR, __VA_ARGS__)
#define KISS_FFT_WARNING(...) KISS_FFT_LOG_MSG(WARNING, __VA_ARGS__)
#define KISS_FFT_INFO(...) KISS_FFT_LOG_MSG(INFO, __VA_ARGS__)
#define KISS_FFT_DEBUG(...) KISS_FFT_LOG_MSG(DEBUG, __VA_ARGS__)
#endif /* kiss_fft_log_h */

View File

@ -0,0 +1,171 @@
/*
* Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#include "kiss_fftr.h"
#include "_kiss_fft_guts.h"
struct kiss_fftr_state{
kiss_fft_cfg substate;
kiss_fft_cpx * tmpbuf;
kiss_fft_cpx * super_twiddles;
#ifdef USE_SIMD
void * pad;
#endif
};
kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem)
{
KISS_FFT_ALIGN_CHECK(mem)
int i;
kiss_fftr_cfg st = NULL;
size_t subsize = 0, memneeded;
if (nfft & 1) {
KISS_FFT_ERROR("Real FFT optimization must be even.");
return NULL;
}
nfft >>= 1;
kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize);
memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2);
if (lenmem == NULL) {
st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded);
} else {
if (*lenmem >= memneeded)
st = (kiss_fftr_cfg) mem;
*lenmem = memneeded;
}
if (!st)
return NULL;
st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */
st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize);
st->super_twiddles = st->tmpbuf + nfft;
kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize);
for (i = 0; i < nfft/2; ++i) {
double phase =
-3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5);
if (inverse_fft)
phase *= -1;
kf_cexp (st->super_twiddles+i,phase);
}
return st;
}
void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata)
{
/* input buffer timedata is stored row-wise */
int k,ncfft;
kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc;
if ( st->substate->inverse) {
KISS_FFT_ERROR("kiss fft usage error: improper alloc");
return;/* The caller did not call the correct function */
}
ncfft = st->substate->nfft;
/*perform the parallel fft of two real signals packed in real,imag*/
kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf );
/* The real part of the DC element of the frequency spectrum in st->tmpbuf
* contains the sum of the even-numbered elements of the input time sequence
* The imag part is the sum of the odd-numbered elements
*
* The sum of tdc.r and tdc.i is the sum of the input time sequence.
* yielding DC of input time sequence
* The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1...
* yielding Nyquist bin of input time sequence
*/
tdc.r = st->tmpbuf[0].r;
tdc.i = st->tmpbuf[0].i;
C_FIXDIV(tdc,2);
CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i);
CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i);
freqdata[0].r = tdc.r + tdc.i;
freqdata[ncfft].r = tdc.r - tdc.i;
#ifdef USE_SIMD
#ifdef HAVE_LASX
freqdata[0].i = (__m256)(__lasx_xvreplgr2vr_w(0));
freqdata[ncfft].i = freqdata[0].i;
#elif defined(HAVE_LSX)
freqdata[0].i = (__m128)(__lsx_vreplgr2vr_w(0));
freqdata[ncfft].i = freqdata[0].i;
#else
freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0);
#endif
#else
freqdata[ncfft].i = freqdata[0].i = 0;
#endif
for ( k=1;k <= ncfft/2 ; ++k ) {
fpk = st->tmpbuf[k];
fpnk.r = st->tmpbuf[ncfft-k].r;
fpnk.i = - st->tmpbuf[ncfft-k].i;
C_FIXDIV(fpk,2);
C_FIXDIV(fpnk,2);
C_ADD( f1k, fpk , fpnk );
C_SUB( f2k, fpk , fpnk );
C_MUL( tw , f2k , st->super_twiddles[k-1]);
freqdata[k].r = HALF_OF(f1k.r + tw.r);
freqdata[k].i = HALF_OF(f1k.i + tw.i);
freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r);
freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i);
}
}
void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata)
{
/* input buffer timedata is stored row-wise */
int k, ncfft;
if (st->substate->inverse == 0) {
KISS_FFT_ERROR("kiss fft usage error: improper alloc");
return;/* The caller did not call the correct function */
}
ncfft = st->substate->nfft;
st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r;
st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r;
C_FIXDIV(st->tmpbuf[0],2);
for (k = 1; k <= ncfft / 2; ++k) {
kiss_fft_cpx fk, fnkc, fek, fok, tmp;
fk = freqdata[k];
fnkc.r = freqdata[ncfft - k].r;
fnkc.i = -freqdata[ncfft - k].i;
C_FIXDIV( fk , 2 );
C_FIXDIV( fnkc , 2 );
C_ADD (fek, fk, fnkc);
C_SUB (tmp, fk, fnkc);
C_MUL (fok, tmp, st->super_twiddles[k-1]);
C_ADD (st->tmpbuf[k], fek, fok);
C_SUB (st->tmpbuf[ncfft - k], fek, fok);
#ifdef USE_SIMD
#ifdef HAVE_LASX
__m256 neg_one = (__m256)__lasx_xvreplgr2vr_w(0xBF800000); // -1.0f
st->tmpbuf[ncfft - k].i = __lasx_xvfmul_s(st->tmpbuf[ncfft - k].i, neg_one);
#elif defined(HAVE_LSX)
__m128 neg_one = (__m128)__lsx_vreplgr2vr_w(0xBF800000); // -1.0f
st->tmpbuf[ncfft - k].i = __lsx_vfmul_s(st->tmpbuf[ncfft - k].i, neg_one);
#else
st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0);
#endif
#else
st->tmpbuf[ncfft - k].i *= -1;
#endif
}
kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata);
}

View File

@ -0,0 +1,54 @@
/*
* Copyright (c) 2003-2004, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#ifndef KISS_FTR_H
#define KISS_FTR_H
#include "kiss_fft.h"
#ifdef __cplusplus
extern "C" {
#endif
/*
Real optimized version can save about 45% cpu time vs. complex fft of a real seq.
*/
typedef struct kiss_fftr_state *kiss_fftr_cfg;
kiss_fftr_cfg KISS_FFT_API kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem);
/*
nfft must be even
If you don't care to allocate space, use mem = lenmem = NULL
*/
void KISS_FFT_API kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata);
/*
input timedata has nfft scalar points
output freqdata has nfft/2+1 complex points
*/
void KISS_FFT_API kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata);
/*
input freqdata has nfft/2+1 complex points
output timedata has nfft scalar points
*/
#define kiss_fftr_free KISS_FFT_FREE
#ifdef __cplusplus
}
#endif
#endif

View File

@ -0,0 +1,22 @@
# gpr3dv-smoke CLI vendored 3DGPRViewer
# : gpr3dv-smoke <lineDir> <linePrefix>
# : gpr3dv-smoke "D:/Downloads/" "_001"
# APIIprhParser::loadImpulseMultiChannel GPRDataModel::buildVolumeData
# RadarProcessor::runPipeline(默认流水线)//
add_executable(gpr3dv_smoke main.cpp)
set_target_properties(gpr3dv_smoke PROPERTIES
OUTPUT_NAME "gpr3dv-smoke"
AUTOMOC OFF AUTOUIC OFF AUTORCC OFF)
target_link_libraries(gpr3dv_smoke PRIVATE geopro_gpr3dv Qt6::Core Qt6::Gui)
target_compile_features(gpr3dv_smoke PRIVATE cxx_std_17)
if(WIN32)
# DLLQt6Core exe 使 PATH
add_custom_command(TARGET gpr3dv_smoke POST_BUILD
COMMAND ${CMAKE_COMMAND} -E copy_if_different
$<TARGET_RUNTIME_DLLS:gpr3dv_smoke> $<TARGET_FILE_DIR:gpr3dv_smoke>
COMMAND_EXPAND_LISTS)
endif()

115
tools/gpr3dv_smoke/main.cpp Normal file
View File

@ -0,0 +1,115 @@
// gpr3dv-smoke —— vendored 3DGPRViewer 数据生成链冒烟入口。
//
// 验证目标(任务 P1 验收):
// ① 读一条多通道测线 → GPRDataModeltraces
// ② GPRDataModel::buildVolumeData → 立方体 volumeData[ch][trace][sample]
// ③ RadarProcessor::runPipeline(默认流水线) → 处理后 GPRDataModel
// 打印:立方体维度 / 处理前后统计量(证明处理生效)/ 通道数(读自数据)。
//
// 用法:
// gpr3dv-smoke <lineDir> <linePrefix>
// 例: gpr3dv-smoke "D:/Downloads/明星路" "明星路_001"
//
// 算法零改动:本文件只调用 vendored 原版静态 API不触碰任何算法体。
#include <QString>
#include <QStringList>
#include <QtGlobal>
#include <cstdio>
#include "GPRDataModel.h"
#include "IprhParser.h"
#include "RadarProcessor.h"
// 全体 traces 的平均绝对幅值——朴素的"处理是否生效"标量证据。
// RadarProcessor 流水线在 model.traces 上做原位变换(非 volumeData故统计取自 traces。
static double meanAbsAmplitude(const GPRDataModel &model, qint64 *sampleCountOut = nullptr)
{
long double sum = 0.0L;
qint64 count = 0;
for (const RadarTrace &tr : model.traces) {
for (short v : tr.amplitudes) {
sum += static_cast<long double>(v < 0 ? -v : v);
++count;
}
}
if (sampleCountOut) {
*sampleCountOut = count;
}
return count > 0 ? static_cast<double>(sum / count) : 0.0;
}
int main(int argc, char **argv)
{
if (argc < 3) {
std::fprintf(stderr,
"用法: gpr3dv-smoke <lineDir> <linePrefix>\n"
"例: gpr3dv-smoke \"D:/Downloads/明星路\" \"明星路_001\"\n");
return 2;
}
const QString lineDir = QString::fromLocal8Bit(argv[1]);
const QString linePrefix = QString::fromLocal8Bit(argv[2]);
std::printf("[gpr3dv-smoke] lineDir=%s linePrefix=%s\n",
lineDir.toLocal8Bit().constData(),
linePrefix.toLocal8Bit().constData());
// ---- ① 读多通道测线(通道数读自数据,不写死)----
GPRDataModel model;
if (!IprhParser::loadImpulseMultiChannel(lineDir, linePrefix, model)) {
std::fprintf(stderr, "[gpr3dv-smoke] 错误: loadImpulseMultiChannel 失败\n");
return 1;
}
const int channels = model.header.numberOfChannels; // 读自 .iprh CHANNELS / _A 文件数
const int tracesPerChannel = model.getTracesPerChannel();
const int samplesPerTrace = model.header.samplesPerTrace;
std::printf("[gpr3dv-smoke] 加载完成: 通道数=%d 每通道道数=%d 每道样本数=%d 总道数=%lld\n",
channels, tracesPerChannel, samplesPerTrace,
static_cast<long long>(model.traces.size()));
// ---- ② 建三维立方体 volumeData[ch][trace][sample] ----
model.buildVolumeData();
const int volCh = model.volumeData.size();
const int volTr = (volCh > 0) ? model.volumeData[0].size() : 0;
const int volSm = (volCh > 0 && volTr > 0) ? model.volumeData[0][0].size() : 0;
std::printf("[gpr3dv-smoke] 立方体维度 [通道 x 道 x 样本] = %d x %d x %d\n",
volCh, volTr, volSm);
short gMin = 0, gMax = 0;
model.getGlobalMinMax(gMin, gMax);
std::printf("[gpr3dv-smoke] 立方体全局幅值范围: min=%d max=%d\n",
static_cast<int>(gMin), static_cast<int>(gMax));
// ---- 处理前统计 ----
qint64 nBefore = 0;
const double meanBefore = meanAbsAmplitude(model, &nBefore);
std::printf("[gpr3dv-smoke] 处理前 平均绝对幅值=%.4f (样本数=%lld)\n",
meanBefore, static_cast<long long>(nBefore));
// ---- ③ RadarProcessor 默认流水线 ----
RadarProcessor::ProcPipeline pipeline;
pipeline.setDefaultFlow();
std::printf("[gpr3dv-smoke] 运行默认流水线(步骤数=%lld...\n",
static_cast<long long>(pipeline.steps.size()));
const GPRDataModel processed = RadarProcessor::runPipeline(model, pipeline);
qint64 nAfter = 0;
const double meanAfter = meanAbsAmplitude(processed, &nAfter);
std::printf("[gpr3dv-smoke] 处理后 平均绝对幅值=%.4f (样本数=%lld 道数=%lld 样本/道=%d)\n",
meanAfter, static_cast<long long>(nAfter),
static_cast<long long>(processed.traces.size()),
processed.header.samplesPerTrace);
const double delta = meanAfter - meanBefore;
std::printf("[gpr3dv-smoke] 处理前后差值=%.4f —— %s\n",
delta,
(qAbs(delta) > 1e-9 ? "处理已生效(统计量变化)" : "警告:统计量未变化"));
std::printf("[gpr3dv-smoke] OK\n");
return 0;
}