Canmv k230 案例3.1——短时傅里叶STFT（一）

短时傅里叶STFT

注：本文可能不严谨，主要是对STFT的简单阐述
Canmv k230中集成了可调用的硬件FFT，但使用过程中存在问题，可以采用ulab中提供的FFT函数

当FFT沿着时间轴滑动得到的频率随时间变化的矩阵，就可以得到信号时频的变化——STFT。

MATLAB中对于STFT的实现方法，自带函数stft，下面对调频信号测试

t = 0:1/1e3:2;
y = chirp(t,0,1,250);
plot(y)
fs = 10e3;
pspectrum(y,1e3,'spectrogram','TimeResolution',0.1, ...
    'OverlapPercent',99,'Leakage',0.85)

stft(y,fs,Window=hamming(64),OverlapLength=32,FFTLength=128);

500Hz内为随时间增加频率增加，大于500Hz为对称的频谱，matlab中对stft的绘图描述如下

需要先对信号进行加窗处理，后续再说明

canmv ide 代码

# 基础示例 # STFT 算法 by 花月
from machine import FFT
import array
import math
from ulab import numpy as np
import time
PI = 3.14159265358979323846264338327950288419716939937510
N=1024 # 输入数据长度
w1=0   # 初始频率
w2=PI  # 停止频率
x = [] # 输入信号
M=64   # 窗的长度
# 产生线性调频信号
def input_data():
    for n in range(N):
        data0 = math.cos(w1 * n + ( w2 - w1 ) * n * n /(2 * (N - 1)))
        x.append(data0)
#初始化需要进行FFT的数据，列表类型
input_data()                                                            
print('输入数据：\n',x)
# hamming 窗
v=[]
def hamming():
    for m in range(M):
        data0 = 0.54 - 0.46 * math.cos(2 * PI * m/M)
        v.append(data0)
hamming()
print('hamming 窗：\n',v)

# 做STFT
STFT1=[]
# 定义为函数
def STFT():
    for n in range(0,N,64):
        STFT=[]
        for m in range(M):
            vx=v[m] * x[n+m]
            STFT.append(vx)

        data = np.array(STFT,dtype=np.float)
        a, b =np.fft.fft(data)
        FA1=[]
        for i in range(M/2):   # 一半的输出为对称结构
            FA=math.pow(a[i], 2)+math.pow(b[i], 2)
            FA=math.sqrt(FA)/32
            FA1.append(FA)

        STFT1.append(FA1)

STFT()
print(STFT1)
STFT_output = np.array(STFT1,dtype=np.float)
print(STFT_output.shape)

将ide中的输出数据复制到MATLAB中

B=reshape(ans,[32,16])

输出结果（mesh命令）
输出结果（surf命令），旋转三维图至俯视图

后续
1)STFT具体参数
2)ADC+STFT