|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
FFT是离散傅立叶变换的快速算法,可以将一个信号变换5 n1 t* u9 @# q! e, d6 _
到频域。有些信号在时域上是很难看出什么特征的,但是如
1 m; I; a, t. L7 U7 a果变换到频域之后,就很容易看出特征了。这就是很多信号
[* ]$ ?5 K- d, C" a c' S分析采用FFT变换的原因。另外,FFT可以将一个信号的频谱
8 u3 v0 d+ `) A+ m提取出来,这在频谱分析方面也是经常用的。
, f+ J4 r/ ]: k% _4 P 虽然很多人都知道FFT是什么,可以用来做什么,怎么去
; p9 _: G c- p9 T3 f: o& n) p做,但是却不知道FFT之后的结果是什意思、如何决定要使用
, f+ x1 B P: \8 Y/ h多少点来做FFT。$ a/ u4 U$ X. a( d: d1 G1 n4 T" x
6 U% ^9 M" q0 E" P* M) p" K; n* V 现在圈圈就根据实际经验来说说FFT结果的具体物理意义。6 `1 l: j( g2 L: \5 R
一个模拟信号,经过ADC采样之后,就变成了数字信号。采样6 h/ g* q7 e' h/ \$ ~
定理告诉我们,采样频率要大于信号频率的两倍,这些我就2 M+ O& d( s$ k; M
不在此罗嗦了。 F1 u# E) v" B6 B! u8 V
8 a9 }( P" b8 T7 N# ? 采样得到的数字信号,就可以做FFT变换了。N个采样点,
5 h. R+ |# T8 p/ j$ S经过FFT之后,就可以得到N个点的FFT结果。为了方便进行FFT4 h' _$ s- a: ~) }% q
运算,通常N取2的整数次方。
: j* a* Z o# r) U6 q: [1 ]" R9 V9 d% u. F/ _) d6 L$ t
假设采样频率为Fs,信号频率F,采样点数为N。那么FFT
& g3 Q/ ]( f0 E9 @$ r3 `% l之后结果就是一个为N点的复数。每一个点就对应着一个频率
* `7 e; X: A4 p( H点。这个点的模值,就是该频率值下的幅度特性。具体跟原始- A6 u' w" A- H9 r1 V- ?3 W1 ~
信号的幅度有什么关系呢?假设原始信号的峰值为A,那么FFT ^7 v, R! R1 t( S
的结果的每个点(除了第一个点直流分量之外)的模值就是A3 D6 V6 p! w3 |8 l, C2 V8 a. N) J
的N/2倍。而第一个点就是直流分量,它的模值就是直流分量& A& [8 K. _+ I
的N倍。而每个点的相位呢,就是在该频率下的信号的相位。" |" O( \ m7 F/ J5 }" x3 T
第一个点表示直流分量(即0Hz),而最后一个点N的再下一个2 ~3 E; q" `( ?. r
点(实际上这个点是不存在的,这里是假设的第N+1个点,也0 l' Y1 N5 c" F7 y( `1 k- S; s6 |/ P
可以看做是将第一个点分做两半分,另一半移到最后)则表示% A( H0 `7 T# {: v
采样频率Fs,这中间被N-1个点平均分成N等份,每个点的频率$ W& j' [% k' I$ p1 e0 x0 `
依次增加。例如某点n所表示的频率为:Fn=(n-1)*Fs/N。: A3 J9 v. f3 k$ Z; u% q
由上面的公式可以看出,Fn所能分辨到频率为为Fs/N,如果/ l- ]8 ~$ y; I2 B0 X- U
采样频率Fs为1024Hz,采样点数为1024点,则可以分辨到1Hz。
- I# a: y: I$ f+ [" {$ P1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒
! ~7 k" g! Q4 r8 C% A时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时
$ W I/ \2 `7 p: O+ V间的信号并做FFT,则结果可以分析到0.5Hz。如果要提高频率
2 G% }6 `4 ]( o分辨力,则必须增加采样点数,也即采样时间。频率分辨率和* }8 m5 E5 f( `: L
采样时间是倒数关系。
7 S! W1 _0 y) _% L 假设FFT之后某点n用复数a+bi表示,那么这个复数的模就是
. P$ W: j7 `) K/ q2 NAn=根号a*a+b*b,相位就是Pn=atan2(b,a)。根据以上的结果,
3 u. ~8 D7 k1 q! s9 U就可以计算出n点(n≠1,且n<=N/2)对应的信号的表达式为:( h( D$ ~. a. h. i" O* j
An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。- h8 m' i$ q+ z0 a0 {
对于n=1点的信号,是直流分量,幅度即为A1/N。7 l* X& {, {& X' a
由于FFT结果的对称性,通常我们只使用前半部分的结果,
( u0 l% q) H: [ m' C, P即小于采样频率一半的结果。
: Y% D- O5 B) r' }9 _: n+ M/ L( F' n$ ~6 N
好了,说了半天,看着公式也晕,下面圈圈以一个实际的
[* ?6 ?# [! i信号来做说明。3 Y# N, C3 l% q) M9 U! l7 m
7 q- {3 `0 Q* W, P x! h8 t
假设我们有一个信号,它含有2V的直流分量,频率为50Hz、* M3 x$ n" T B. H1 z
相位为-30度、幅度为3V的交流信号,以及一个频率为75Hz、
. ?. w* J V# o! O- m4 C8 m2 ~相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下:; c) A1 _$ q- B$ V& {% X! X+ y
0 w; K# c- d5 D3 @ cS=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)5 D0 b8 P5 V( B& y5 b0 K- d2 a
' Y8 ?; Z; m" m% k) ? 式中cos参数为弧度,所以-30度和90度要分别换算成弧度。
- ~9 H$ a9 P* S, u+ V& C+ U5 F b$ e; G我们以256Hz的采样率对这个信号进行采样,总共采样256点。
( p9 |4 l$ @, X* x; }. I0 b按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个# h6 T- v* l9 X2 O6 v
点之间的间距就是1Hz,第n个点的频率就是n-1。我们的信号% M% M- C/ o/ i# o
有3个频率:0Hz、50Hz、75Hz,应该分别在第1个点、第51个点、7 W3 K* g4 M8 x( \; n4 z* L2 u. g
第76个点上出现峰值,其它各点应该接近0。实际情况如何呢?6 O/ Q! x; p- m. y D
我们来看看FFT的结果的模值如图所示。
6 |% r$ x( K% u" v5 q
W' m: T3 [0 ~0 E, P8 }8 W5 h- W! P8 I( r2 w3 d1 h
2 n9 N' j( @7 A" v; S& h
图1 FFT结果& `2 b: ]; J, N" B3 S' ~2 M
从图中我们可以看到,在第1点、第51点、和第76点附近有+ T9 O% `) ^6 U+ u+ T I+ l& ^$ J
比较大的值。我们分别将这三个点附近的数据拿上来细看:
% S4 c: Q- F" j4 r$ n0 w1点: 512+0i
8 p, Y& C6 x5 F2点: -2.6195E-14 - 1.4162E-13i, G9 h6 p8 C- E# d) G
3点: -2.8586E-14 - 1.1898E-13i
( L1 b0 D; A: q' M
9 {( g0 A7 \8 n. q* S9 C50点:-6.2076E-13 - 2.1713E-12i
, g( D7 ~ b/ g8 ?" I8 w: c51点:332.55 - 192i
( k- D$ ~4 _2 Q52点:-1.6707E-12 - 1.5241E-12i7 e* n/ n) g3 T* u- Y
- Q) t5 f& Z' ?9 \1 \4 w$ o75点:-2.2199E-13 -1.0076E-12i6 Q& j/ w5 o/ f- E0 h
76点:3.4315E-12 + 192i
% n( }# }. z# C9 H2 R% ? ~9 l77点:-3.0263E-14 +7.5609E-13i4 w; t: B0 g' p* A8 {) E
, P9 N# L6 V$ @! ?2 V
很明显,1点、51点、76点的值都比较大,它附近的点值
3 d2 D+ r3 @/ b8 x都很小,可以认为是0,即在那些频率点上的信号幅度为0。
1 s$ K6 U, O% C) A# N接着,我们来计算各点的幅度值。分别计算这三个点的模值,
8 ]! e0 E: f% a7 j. @结果如下:7 c* k/ i' n( M+ u
1点: 512
8 v* ]4 ]# f. _' p51点:384' n- l5 M X o0 ]. s; s8 k
76点:1927 o) ~+ V: {/ q& n5 e' K
按照公式,可以计算出直流分量为:512/N=512/256=2;8 K7 G' R. s) T: V; v d; ^* }
50Hz信号的幅度为:384/(N/2)=384/(256/2)=3;75Hz信号的
& l5 P& S \& {$ P2 A# E4 w幅度为192/(N/2)=192/(256/2)=1.5。可见,从频谱分析出来2 ~0 b: p9 L& x
的幅度是正确的。! ]2 E3 r, T, I8 \4 `+ ~
然后再来计算相位信息。直流信号没有相位可言,不用管
- B& D3 l) ~0 o+ x" }5 x它。先计算50Hz信号的相位,atan2(-192, 332.55)=-0.5236,
]8 {2 r! K% g- l3 o! o$ b2 g结果是弧度,换算为角度就是180*(-0.5236)/pi=-30.0001。再% G0 |+ J- `9 ^* w/ D- r
计算75Hz信号的相位,atan2(192, 3.4315E-12)=1.5708弧度,/ w( s3 v* \) z; C7 v: G7 n. r
换算成角度就是180*1.5708/pi=90.0002。可见,相位也是对的。( }- p0 i0 o% J) p$ i
根据FFT结果以及上面的分析计算,我们就可以写出信号的表达# x3 U3 G8 q) L4 X
式了,它就是我们开始提供的信号。 k3 C" ~! [1 ~$ b/ s
: ^, |; Y: ~0 y0 k+ D
总结:假设采样频率为Fs,采样点数为N,做FFT之后,某% J3 Q) J! R) x5 q
一点n(n从1开始)表示的频率为:Fn=(n-1)*Fs/N;该点的模值
4 H0 O( C& P% h, W除以N/2就是对应该频率下的信号的幅度(对于直流信号是除以
, E: Y/ y8 w/ g/ Z+ i* tN);该点的相位即是对应该频率下的信号的相位。相位的计算
) e; a# y5 i1 S# ~. j可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角/ R$ f2 o( G% e A' o2 b8 }! v
度值,范围从-pi到pi。要精确到xHz,则需要采样长度为1/x秒
# i6 a9 K6 z5 E2 ? w! T的信号,并做FFT。要提高频率分辨率,就需要增加采样点数,$ V4 y" n7 w+ i5 C2 E; d4 t! W) K
这在一些实际的应用中是不现实的,需要在较短的时间内完成
: K( k5 d4 z" e. p( { J- M分析。解决这个问题的方法有频率细分法,比较简单的方法是$ w$ p- |. }8 A7 G r0 x
采样比较短时间的信号,然后在后面补充一定数量的0,使其长度
. j) L5 ]! o1 v达到需要的点数,再做FFT,这在一定程度上能够提高频率分辨力。3 y; R; u: w" ], b
具体的频率细分法可参考相关文献。% Q8 Z5 z3 r) }- N" t: \3 L
5 m* y& n0 x( s' _+ B( r
[附录:本测试数据使用的matlab程序]
! L3 z: |' |1 P! H9 ]close all; %先关闭所有图片0 e; i3 U4 d7 U9 s
Adc=2; %直流分量幅度
/ p+ N' H5 R! S: O$ fA1=3; %频率F1信号的幅度
: h- G! y2 n% F" |: }3 DA2=1.5; %频率F2信号的幅度- d7 Y. F5 _0 H: A. A/ [* C- G
F1=50; %信号1频率(Hz)! R* W; J8 C u6 s7 b
F2=75; %信号2频率(Hz)
; a# P# z9 }2 R0 r, q3 m- H4 Z2 r% ^Fs=256; %采样频率(Hz)! T3 g a8 E$ f+ L7 i. H
P1=-30; %信号1相位(度)- _8 F8 ]0 o! ]5 ~6 I
P2=90; %信号相位(度)
6 W; ^1 e+ q9 j( W- lN=256; %采样点数
: M; b( W8 y! o( d; U9 z* ^t=[0:1/Fs:N/Fs]; %采样时刻
: h# w8 |7 X% z5 |/ \0 f% G, `2 R
3 g( ~6 Y9 B- N% \6 P3 T%信号; \: A) a2 p& J' Q$ ~
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
) ^1 D: s" {: M%显示原始信号
! E; V4 |4 H" |" F2 wplot(S);
: _" C. x1 b, |1 z6 \title('原始信号');
" e4 K2 J0 a$ N
& l# l, i( v. ?7 {! V& E, _& v) d6 [figure;
( v! R/ @9 `& ` A( }/ ?7 @Y = fft(S,N); %做FFT变换* g! c$ T" Y4 R" p$ m% K5 D+ z# E( O0 b
Ayy = (abs(Y)); %取模0 p' y, y! `7 c6 o1 F& c
plot(Ayy(1:N)); %显示原始的FFT模值结果$ E$ n) Y4 V! l& I
title('FFT 模值');' Y( v& R1 d# o4 d7 `
: w1 {3 W& t7 x5 e6 E% M" j
figure;
$ L: z/ J1 [! Z- e, Y- EAyy=Ayy/(N/2); %换算成实际的幅度3 H9 b% K/ Y, U- W& z' s
Ayy(1)=Ayy(1)/2;
% f- n! Z3 C, n; H1 p) R3 K) JF=([1:N]-1)*Fs/N; %换算成实际的频率值- ?% S9 h6 Z! b" r% ^
plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果
0 @& |. E/ u9 ]/ U- B: z/ s8 ititle('幅度-频率曲线图');9 v9 w1 Q8 H X" m2 I. M! ~4 u
0 L/ ]7 V' N- Z* d0 Ffigure;5 h" s! {1 a1 r" b6 V) Y7 T
Pyy=[1:N/2];: \$ B* J4 Y% V1 f! V
for i="1:N/2"
+ }2 N3 w5 N( y* ]* s" |* i' |Pyy(i)=phase(Y(i)); %计算相位8 ?: y5 I6 s f& J S( D1 D2 T
Pyy(i)=Pyy(i)*180/pi; %换算为角度
! _8 ~" O4 N# E. Z$ J' \* wend;! b% u/ D$ w x
plot(F(1:N/2),Pyy(1:N/2)); %显示相位图- k9 Q. Y9 g- Z$ ~$ D
title('相位-频率曲线图');. R N5 [4 j, o+ O$ f- a \
- v3 k. b# j9 O6 J; t) L" L看完这个你就明白谐波分析了。
6 Y& ]6 A" [8 e q/ ]( Q |
|