Canmv k230 案例6.2—— S变换代码Python版

S变换(Stockwell Transform)代码

由于S变换的原理解释需要花费更多的时间，而python代码列出来又会过于看不懂，但还是先放出来，总比一直不做的要好。

# 实现ST变换- By: 花月mmc
# 初步实现 基本无错误，误差较小

import array  # 用于转换成数组
import math   # 提供一些基本的运算
from ulab import numpy as np # 这是一个支持部分numpy功能的库
# 定义PI常数
PI = 3.14159265358979323846264338327950288419716939937510
# 这是一个要生成的数据，待测试数据
rx = []
def input_data():
    # 这是一个基波
    for i in range(8):
        data0 = 1 * math.sin(2 * 1 * PI * (i+1) / 8)
        rx.append((float(data0)))
input_data()

rx1=[]
# 这是一个三次谐波
for i in range(8):
    data1 = 1 * math.sin(2 * 3 * PI * (i+1) / 8)
    rx1.append((float(data1)))
# 让数据的后半部分是存在谐波，频谱上存在两种频率
rx[4:8]=rx1[4:8]
print('列表输出',rx)
x=np.array(rx,dtype=np.complex) # 待处理数据
print('复数输出\n',x)


# 信号长度
N=len(rx)
# 信号长度的一半，用于频率计算
nhaf=int(N/2)
# 这是一个频率索引
f=[]
# 前半部分，符合奈奎斯特采样频率，采样率的一半
f1=list(range(nhaf+1))
# 后半部分是负数，即表示对称结构的
f2=list(range(-nhaf+1,0))
for i in range(len(f1)):
    f.append(f1[i]/N)
for i in range(len(f2)):
    f.append(f2[i]/N)

#print('频率分割\n',f)

# 执行DFT 前面的文章由分享
def DFT(x):
    Hw=[] # Hw H(w) 频域分布
    xl=len(x)
    for m in range(xl):
        ww=[]
        for n in range(xl):
            W=x[n]*np.exp(-1j*2*math.pi*n*m/xl)
            ww.append(W)
        Hw.append(sum(ww))
    Hw=np.array(Hw,dtype=np.complex)
    return Hw
# 调用 DFT
# 计算频率分布的倒数
invfk=[]
for i in range(nhaf):
    invfk.append(1/f[i+1])

#print('1/频率分割\n',invfk)
# 按行生成nhaf个频率矩阵 可以预定义

f_nhaf=np.zeros((nhaf, N))
#print(f_nhaf)

for i in range(nhaf):
    f_nhaf[i,:]=f

#print('repmat(f,nhaf,1)\n',f_nhaf)
# 按行生成 1/f个频率的倒数矩阵 可以预定义
invfk_N=np.zeros((nhaf, N))
#print(f_nhaf)
for i in range(N):
    invfk_N[:,i]=invfk

#print('repmat(invfk,1,N);\n',invfk_N)
W=np.zeros((nhaf, N))
# 它们是构建G(m,n)矩阵的基本单元
# 定义两个矩阵
for i in range(nhaf):
    for j in range(N):
        W[i,j]=2*PI*f_nhaf[i,j]*invfk_N[i,j]

#print('高斯窗系数\n',W)
# G=exp((-W.^2)/2);
G=np.zeros((nhaf, N),dtype=np.float)  # 和复数计算不同哦
for i in range(nhaf):
    for j in range(N):
        G[i,j]=np.exp(-(W[i,j]*W[i,j])/2)

#print('高斯窗G(m,n)\n',G)
# 以上的都是预定义的内容


# 开始计算 ST
Hft=DFT(x)
#print('DFT 输出 Hft \t\t:\n',Hft)

HW=np.zeros((nhaf+1, N),dtype=np.complex)  # 和复数计算不同哦
#print('频域H(m,n)\n',HW)

for i in range(nhaf+1):
    for j in range(N):
        if(i-j)<0:
            HW[i,j]=Hft[abs(i-j)]
        else:
            HW[i,j]=np.conjugate(Hft[abs(i-j)])

#print('HW 5*8 输出矩阵\n',HW)

HW1=np.zeros((nhaf, N),dtype=np.complex)
for i in range(nhaf):
    HW1[i,:]=HW[i+1,:]

#print('HW1 4*8 输出矩阵\n',HW1)

# 计算信号与高斯窗的卷积？
HW1G=np.zeros((nhaf, N),dtype=np.complex)
for i in range(nhaf):
    for j in range(N):
        HW1G[i,j]=HW1[i,j]*G[i,j]

#print('HW1G 4*8 输出矩阵\n',HW1G)
#print(abs(HW1G))


# IDFT 傅里叶的逆变换
def IDFT(Hw):
    Hw=np.array(Hw,dtype=np.complex)
    x=[]
    xl=len(Hw)
    for m in range(xl):
        ww=[]
        for n in range(xl):
            W=Hw[n]*np.exp(1j*2*math.pi*n*m/xl)
            ww.append(W)
        x.append(sum(ww)/xl)
    x=np.array(x,dtype=np.complex)
    return x
# 调用 IDFT
#x=IDFT(Hft)
#print('IDFT 输出 x \t\t:\n',x.real)

# 需要进行按行逆变换输出S变换矩阵
ST=np.zeros((nhaf, N),dtype=np.complex)
for i in range(nhaf):
    ST[i,:]=IDFT(HW1G[i,:])

print('ST 中间 输出\n',ST)
# 计算直流分量
ST0=np.ones(N,dtype=np.float)*np.mean(np.array(rx,dtype=np.float))  # 和复数计算不同哦
print('均值初始化ST0\n',ST0)

# 此处应为矩阵拼接
# 顺序问题
#ST0=np.array(ST0,dtype=np.complex) # 转换为复数后数据变为0
print('查看ST0类型\n',ST0)
ST0=ST0.tolist()
print('ST0转换为list\n',ST0)
ST0=np.array(ST0,dtype=np.complex)
print('ST0转换为复数\n',ST0)
print('赋值转换\n',ST)
ST[1:nhaf,:]=ST[0:nhaf-1,:]
ST[0,:]=ST0

print('赋值前的ST[0,:]转换\n',ST[0,:])

# 复数的实部和虚部不支持迭代
# 通过转换为列表再转换为数据实现赋值

print('赋值后的ST[0,:]\n',ST[0,:])
#print('ST 最终 输出  \t\t:\n',ST)
print('ST 最终 输出  \t\t:\n',ST)

print('ST 最终 输出ABS幅值  \t\t:\n',abs(ST))

以后后完善，到终版可能还有几版

Canmv k230

实测待续。。。