【MicroPython学习笔记】02 添加自定义模块（使用C语言实现）

Fyra-BH

已于 2024-04-21 01:09:20 修改

阅读量3.7k

点赞数 6

文章标签： python c语言嵌入式 fft

于 2021-03-05 12:28:13 首次发布

本文链接：https://blog.csdn.net/Slade99X/article/details/114384004

版权

系列文章目录

【MicroPython学习笔记】01 使用块设备

前言

在使用MicroPython的过程中，我们可以通过import导入各种各样的库，非常的方便。但有的时候我们想要使用的功能在这些库中找不到，又或者是某些功能因MicroPython的性能受限，这个时候就需要为MicroPython固件添加自定义模块了。

这篇博文要实现的目标就是使用C语言实现一个简单的FFT模块，编译到固件中，最后在ESP32开发板上进行测试，并与在matlab中的计算结果对比。

一、FFT是什么？

FFT是快速傅里叶变换，不懂的童鞋请移步《信号与系统》《数字信号处理》

二、实现过程

1.编写C源码

代码如下

ez_fft.c

#include "ez_fft.h"

/**
 * @brief bit reverse
 * 
 * @param x the number to be reversed
 * @param bit_len .eg: for N in 1~7, bit len = 3. (as 2^3 = 8)
 * @return int result
 */
int bit_reverse(int x, int bit_len)
{
    int ret = 0;
    for(int i = 0; i < bit_len; i ++)
    {
        if(x & (0x1 << i))
            ret |= 0x1 << (bit_len - 1 - i);
    }
    return ret;
}

/**
 * @brief complex multiply
 * 
 * @param a 
 * @param b 
 * @return complex_t 
 */
complex_t multi_c(complex_t a, complex_t b)
{
    complex_t ret;
    ret.re = a.re * b.re - a.im * b.im;
    ret.im = a.re * b.im + a.im * b.re;
    return ret;
}

/**
 * @brief complex add
 * 
 * @param a 
 * @param b 
 * @return complex_t 
 */
complex_t add_c(complex_t a, complex_t b)
{
    complex_t ret;
    ret.re = a.re + b.re;
    ret.im = a.im + b.im;
    return ret;
}

/**
 * @brief complex subtraction
 * 
 * @param a 
 * @param b 
 * @return complex_t a - b
 */
complex_t sub_c(complex_t a, complex_t b)
{
    complex_t ret;
    ret.re = a.re - b.re;
    ret.im = a.im - b.im;
    return ret;
}

/**
 * @brief 旋转因子
 * 
 * @param N 下标
 * @param k 上标
 * @return complex_t 
 */
complex_t W(int N, int k)
{
    complex_t ret;
    ret.re =  cos(2 * PI / N * k);
    ret.im = -sin(2 * PI / N * k);
    return ret;
}

/**
 * @brief 蝶形运算
 * 
 * @param a 
 * @param b 
 * @param w 旋转因子 
 */
void butterfly_compu(complex_t* a, complex_t* b, complex_t w)
{
    complex_t tmp1, tmp2;
    tmp1 = add_c(*a, multi_c(*b, w));
    tmp2 = sub_c(*a, multi_c(*b, w));
    *a = tmp1;
    *b = tmp2;
}

/**
 * @brief 计算复数序列的摸，使用原缓冲区的前一半储存
 * 
 * @param buff 缓冲区，实部虚部交错
 * @param N 复数列长度，为缓冲区长度的一半
 */
void ezfft_abs(float* buff, int N)
{
    for (int i = 0; i < N; i++)
        *(buff + i) = pow(*(buff + i * 2) * *(buff + i * 2) + *(buff + i * 2 + 1) * *(buff + i * 2 + 1), 0.5);
}

/**
 * @brief 求N对2的对数
 * 
 * @param N 2的整次幂
 * @return int 当N不是2的整次幂时返回0
 */
int ezlog2(int N)
{
    int ret = 0;
    while (!(N&0x01))
    {
        N >>= 1;
        ret ++;
    }
    int bitsum = 0;
    for (int i = 0; i < sizeof(N) * 8; i++)
        if(N & (1 << i)) bitsum ++;
    if(bitsum != 1) ret = 0;
    return ret;
}

/**
 * @brief N点fft，注意N为2的整次幂
 * 
 * @param N fft点数
 * @param buff 缓冲区，实部虚部交错
 * @return int 正常返回0, 否则返回1
 */
