Canmv k230 案例2.2——FFT测试（二）

最近关注的人在增加，因此快速更新一期
1）现有FFT的问题
2）FFT的一些基本概念
3）计算FFT的结果

首先看下Canmv IDE 中FFT的函数部分及输出结果，原程序可能是为了演示效果，因此以下程序进行了部分修改：1、修改数组转换为int16，硬件FFT支持的数据类型；2、补偿了原始FFT存在系数上的差异，将幅值除以2

# 基础示例
#
# 欢迎使用CanMV IDE, 点击IDE左下角的绿色按钮开始执行脚本

from machine import FFT
import array
import math
from ulab import numpy as np
PI = 3.14159265358979323846264338327950288419716939937510

rx = []   # 空列表，用于储存生成的输入数据
def input_data():    # 输入列表，5种频率
    for i in range(64):
        data0 = 10 * math.cos(2 * PI * i / 64)
        data1  = 20 * math.cos(2 * 2 * PI * i / 64)
        data2  = 30 * math.cos(3 * 2 * PI * i / 64)
        data3  = 0.2 * math.cos(4 * 2 * PI * i / 64)
        data4  = 1000 * math.cos(5 * 2 * PI * i / 64)
        rx.append((int(data0 + data1 + data2 + data3 + data4)))
input_data()                                                            #初始化需要进行FFT的数据，列表类型
#print(rx)#
data = np.array(rx,dtype=np.int16)                                     #把列表数据转换成数组，有符号整型
#print(data)
fft1 = FFT(data, 64, 0)                                                 #创建一个FFT对象,运算点数为64，偏移是0
res = fft1.run()                                                        #获取FFT转换后的数据
print(res)

res = fft1.amplitude(res)                                              #获取各个频率点的幅值
print(res)
print(res[1]/2,res[2]/2,res[3]/2,res[4]/2,res[5]/2)                    #为了串口输出
resA = np.array(res,dtype=np.float)/2   # 转换成数据并补偿幅值参数除以2
print(resA)

res = fft1.freq(64,64)                                               #获取所有频率点的频率值
print(res)

输出的计算结果如下所示

#FFT的复数结算结果
[(5, 0), (318, 3), (642, 0), (958, -1), (6, 0), (31980, -1), (-2, -1), (0, 1), (3, 0), (4, 1), (-2, 0), (4, 1), (1, 0), (5, 1), (0, 1), (3, 2), (5, 0), (1, 0), (2, 1), (-1, 3), (3, 0), (0, 1), (-2, 0), (-2, -1), (3, 0), (-6, -1), (-6, -1), (-4, -1), (2, 0), (2, 1), (0, 0), (2, 1)]
#FFT的复数结算结果
[0, 19, 40, 59, 0, 1998, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
# 幅值/2的结果，因为原始数据为10、20、30、0.2、1000
9.5 20.0 29.5 0.0 999.0
# 幅值数组形式，串口数据显示不全，因此采用上述的列表形式
array([0.0, 9.5, 20.0, ..., 0.0, 0.0, 0.0], dtype=float32)
# 每个点所对应的频率
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]

程序中存在一段代码，具体内容可以参阅信号处理相关书籍

res = fft1.freq(64,64)   # fs/N 为FFT的最小分辨率

从计算频率的结果可以看出，因此计算采用int16，因此在幅值计算时可能出现截断，幅值为太小则无法就显示正确结果，因此调整输入数据的范围可以得到接近理论值的结果。
修改输入数据幅值均为为10，结果如下

9.5 9.5 10.0 9.5 9.5

看起来精度还是不太高，那么如果自定义FFT的结果怎样呢？
幸运地是，ulab中numpy中提供fft

# 方法1 k230 FFT
fft1 = FFT(data, 64, 0)                                                 #创建一个FFT对象,运算点数为64，偏移是0
res = fft1.run()                                                        #获取FFT转换后的数据
res = fft1.amplitude(res)
print('fft1.run() :',res[0:10])
resA = np.array(res,dtype=np.float)/2
print('fft1.run()/2 :',resA[0:10])

# 方法2 ulab的numpy
a, b =np.fft.fft(data)  # a实部 b虚部
# 求幅值
FA1=[]
for i in range(10):
    FA=math.pow(a[i], 2)+math.pow(b[i], 2)
    FA=math.sqrt(FA)/32
    FA1.append(FA)
print('fft.fft :',FA1)

结果

fft1.run() : [0, 19, 19, 20, 19, 19, 0, 0, 0, 0]
fft1.run()/2 : array([0.0, 9.5, 9.5, 10.0, 9.5, 9.5, 0.0, 0.0, 0.0, 0.0], dtype=float32)
fft.fft : [0.4375, 9.661581, 9.865005, 10.07208, 9.709345, 9.738851, 0.09873316, 0.1562796, 0.08838835, 0.02080378]

计算结果相近浮点更接近实际值

下期可能更新FFT的一些实际应用问题、与ADC结合的操作