文章目录
文章配套资源
背景
使用 STM32 进行 FFT 时，由于官方函数的 FFT 长度仅为 4096，使 FFT 后的分辨率达不到要求，同时，MATLAB 可以被用来生成算法，于是通过 MATLAB 进行 FFT，并转为 C 语言的想法便由此而生，同时记录一下学习经历。
介绍
从一个初学者的角度来讲，我们一定知道频谱图，比如：
图1 的两个注意点
x轴
x 轴代表的是频率，该图的 x有效长度为 0->4096，图中标记点的 x轴是1025，那么这两个参数的物理意义是什么呢？
了解频谱图的人都知道，我们可以很直观的从频谱图中得知采集的信号的频率，那么这个频率是怎么算的呢？
首先引入三个概念，采样频率和信号频率与傅里叶点数，采样频率指的是每秒从连续信号中提取并组成离散信号的采样个数，信号频率指的是信号源产生的信号，通俗来讲，大家都玩过水果忍者这个游戏，信号就是游戏中的西瓜，采样就是游戏中的刀，采样频率就是挥刀的速度，挥刀速度越快，西瓜切的份数就越多。这样，便从连续到了离散，类似ADC。采样频率越快，信号频率越慢，离散的效果也就更好。傅里叶点数是指采样频率能够分成多少份，当采样频率除以傅里叶点数后，得到的就是频谱图中 x 轴的频率分辨率。
上面我们初步了解了采样频率和信号频率与傅里叶点数，现在来回答上面的两个问题。
x的有效长度指的是傅里叶点数，也就是采样频率的均摊点的个数（个人语言），图一的标记点的 x轴数据为1025 ，x轴的有效长度为4096 ，通过程序可知，采样频率为200e3 ，所以
频谱图的频率分辨率 = 采样频率 / 傅里叶点数 = 200e3 / 4096 = 48.828
标记点的频率 = x数据 * 分辨率 = 1025 * 48.828 = 50048
由于频率与分辨率不是整数倍，所以会与理论由偏差。
由程序可得，信号频率为：50e3。
以上便是频率的基础知识。
% 图一的程序
clear;
clc;
clear;
clc;
f = 50e3;
fs = 200e3;
fd = 1 / fs;
flen = 4096;
t = 0 : fd : fd * (flen-1);
tdate = sin(2*pi*f*t);
fft_date = zeros(1, 2*flen);
abs_date = zeros(1, flen);
fft_date(1 : 2 : 2*flen - 1) = real(fft(tdate));
fft_date(2 : 2 : 2*flen) = imag(fft(tdate));
for i = 1 : length(abs_date)
abs_date(i) = sqrt(power(fft_date(2*i - 1), 2) + power(fft_date(2*i), 2));
end
% 图二的程序
clear;
clc;
f = 50e3;
flen = 4096;
fs1 = 200e3;
fd1 = 1 / fs1;
t1 = 0 : fd1 : fd1*flen;
fdate1 = sin(2*pi*f*t1);
fs2 = 500e3;
fd2 = 1 / fs2;
t2 = 0 : fd2 : fd2*flen;
fdate2 = sin(2*pi*f*t2);
subplot(2, 1, 1);
scatter(t1, fdate1);
title(['采样频率：', num2str(fs1)], ['信号频率：', num2str(f)]);
xlim([0, fd1*10]);
subplot(2, 1, 2);
scatter(t2, fdate2);
title(['采样频率：', num2str(fs2)], ['信号频率：', num2str(f)]);
xlim([0, fd2*10]);
% result = fft_my_fun(fdate1, flen);
% plot(result);
y轴
y 轴的数据代表的是幅值乘傅里叶点数，当x = 0时，y轴的数据表示直流分量的幅值乘傅里叶点数，当x > 0时，y轴的数据的二倍表示该频率分量的幅值乘傅里叶点数。
这是 FFT 的特性，按照个人理解而言，将频谱图的 x轴，首尾相连，形成一个环状，x = 0的位置表示的是直流分量，左边是频率的负分量，右边是频率的正分量，所以一般只看其中的一个尖峰即可。这也就是为什么当信号频率为采样频率的一半时，频谱图什么都没有的原因。
对于图一而已，标记点的y轴数据为2048，所以
该频率幅值为 = y轴数据 / 傅里叶点数 * 2 = 2048 / 4096 * 2 = 1；
由程序可得，信号幅值为：1。
以上便是幅值的基础知识。
MATLAB
创建文件
首先创造两个 .m 文件，第一个叫 fun.m，用来编写 FFT 的主要功能，类似于 C 语言的功能函数，第二个叫 test.m，用来测试编写的功能函数；（建议文件路径不要有空格！）
% fun.m
function abs_date = fun(tdate, fft_number)%#codegen
fft_date = zeros(1, 2 * fft_number); % 初始化FFT的输入数组
abs_date = zeros(1, fft_number); % 初始化abs的数组
fft_date(1:2:2 * uint16(fft_number) - 1) = real(fft(tdate, fft_number)); % 将实部放在奇数索引中
fft_date(2:2:2 * uint16(fft_number)) = imag(fft(tdate, fft_number)); % 将虚部放在偶数索引中
for i = 1 : fft_number
% abs_date(i) = abs(fft_date(2*i - 1), fft_date(2*i));
abs_date(i) = sqrt(power(fft_date(2*i - 1), 2) + power(fft_date(2*i), 2)); % 计算实部与虚部的模长
end
end
% test.m
clear;
clc;
fs = 200e3; % 采样频率
f = 50e3; % 信号频率
fd = 1 / fs; % 步长
flen = 4096; % 傅里叶长度
t = 0 : fd : fd * (flen - 1); % 时间
tdate = sin(2*pi*f*t) + 1; % 时域数据
abs = zeros(1, flen); % 最后的频谱数据
abs = fun(tdate, flen); % 通过功能函数得到频谱数据
plot(abs);
测试两个文件
注意，由于test.m文件调用了fun.m，所以两个文件应该放在同一个文件夹中，比如：
运行程序，得到 图一，运行该程序时，应当遵循 奈奎斯特定理，如下是一组测试数据，当频率小于采样频率的一半时，随着频率的增加，频谱图的尖峰会朝着中间走，当频率等于采样频率的一半，频谱图没有尖峰，当频率大于采样频率的一半，尖峰会越过中间轴继续前进，如图5.
采样频率 FFT点数 理论分辨值 信号频率 FFT后数据 100 256 0.3906 5 14 100 256 0.3906 10 27 100 256 0.3906 20 52 100 256 0.3906 30 78 100 256 0.3906 40 103 100 256 0.3906 45 116 100 256 0.3906 48 124 100 256 0.3906 49 126 100 256 0.3906 49.5 130 100 256 0.3906 49.6 130 100 256 0.3906 49.9 129 100 256 0.3906 50 无 100 256 0.3906 >50 从头开始
开始 MATLAB转C
MATLAB Coder
依次点击 APP -> MATLAB Coder （如图6），点击后，进入图7 界面，点击红色部分的...，选择fun.m，并且选择单精度（如图8）。点击Next。
MATLAB Coder 选择入口函数的数据类型，并进行检查
点击Let me enter input or global types directly，根据需要确定入口函数的数据类型（如图9），点击Next，选择测试文件test.m，该目的是测试功能函数是否满足要求（如图10），点击Check for Issues，测试成功后，进入图11.
MATLAB Coder 进行编译
点击Next后，进入图12，建议Toolchain选项卡直接选择自己的编译器，注意：建议支持 MATLAB 正版，这里可能需要安装 MATLAB 对对应编译器的支持包，盗版应该没有安装入口。这里我选择MinGW64 编译器。如图配置后，点击Generate。如图13、图14为成功编译界面。
恭喜大家，总工程已经完成四分之三。
Visual Studio
建立工程
建议大家，打开图14中 C Code对应的文件夹，将其中所有的.c .h文件复制到新的文件夹，比如 "C:\Users\14104\Desktop\temp"文件夹中，使用完毕后，直接删掉。别忘了还有\example中的.c .h文件。
打开Visual Studio，新建一个项目，并将刚刚复制的.c .h文件添加到工程里。如图15、图16。
点击编译，发现出现报错，如图17。经过搜索网上资料，可知：tmwtypes.h文件是 MATLAB 转为 C 语言后，不会自动添加的文件，通过查资料可以找到该文件的地址为MATLABL安装路径\extern\include\中（图18）。
随后，按照上述步骤，将tmwtypes.h文件添加到该项目后，点击编译。出现图19，表示该工程没有问题。
观察fun.c文件
找到与 MATLAB 文件名相同的fun.c文件，并找到入口函数。如图20。
void fun(const emxArray_real32_T *tdate, float fft_number, emxArray_real32_T *abs_date)
现在，我们已经找到了 MATLAB 生成的入口函数了，我们的目的是调用这个入口函数，我们观察入口函数及其参数的数据类型。
void fun(const emxArray_real32_T *tdate, float fft_number, emxArray_real32_T *abs_date)
// 该结构体与 C 语言中的 calloc 十分相似
struct emxArray_real32_T
{
<type> *data; // 数据数组，要求填入地址。
int *size; // 数据数组的大小，要求填入地址。
int allocatedSize; // 实际需要分配的内存大小，要求填入数值。
int numDimensions; // 数据数组的维度
boolean_T canFreeData; // 是否需要清除该数组
};
修改main.c文件
打开main.c文件，将main函数进行修改。代码已全部写上注释，大家应该可以看懂。注意：请注意程序中的 头文件 与 数据类型。
#include <stdio.h>
#include <stdlib.h>
#include <math.h> // 一定要加，否则 tdate.data的数据全部为整数
int main()
{
uint16_T num_t = 4096; // 定义时间长度
uint16_T flen = 4096; // 定义FFT的长度
uint32_T fs = 200e3, f = 50e3; // 定义采样率和信号频率
uint16_T size[2] = { num_t, flen }; // 定义长度
float fmin = fs / flen; // 定义最小分辨率
float f_t_point = f / fmin; // 定义频率点
float fd = (float)(1.0 / fs); // 定义步长
float* t = (float*)calloc(num_t, sizeof(float)); // 定义时间
emxArray_real32_T tdate, abs_date; // 定义时域和频域的变量，用于MATLAB测试
tdate.data = (float*)calloc(num_t, sizeof(float)); // 放时域数据
tdate.size = &size[0]; // 实际存放的元素数量
tdate.allocatedSize = size[0]; // 分配的内存大小，与上面不一样，这个是需要的个数
tdate.numDimensions = 1; // 维度
tdate.canFreeData = 1; // 可以清除
abs_date.data = (float*)calloc(flen, sizeof(float)); // 放时域数据
abs_date.size = &size[1]; // 实际存放的元素数量
abs_date.allocatedSize = size[1]; // 分配的内存大小，与上面不一样，这个是需要的个数
abs_date.numDimensions = 1; // 维度
abs_date.canFreeData = 1; // 可以清除
FILE* file = fopen("yang.csv", "w+"); // 文件初始化
printf("全部已经初始化完毕\n");
if (t == NULL || tdate.data == NULL || abs_date.data == NULL)
{
printf("内存分配出错\n");
}
if (t != NULL && tdate.data != NULL && abs_date.data != NULL)
{
printf("内存分配成功\n");
}
for (uint16_T i = 0; i < num_t; i++)
{
t[i] = i * fd; // 初始化时间函数
fprintf(file, "%f", t[i]); // 写入时间数据
fprintf(file, ","); // 换列
}
fprintf(file, "\n"); // 换行
printf("时间已填充完毕\n");
for (uint16_T i = 0; i < num_t; i++)
{
tdate.data[i] = sin(2 * 3.1415926 * f * t[i]); // 填写时域数据
fprintf(file, "%f", tdate.data[i]); // 写入时域数据
fprintf(file, ","); // 换列
}
fprintf(file, "\n"); // 换行
printf("时域数据已填充完毕\n");
//fclose(file); // 关闭对文件的权限使用
//printf("正在打开 yang.csv\n");
//system("yang.csv");
printf("MATLAB 算法开始进行\n");
fun(&tdate, (float)flen, &abs_date);
printf("MATLAB 算法结束进行\n");
for (uint16_T i = 0; i < flen; i++)
{
fprintf(file, "%f", abs_date.data[i]); // 写入频域数据
fprintf(file, ","); // 换列
}
fprintf(file, "\n"); // 换行
printf("频域数据已填充完毕\n");
fclose(file); // 关闭对文件的权限使用
printf("理论步长为 %f\n", fd);
printf("理论上分辨率为 %d / %d = %f\n", fs, flen, fmin);
printf("理论上频率点为 %d / %f = %f\n", f, fmin, f_t_point);
printf("正在打开 yang.csv\n");
system("yang.csv");
return 0;
}
查看效果
如图21，C文件几乎与 MATLAB 一模一样。至此，艺术已成。
结束语
至此，该文章已全部结束，感谢大家观看，以上理论基于个人理解，如有错误，欢迎指正。
927
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