int fft_N(int N, float* buff)
{
    if(ezlog2(N) == 0) return 1;
    complex_t* buff_c = (complex_t*)malloc(sizeof(complex_t) * N);
    //将原序列的序号比特翻转，并存入临时数组
    for (size_t i = 0; i < N; i++)
    {
        buff_c[i].re = *(buff + 2 * bit_reverse(i, ezlog2(N)));
        buff_c[i].im = *(buff + 2 * bit_reverse(i, ezlog2(N)) + 1);
    }
    complex_t *w = (complex_t *)malloc(sizeof(complex_t) * N / 2);
    for (int i = 0; i < N / 2; i++)
        *(w + i) = W(N, i);
    //蝶形运算
    for (int i = 0; i < ezlog2(N); i++)
    {
        for (int j = 0; j < N / (1<<(i + 1)); j++)
        {
            for (int k = 0; k < 1<<i; k++)
            {
                int x = (1<<(i + 1)) * j + k;
                int y = x + (1 << i);
                int z = (i != 0) * k * N / (1 << (i + 1));
                // printf("butterfly_compu(&buff_c[%d], &buff_c[%d], w%d);\n",x, y, z);
                butterfly_compu(&buff_c[x], &buff_c[y], w[z]);
            }
        }
    }

    for (size_t i = 0; i < N; i++)
    {
        *(buff + i * 2) = buff_c[i].re;
        *(buff + i * 2 + 1) = buff_c[i].im;
    }
    free(buff_c);
    free(w);
    return 0;
}

ez_fft.h

#ifndef __EZ_FFT__H
#define __EZ_FFT__H

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI (3.1415926535)

typedef struct complex
{
    float re;
    float im;
}complex_t;

void butterfly_compu(complex_t* a, complex_t* b, complex_t w);
int fft_N(int N, float* buff);
int bit_reverse(int x, int bit_len);
void ezfft_abs(float* buff, int N);
int ezlog2(int N);
complex_t W(int N, int k);
complex_t add_c(complex_t a, complex_t b);
complex_t sub_c(complex_t a, complex_t b);
complex_t multi_c(complex_t a, complex_t b);

#endif

mod_ezfft.c （此文件是通向MicroPython的大门）

#include "ez_fft.h"
#include "py/runtime.h"
#include "py/objarray.h"


STATIC mp_obj_t ez_memtest(mp_obj_t len_obj)
{
    int len = mp_obj_get_int(len_obj);
    char* p = (char*)malloc(len);
    if(p == NULL) return mp_obj_new_bool(0);
    for (size_t i = 0; i < len; i++)
    {
        *(p + i) = 0x33;
    }
    free(p);
    return mp_obj_new_bool(1);
}
STATIC MP_DEFINE_CONST_FUN_OBJ_1(ez_fft_memtest_obj, ez_memtest);

STATIC mp_obj_t ez_fftN(mp_obj_t buffer)
{
    mp_obj_array_t* array = (mp_obj_array_t *)buffer;
    if(array->typecode != 'f')
    {
        printf("please input float array");
        return mp_obj_new_bool(0);
    }
    return mp_obj_new_bool(!fft_N(array->len, (float*)(array->items)));
}

STATIC MP_DEFINE_CONST_FUN_OBJ_1(ez_fft_fftN_obj, ez_fftN);

STATIC mp_obj_t ez_fft_abs(mp_obj_t buffer)
{
    mp_obj_array_t* array = (mp_obj_array_t *)buffer;
    if(array->typecode != 'f')
    {
        printf("please input float array");
        return mp_obj_new_bool(0);
    }
    ezfft_abs(array->items, array->len / 2);
    array->items = (float *)m_realloc(array->items, sizeof(float) * array->len /2);
    array->len = array->len / 2;
    return buffer;
}

STATIC MP_DEFINE_CONST_FUN_OBJ_1(ez_fft_abs_obj, ez_fft_abs);

STATIC const mp_rom_map_elem_t ezfft_module_globals_table[] = {
    { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_ezfft) },
    { MP_ROM_QSTR(MP_QSTR_fft), MP_ROM_PTR(&ez_fft_fftN_obj) },
    { MP_ROM_QSTR(MP_QSTR_memtest), MP_ROM_PTR(&ez_fft_memtest_obj) },
    { MP_ROM_QSTR(MP_QSTR_ezabs), MP_ROM_PTR(&ez_fft_abs_obj) },
};
STATIC MP_DEFINE_CONST_DICT(ezfft_module_globals, ezfft_module_globals_table);

// Define module object.
const mp_obj_module_t ezfft_user_cmodule = {
    .base = { &mp_type_module },
    .globals = (mp_obj_dict_t *)&ezfft_module_globals,
};

