信号特征之相空间重构（Python、C++、MATLAB实现）

最新推荐文章于 2024-07-26 21:39:10 发布

Vistyxio

最新推荐文章于 2024-07-26 21:39:10 发布

阅读量2.7k

点赞数 31

文章标签： python c++ 算法 matlab 信息与通信信号处理

本文链接：https://blog.csdn.net/weixin_45715390/article/details/136527188

版权

本文介绍了相空间重构的基本概念，包括利用序列数据构建相空间轨迹的方法，着重于时间延迟参数τ的选择策略，如通过自相关函数确定，以及如何在Python、C++和MATLAB中实现相空间重构函数。通过二维轨迹图展示信号特征的研究实例。

摘要由CSDN通过智能技术生成

相空间重构

1 特征描述

相空间重构的基本思想为：系统中任意变量的变化过程由与之相互作用的分量共同决定，其发展过程也隐含着其他变量的变化规律，故可从变量的变化过程中构建和恢复整个系统的变化规律。针对序列 $x_1,x_2,\cdots,x_N$ ，通过坐标延迟重构法，可以将序列映射到一个 $d$ 维的相空间中，重构后的轨迹矩阵如下。

$\begin{bmatrix} x_1 & x_{1+\tau} & \cdots & x_{1+(d-1)\tau} \\ x_2 & x_{2+\tau} & \cdots & x_{2+(d-1)\tau} \\ \vdots & \vdots & & \vdots \\ x_P & x_{P+\tau} & \cdots & x_{P+(d-1)\tau} \end{bmatrix}$

上式中， $\tau$ 为时间延迟参数； $d$ 为嵌入维数，每一行可视为相空间中的一个相点的坐标；重构的相空间轨迹由 $P$ 个相点组成。 $P$ 的大小取决于 $d$ 和 $\tau$ ， $P=N-(d-1)\tau$ 。

相空间重构的核心问题在于如何选择合适的时间延迟参数 $\tau$ 和嵌入维数 $d$ 。 $\tau$ 太小会使相点分布过于密集，太大则过于分散； $d$ 过大会增大计算量，过小则难以显示隐藏特征。

一般嵌入维数 $d$ 可根据具体课题选定，例如使用二维轨迹图来提取特征进行研究，则直接选择 $d$ 为 $2$ 。而时间延迟参数 $\tau$ 可以使用自相关函数进行计算，当自相关函数值达到第一次极小值或近似为零时，此时的 $\tau$ 比较合适。若有多个信号序列，为简化计算，可取所有序列的合适 $\tau$ 值的平均值作为设定的时间延迟参数 $\tau$ 。定义自相关函数如下。

$C(\tau)=\frac{\sum\limits_{i=1}^{n-\tau}(x_i-\overline{x})(x_{i+\tau}-\overline{x})}{\sum\limits_{i=1}^{n-\tau}(x_i-\overline{x})^2}$

其中， $\overline{x}$ 为时间序列的均值。

2 数据来源

提供一串信号的数据:
这组数据存在三个频率分量，分别为100Hz，200Hz和300Hz, 用MATLAB生成。

# 定义三个频率分量的参数
t = 0:0.0001:2/100;

freq1 = 100; # 第一个频率分量的频率(Hz)
freq2 = 200; # 第二个频率分量的频率(Hz)
freq3 = 300; # 第三个频率分量的频率(Hz)

# 生成三个频率分量的正弦波信号
signal1 = sin(2 * pi * freq1 * t);
signal2 = sin(2 * pi * freq2 * t);
signal3 = sin(2 * pi * freq3 * t);

# 将三个频率分量的信号相加
result = signal1 + signal2 + signal3;

3 函数的Python代码

3.1 Python库要求

import numpy as np

3.2 PhaseSpaceReconstruction(inputS, dd, tt)相空间重构函数

其中，inputS为输入的信号序列，dd为嵌入维度，tt为延迟参数。当tt为-1时，相空间重构函数会按照自相关函数计算合适的延迟参数，若tt为任意正整数，则相空间重构函数会根据用户输入的延迟参数进行相空间重构。函数的返回值matrixX为numpy类型的轨迹矩阵，每一行是一个相空间中的相点。

def PhaseSpaceReconstruction(inputS, dd, tt):
    data = []
    for i in inputS:
        data.append(float(i))
    length = len(data)
    xave = np.mean(data)
    Cup, Cdown, Clast = 0x3f3f3f3f, 1, 0
    besttau = 0
    # 计算合适的tau
    if tt == -1:
        for tau in range(1, length):
            Clast = (Cup / Cdown)
            Cup, Cdown = 0, 0
            for i in range(0, length - tau):
                Cup = Cup + (data[i] - xave) * (data[i + tau] - xave)
                Cdown = Cdown + (data[i] - xave) * (data[i] - xave)
            if (Cup / Cdown) > Clast:
                besttau = tau - 1
                break
    else:
        besttau = tt
    P = length - besttau * (dd - 1)
    maxdata, mindata = max(data), min(data)
    # 归一化
    for i in range(0, length):
        data[i] = (data[i] - mindata) / (maxdata - mindata)

    matrixX = np.zeros((P, dd))
    for i in range(0, P):
        for j in range(0, dd):
            matrixX[i, j] = data[i + besttau * j]
    
    return matrixX

3.3 相空间重构函数PhaseSpaceReconstruction(inputS, dd, tt)验证

dd = 2
inputS = input().split()


matrixX = PhaseSpaceReconstruction(inputS, dd, 10)


for i in range(0, matrixX.shape[1]):
    for j in matrixX[:,i]:
        print("{} ".format(j), end='')
    print()
    print()

