快速傅里叶变换(FFT)是一种数字信号处理中常用的技术,用于将快速序列转换为频域表示。在嵌入式系统中,如基于STM32的微控制器,实现FFT可以帮助解决信号处理的需求,例如声音处理、图像处理等。本文将介绍基于STM32的离散傅里叶变换的原理、实现方法和应用。
8 |7 H' q$ @& Z) z& l3 P- b( @
1 ]7 ?, g. c2 c. U
7 ~( {, }! f& c3 t+ d4 `
* ^* Y9 v: D, K0 \6 c1 W- C
FFT是一种将时域序列转换为频域表示的技术,它将一个序列的N个采样点映射到频域中N个频率分量。其数学表达式如下:6 O+ ]0 O1 @) U
# c ?# W4 g- G0 I5 F
# X: E3 {0 j8 |) F# ~% A% S7 a# _5 p) ?" ^: M# u
其中,x(n) 是输入序列,X(k) 是输出的频域表示。' o9 ]8 J8 E g6 z& S
3 |' ^( e/ A; [4 m/ U准备工作:
) U) k/ x2 y* s w) Q; ?& m( D K! L e8 {0 y+ |1 l
+ d8 T7 j4 w8 g" h8 c
$ U1 r. b4 s0 r5 n& G8 _; [$ |* I& m" F. r: J& Z' _. {
- G9 m: g! K: N' l2 f9 @
' c/ p! k" X C* W$ PKeil中的DSP库(Digital Signal Processing Library,数字信号处理库)是针对ARM Cortex-M处理器系列的一组软件库,用于提供各种数字信号处理功能的支持。这些库提供了一系列优化过的算法,可以帮助开发人员在嵌入式系统中高效地实现音频处理、图像处理、通信系统等各种信号处理应用。' }5 v9 h# p; Y, v$ m; ^
9 R! q/ E" T' P$ E! P( n
因此我们需要在Keil中安装我们的DSP库。! {. ?1 m. s# Q# |! q
- #include "arm_math.h" // 包含DSP库
复制代码 ! f$ h6 s0 m+ I( O( L
首先包含我们的DSP库。
$ d7 n7 t" u& t; \- h6 V4 N5 [7 e2 y* y
! G+ [0 |7 u0 }* M, y- ^; a) h5 f2 F- M. a$ l! f* U2 ]
- #define FFT_LENGTH 100
( A3 x, H! e" S - // 输入序列- Y" `( Y0 [+ Y7 P) ]& X7 Y+ t
- float32_t inputSignal[FFT_LENGTH*2];
+ t' G. I0 n" _1 M, r - 6 f E" Z. a) Q* Q
- // 输出序列,存储变换后的结果9 N. a- l% t- A
- float32_t outputSignal[FFT_LENGTH];
复制代码
4 r0 i2 E; T' F5 y' _, t1 m6 }8 B定义FFT的的输入和输出数组还有数组长度
; v! o, [, q9 t8 o1 u& t( g- arm_status status;/ l" A' P4 \$ D. g' L5 b
- arm_cfft_radix4_instance_f32 fft_inst;
! J; f' h0 P' ?2 W* K - status = arm_cfft_radix4_init_f32(&fft_inst, FFT_LENGTH,0,1);
复制代码- void arm_cfft_radix4_init_f32(. I, i y ?; c) L' W( v# m
- arm_cfft_radix4_instance_f32 * S,
! N: J8 f# u( i; @0 w - uint16_t fftLen," W. t; r. ^' L# f: y: N3 v; x) O. t
- uint8_t ifftFlag,. x" K% l& z5 F' U5 [5 B( Y
- uint8_t bitReverseFlag; w, L# P8 E& Q# Z0 J" j9 u
- );
复制代码
- P0 U% E0 }$ V3 p+ Z定义一个状态变量用来显示FFT的初始化是否成功。
' k5 T+ H7 I. L3 u j9 P定义一个FFT的配置变量。) Q5 m% o" K Y( Z
初始化FFT。8 X6 F" h7 f; Q) S A+ \
S:指向 arm_cfft_radix4_instance_f32 结构体的指针,该结构体定义了 FFT 实例的状态信息。
& N+ O6 ]( t9 u* Q, w5 x; }! Q9 ]. Z6 N
fftLen:FFT 的长度。
' v2 n4 c- o0 c6 {
$ y% B( \+ H8 F& yifftFlag:指定是否进行逆变换。如果为 1,则表示初始化的是逆变换的 FFT;如果为 0,则表示初始化的是正变换的 FFT。% r& ^- W. u" J- v: K
& v1 m, ^8 M( t0 L+ B% F
bitReverseFlag:指定是否进行比特翻转。如果为 1,则表示进行比特翻转;如果为 0,则表示不进行比特翻转。
+ ^* l8 [- m( C. M; `8 u8 g5 C8 _+ n6 y2 [" u
在FFT算法中,比特(bit)反转是一种关键的步骤,用于将输入数据重新排列为正确的顺序,以便在后续的计算中进行有效处理。7 O3 n4 t0 I A5 L3 Y& v
4 N8 S9 m- k* P1 K$ P. s( D. T
当进行快速傅立叶变换时,算法要求输入数据的顺序是按照特定的方式排列的。特别是在使用基于分治法的算法(如Cooley-Tukey算法)时,输入数据的顺序必须满足按照一定规律的排列。
7 \2 u3 [" x- t! k# g0 C7 h! ]: K, n/ Z
在实际的FFT实现中,最常见的方式是通过比特反转来重新排列输入数据。比特反转就是将输入数据的比特位(二进制位)的顺序进行颠倒。这是因为在FFT算法中,数据会被分组,并按照一定规则进行反转,以便在每个阶段的运算中,数据可以正确地与其它组合进行配对。! h# t" V- o0 O) A' C1 u# \5 l
2 ~% J' l$ h4 x举个简单的例子,假设有一个长度为8的数据序列,按照0到7的顺序排列: m) X6 l" `, [: [
, `/ h# v/ C. Z1 x
0 1 2 3 4 5 6 7
6 v9 R* b, Z& y- M& q
7 U$ I; J; e; v3 o4 h* G* b. v在进行FFT时,需要按照一定规则重新排列这些数据。比特反转操作将会对这个数据序列进行如下的重新排列:4 R9 j: g" g9 I0 q8 U) C5 B
% B3 G+ }. ]5 h2 k0 4 2 6 1 5 3 77 l* c2 g/ a" Y9 K8 e" K" @3 S% c
" S( H2 l/ i) r7 G9 u在FFT算法的每个阶段中,这种重新排列都会使得数据正确地与其它组合进行配对,从而实现快速傅立叶变换的计算。
' x+ D1 U, m- @2 O: y
/ T4 }% e0 L; [6 f2 |& ]. o' w进行FFT并转换为模值
1 ~; X: N) p' o' D3 V- arm_cfft_radix4_f32(&fft_inst,inputSignal); //FFT计算
* p' }8 `9 p) c: C' l - arm_cmplx_mag_f32(inputSignal,outputSignal,FFT_LENGTH); //取模得幅值
复制代码
# V! X+ n" F. S4 k2 s# y对输入数组进行FFT变换,并将FFT的结果转化为模值。
" |- J0 k/ {. W0 }/ M, @
9 _0 T% \7 O* Q+ u- U* f测试
2 F: S$ Q% q. ?1 q/ l我们进行一个简单的测试
5 s7 X! _8 T9 k) y4 H4 }# E- #define FFT_SIZE 1024
6 R0 t' i# z* q; W; n - #define SAMPLE_RATE 1000
H, A) }! d) N$ }9 `2 ?# N( D - #define NUM_SAMPLES 10001 J/ d) J" B" }6 o: P8 q" S
- #define FREQ_OF_INTEREST 100
, m: A) j( I- G; Z4 ~3 E! Z+ b - for (int i = 0; i < NUM_SAMPLES; i++) {/ d j" j, G9 _, g' E1 n
- float32_t t = (float32_t)i / SAMPLE_RATE;
9 i- ^; d$ f. b \$ y6 \0 ~! C% ? - float32_t sin_value = sinf(2 * PI * FREQ_OF_INTEREST * t); // 计算正弦波值
+ H" {0 l8 a0 h/ h - inputSignal[i * 2] = sin_value; // 实部( t" O/ _5 v6 S% v
- inputSignal[i * 2 + 1] = 0; // 虚部- X0 t; i, f$ Y5 f* N. X6 u
- }
复制代码
9 J- G8 l& U9 x( ^; b8 n+ ~一千个点的采样值,频率假设为100HZ作为输入信号。' q, R4 ~! B1 }$ c* v
- for (int i = 0; i < FFT_SIZE; i++) {
! t, q0 r4 }. Q b - // 计算复数的模值1 s; T% c8 N! j1 G; Q; }' u
- float32_t real = inputSignal[2 * i];
: b8 H- `5 L2 J, I: o% b - float32_t imag = inputSignal[2 * i + 1];
8 t# i) ]- ^6 b- E8 [+ O% ` - float32_t magnitude = sqrtf(real * real + imag * imag);
+ P% f |* ~* l8 ?2 ^8 \( r1 U/ ] -
0 @; n+ d" m4 g- v5 l - // 打印每个频率分量的模值
! v5 i6 \$ M; d/ A, L8 z - printf("Magnitude: %f\n", magnitude);
* E$ ?7 L; a9 s2 D - }
复制代码 ) J9 D0 n6 w1 o/ Z+ i
进行傅里叶变换后打印模值。% g- b9 y! h; D' s# m
/ O& C. ]! J$ F6 k
3 X8 \" Y% l4 D: S1 E0 M
4 R" @! T$ G9 H' K7 K: W
可以看到傅里叶变换执行成功。
/ \( [7 Z) P) y- for (int i = 0; i < NUM_SAMPLES; i++) {1 u" `. L5 C. d+ T
- float32_t t = (float32_t)i / SAMPLE_RATE;
# T8 A6 ]% \% w. _9 \- R6 n. ?! p - float32_t sin_value = sinf(2 * PI * FREQ_OF_INTEREST * t)+sinf(2 * PI * FREQ_OF_INTEREST * t*2)+sinf(3*2 * PI * FREQ_OF_INTEREST * t); // 计算正弦波值" k/ p5 u9 j9 B+ O
- inputSignal[i * 2] = sin_value; // 实部& a" P3 z: t% J. N4 t( ]( `
- inputSignal[i * 2 + 1] = 0; // 虚部
- f& x0 D- a' m, z9 q - }
复制代码
% C* i1 o3 [& O8 R我们将信号制作成100HZ+200HZ+300HZ的信号。
/ N: L) S+ ^: E( `8 ~6 ~5 k- f2 l- [( j' H1 k, l# t, a- U
9 M( u% W% \/ y3 M6 w
' w5 Y8 _# J' Y$ T& T# m3 m% C
) g* U+ {! q1 F8 V, C
' y0 W/ C \9 Q# h- i3 c2 V0 E转载自:电路小白$ x# l9 V! ~* R
如有侵权请联系删除+ X( M4 s# Y, f7 \! B" \7 r9 S
B% @9 w( y0 m& p6 b1 R
|
这个FFT不错,学习参考一下