// Register the module to make it available in Python.
MP_REGISTER_MODULE(MP_QSTR_ezfft, ezfft_user_cmodule, 1);

2.编译固件

好消息，MicroPython项目现在支持esp-idf了！！

也就是说我们现在可以通过编辑CMakeLists.txt来添加自己的源代码了，十分的方便

先把MicroPython项目clone到本地

git clone --recurse-submodules https://github.com/micropython/micropython.git

如果嫌太慢可以在这里免费下载：> micropython.tar.gz

下载完毕之后，先进入micropython目录

 make -C mpy-cross/

等待编译完成，然后创建一个目录用于存放我们自己的代码

cd ports/esp32
mkdir -p usermod/ezfft

将ez_fft.c ; ez_fft.h ; mod_ezfft.c这三个文件复制到usermod/ezfft目录下

然后修改 micropython/ports/esp32/main/CMakeLists.txt 文件，修改过的文件请自行下载，已经注释好了

CMakeLists

最后进入micropython/ports/esp32目录，执行

make deploy

这样就完成了，下面给出最后的固件，方便大家使用

点击这里下载👇

固件

用下面的命令烧录

esptool  -b 460800 --before default_reset --after hard_reset write_flash --flash_mode dio --flash_freq 40m --flash_size 4MB 0x8000 partition-table.bin 0x1000 bootloader.bin 0x10000 micropython.bin

3.测试

通过串口终端连接ESP32开发板，尝试能否导入模块：

import ezfft

没报错即为导入成功

接下来编写main.py脚本测试。以8点FFT为例：

main.py

import ezfft, array

# 8 点FFT测试
# 待测函数: 实函数 x[N] = a^N
def test():
    x = array.array("f", [])
    # 填充数据
    a = 0.8
    for i in range(8):
        # 实部虚部交错
        x.append(a ** i)
        x.append(0)
        
    # 注意ezfft.fft(x) 函数会返回一个bool类型，当x的长度不是2的整次幂是返回False
    # 否则返回True，结果会储存在x中，也就是说原数据会被破坏
    if ezfft.fft(x) :
        for i in range(8):
            print("{0} \t+\t i {1}".format(x[i*2], x[i*2+1]))

测试结果如下：

与matlab结果对比

感兴趣的朋友可以试试更多点数的FFT，在ESP32上运行1024点FFT还是没什么压力的: )

总结

MicroPython项目可移植性强，方便扩充，我们可以用C语言实现各种功能并作为模块编译进固件中，可以获得比纯Python代码更佳的性能

2021/3/5

2021/5/12更新，加入IFFT代码


/**
 * @brief N点ifft，注意N为2的整次幂
 * 将FFT的旋转因子的kn改为-kn，再整体乘以1/N就是IFFT了
 * @param N ifft点数
 * @param buff 缓冲区，实部虚部交错
 * @return int 正常返回0, 否则返回1
 */
int ifft_N(int N, float* buff)
{
    if(ezlog2(N) == 0 ) return 1;
    complex_t* buff_c = (complex_t*)malloc(sizeof(complex_t) * N);
    //将原序列的序号比特翻转，并存入临时数组
    for (size_t i = 0; i < N; i++)
    {
        buff_c[i].re = *(buff + 2 * bit_reverse(i, ezlog2(N)));
        buff_c[i].im = *(buff + 2 * bit_reverse(i, ezlog2(N)) + 1);
    }
    //计算旋转因子
    complex_t *w = (complex_t *)malloc(sizeof(complex_t) * N / 2);
    for (int i = 0; i < N / 2; i++)
        *(w + i) = W(-N, i);    //注意因为是IFFT,这里的参数是-N
    //蝶形运算
    for (int i = 0; i < ezlog2(N); i++)
    {
        for (int j = 0; j < N / (1<<(i + 1)); j++)
        {
            for (int k = 0; k < 1<<i; k++)
            {
                int x = (1<<(i + 1)) * j + k;
                int y = x + (1 << i);
                int z = (i != 0) * k * N / (1 << (i + 1));
                // printf("butterfly_compu(&buff_c[%d], &buff_c[%d], w%d);\n",x, y, z);
                butterfly_compu(&buff_c[x], &buff_c[y], w[z]);
            }
        }
    }
    for (size_t i = 0; i < N; i++)
    {
        *(buff + i * 2) = buff_c[i].re / (float)N;      //注意因为是IFFT,这里要除以N,下同
        *(buff + i * 2 + 1) = buff_c[i].im / (float)N;
    }
    free(buff_c);
    free(w);
    return 0;
}