4 函数的C++代码

4.1 C++库要求

#include<iostream>
#include<fstream>
#include<vector>
#include<algorithm>

4.2 void PhaseSpaceReconstruction(vector<double> sequenceX, vector<vector<double>> &matrixX, int d, int tau)相空间重构函数

其中，sequenceX为输入的信号序列，d为嵌入维度，tau为延迟参数。当tau为-1时，相空间重构函数会按照自相关函数计算合适的延迟参数，若tau为任意正整数，则相空间重构函数会根据用户输入的延迟参数进行相空间重构。matrixX为vector<vector<double>>类型的轨迹矩阵，每一行是一个相空间中的相点。

void PhaseSpaceReconstruction(vector<double> sequenceX, vector<vector<double>> &matrixX, int d, int tau) {

    int length = sequenceX.size();

    if (tau == -1) {
        auto besttau = [length](vector<double> sequenceX) -> int {
            double xsum = 0., xave = 0.;
            double Cup = double(0x3f3f3f3f), Cdown = 1., Clast = 0.;
            for (auto i : sequenceX) xsum += i;
            xave = xsum / double(length);
            for (size_t tau = 1; tau < length; tau++) {
                Clast = Cup / Cdown;
                Cup = 0., Cdown = 0.;
                for (size_t i = 0; i < (length - tau); i++) {
                    Cup = Cup + (sequenceX[i] - xave) * (sequenceX[i + tau] - xave);
                    Cdown = Cdown + (sequenceX[i] - xave) * (sequenceX[i] - xave);
                }
                if ((Cup / Cdown) > Clast) {
                    return (tau - 1);
                }
            }
            return -1;
        };

        tau = besttau(sequenceX);
    }

    int P = length - tau * (d - 1);

    double maxX = *max_element(sequenceX.begin(), sequenceX.end());
    double minX = *min_element(sequenceX.begin(), sequenceX.end());
    for (auto & i : sequenceX) i = (i - minX) / (maxX - minX);

    vector<double> vtemp;

    for (size_t i = 0; i < P; i++) {
        for (size_t j = 0; j < d; j++) {
            vtemp.emplace_back(sequenceX[i + tau * j]);
        }
        matrixX.emplace_back(vtemp);
        vtemp.clear();
    }

    return;
}

4.3 相空间重构函数void PhaseSpaceReconstruction(vector<double> sequenceX, vector<vector<double>> &matrixX, int d, int tau)验证

int main() {
    
    vector<double> sequenceX; //N个
    vector<vector<double>> matrixX; //P行d列

    int NN, dd, tt;

    NN = 201; //数据个数
    dd = 2; //嵌入维数
    tt = 10; //时间延迟参数

    double temp;
    fstream fs;

    //从文件读
    fs.open("PhaseSpaceReconstructionINPUTDATA.txt", ios::in);
    if (!fs.is_open()) {
        cout << "no txt file" << endl;
        return -1;
    }
    while (fs >> temp) {
        sequenceX.emplace_back(temp);
    }
    fs.close();

    PhaseSpaceReconstruction(sequenceX, matrixX, dd, tt);

    //输出至文件
    fs.open("PhaseSpaceReconstructionOUTPUTDATA.txt", ios::out);
    for (size_t i = 0; i < dd; i++) {
        for (auto j : matrixX) fs << j[i] << " ";
        fs << endl << endl;
    }
    fs.close();

    return 0;
}

5 函数的MATLAB代码

5.1 PhaseSpaceReconstruction相空间重构函数

其中，result为输入的信号序列，d为嵌入维度，tau为延迟参数。当tau为-1时，相空间重构函数会按照自相关函数计算合适的延迟参数，若tau为任意正整数，则相空间重构函数会根据用户输入的延迟参数进行相空间重构。函数的返回值matrixX为相空间重构的轨迹矩阵，每一行是一个相空间中的相点。

function matrixX = PhaseSpaceReconstruction(result, d, tau)

% 相空间重构
% 
% matrixX为输出相点矩阵
% result为信号序列
% d为嵌入维数
% tau为时间延迟参数
% tau为-1时由程序计算合适的tau

Cup = 1061109567;
Cdown = 1;
Clast = 0;

len = length(result);
xave = mean(result);

if tau == -1
    for tau = 1 : (len - 1)
        Clast = Cup / Cdown;
        Cup = 0;
        Cdown = 0;
        for i = 1 : (len - tau)
            Cup = Cup + (result(i) - xave) * (result(i + tau) - xave);
            Cdown = Cdown + (result(i) - xave) * (result(i) - xave);
        end
        if (Cup / Cdown) > Clast
            tau = tau - 1;
            break
        end
    end
end

maxR = max(result);
minR = min(result);

resultONE = (result - minR) / (maxR - minR);

P = len - tau * (d - 1);

matrixX = zeros(P, d);

for i = 1 : P
    for j = 1 : d
        matrixX(i, j) = resultONE(i + tau * (j - 1));
    end
end

end

5.2 相空间重构函数PhaseSpaceReconstruction验证

result = \\* 输入样本数据 *\\;
matrixX = PhaseSpaceReconstruction(result, 2, 10);
plot(matrixX(:, 1), matrixX(:, 2))

6 结果

6.1 Python绘制的相空间重构图像

在这里插入图片描述

6.2 C++绘制的相空间重构图像

在这里插入图片描述

6.3 MATLAB绘制的相空间重构图像

在这里插入图片描述

Vistyxio

关注

31
点赞
踩
48

收藏

觉得还不错? 一键收藏
3
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