加速度积分得到速度的Omega算法

       在工程中,有时候很难直接测量机械设备部件的速度或位移,特别是在运动着的设备上。即使在一般的场景下可以布置速度传感器,但是有时候需要同时监测速度和加速度,为了节约成本,往往通过将加速度信号直接进行积分得到速度信号。

       将加速度信号积分得到速度信号有两种方法,一是直接在时域进行积分,这往往会产生趋势项,积分的结果需要进行去趋势处理。二是在频域进行积分,典型的频域积分算法是omega算法,这里做一下简单的介绍:

        如果我们已经获得了加速度信号x''(t) , 其傅里叶变换为X''(f), 那么有:

                                    20200818164254175.png

        记速度信号为x'(t),这是我们所要求的,设其傅里叶变换为X'(f), 则有:

                                     20200818164509496.png

        因为加速度可以由速度求导得到,所以他们之间的关系是:

                                        20200818164644835.png

        把速度的表达式带入上式,有:

                                       20200818164820159.png

        比较上式与第一个式子,可以得出:

                                              20200818164952768.png

        因为时域加速度信号是传感器获取到的,可以通过计算得到其频域表示,进而通过计算得到速度的频域表示,即为:

                                             20200818165257245.png

        因此再通过傅里叶逆变换就可以得到速度的时域波形。这里要注意的是,在计算速度的频域表示时,在f=0时刻,上述公式无法计算,因此直接令X'(0) = 0 + 0j ,这样得到的速度信号没有了直流值,因此无需去趋势就可以得到想要的速度波形。下面附上我在matlab上通过时域积分(最小二乘算法去趋势后)和频域积分方法得到的速度波形:

              watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0QzNDQ4OTYyMjQ=,size_16,color_FFFFFF,t_70

  可以看到,基本上是一致的。

注:源代码有偿获取。

 

 

 

 

 

 

  • 14
    点赞
  • 80
    收藏
    觉得还不错? 一键收藏
  • 51
    评论
以下是一个基于三轴加速度数据的计步算法示例,使用了低通滤波、华东窗口和波峰波谷检测。请注意,这是一个简单的示例,实际应用中可能需要根据具体情况进行调整和优化。 ```c #include <stdio.h> #include <math.h> #define SAMPLE_RATE 50 // 采样率,单位 Hz #define WINDOW_SIZE 20 // 窗口大小,单位采样点 #define THRESHOLD_FACTOR 0.25 // 波峰波谷阈值系数 #define MIN_STEP_INTERVAL 300 // 最小步间隔,单位毫秒 // 低通滤波器 typedef struct { double a1, a2, b0, b1, b2; double x1, x2, y1, y2; } LowPassFilter; void initLowPassFilter(LowPassFilter* filter, double cutoffFreq) { double omega = 2.0 * M_PI * cutoffFreq / SAMPLE_RATE; double alpha = sin(omega) / (2.0 * 0.707); double beta = sqrt(2.0) * alpha; double a0 = 1.0 + beta + alpha * alpha; filter->a1 = (2.0 * (alpha * alpha - 1.0)) / a0; filter->a2 = (1.0 - beta + alpha * alpha) / a0; filter->b0 = (alpha * alpha) / a0; filter->b1 = (2.0 * alpha * alpha) / a0; filter->b2 = (alpha * alpha) / a0; filter->x1 = 0.0; filter->x2 = 0.0; filter->y1 = 0.0; filter->y2 = 0.0; } double lowPassFilter(LowPassFilter* filter, double x) { double y = filter->b0 * x + filter->b1 * filter->x1 + filter->b2 * filter->x2 - filter->a1 * filter->y1 - filter->a2 * filter->y2; filter->x2 = filter->x1; filter->x1 = x; filter->y2 = filter->y1; filter->y1 = y; return y; } // 华东窗口 typedef struct { double* buffer; int size; int head; double sum; } MovingWindow; void initMovingWindow(MovingWindow* window, int size) { window->buffer = (double*) malloc(size * sizeof(double)); window->size = size; window->head = 0; window->sum = 0.0; } void destroyMovingWindow(MovingWindow* window) { free(window->buffer); } void addToMovingWindow(MovingWindow* window, double x) { window->sum += x; window->sum -= window->buffer[window->head]; window->buffer[window->head] = x; window->head = (window->head + 1) % window->size; } double getMovingWindowAvg(MovingWindow* window) { return window->sum / window->size; } // 波峰波谷检测状态 typedef enum { STATE_IDLE, STATE_RISING, STATE_FALLING } PeakDetectionState; // 波峰波谷检测 typedef struct { PeakDetectionState state; double threshold; double lastValue; double lastPeak; double lastValley; } PeakDetector; void initPeakDetector(PeakDetector* detector, double threshold) { detector->state = STATE_IDLE; detector->threshold = threshold; detector->lastValue = 0.0; detector->lastPeak = 0.0; detector->lastValley = 0.0; } int updatePeakDetector(PeakDetector* detector, double value, double* peak, double* valley) { switch (detector->state) { case STATE_IDLE: if (value > detector->lastValue + detector->threshold) { detector->state = STATE_RISING; detector->lastPeak = value; *peak = value; return 1; } else if (value < detector->lastValue - detector->threshold) { detector->state = STATE_FALLING; detector->lastValley = value; *valley = value; return 1; } break; case STATE_RISING: if (value < detector->lastPeak) { detector->state = STATE_FALLING; detector->lastValley = value; *valley = value; return 1; } break; case STATE_FALLING: if (value > detector->lastValley) { detector->state = STATE_RISING; detector->lastPeak = value; *peak = value; return 1; } break; } detector->lastValue = value; return 0; } // 计步器 typedef struct { LowPassFilter filter; MovingWindow window; PeakDetector peakDetector; int stepCount; int lastStepTime; } StepCounter; void initStepCounter(StepCounter* counter) { initLowPassFilter(&counter->filter, 5.0); initMovingWindow(&counter->window, WINDOW_SIZE); initPeakDetector(&counter->peakDetector, THRESHOLD_FACTOR * getMovingWindowAvg(&counter->window)); counter->stepCount = 0; counter->lastStepTime = 0; } void destroyStepCounter(StepCounter* counter) { destroyMovingWindow(&counter->window); } int updateStepCounter(StepCounter* counter, double x, double y, double z) { double norm = sqrt(x * x + y * y + z * z); double filtered = lowPassFilter(&counter->filter, norm); addToMovingWindow(&counter->window, filtered); double avg = getMovingWindowAvg(&counter->window); double peak = 0.0, valley = 0.0; if (updatePeakDetector(&counter->peakDetector, filtered - avg, &peak, &valley)) { int interval = SAMPLE_RATE * (peak - valley) / norm; if (interval > MIN_STEP_INTERVAL && counter->lastStepTime + interval < SAMPLE_RATE * 2) { counter->stepCount++; counter->lastStepTime += interval; return 1; } } return 0; } int main() { // 假设有一组加速度数据 ax, ay, az,每秒50个样本 double ax[1000], ay[1000], az[1000]; // 初始化计步器 StepCounter counter; initStepCounter(&counter); // 处理加速度数据 for (int i = 0; i < 1000; i++) { int step = updateStepCounter(&counter, ax[i], ay[i], az[i]); if (step) { printf("Step %d detected\n", counter.stepCount); } } // 销毁计步器 destroyStepCounter(&counter); return 0; } ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 51
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值