MATLAB之Fourier Transforms ( fft )

mytomorrow

: b. ]* c) j! L# F) _/ w- g7 U( z
Fourier Transforms


. P! \$ W  E6 y傅立叶变换是将在时间或空间中采样的信号与频率采样的相同信号相关联的数学公式。 在信号处理中，傅里叶变换可以揭示信号的重要特征，即其频率分量。
/ c: i2 r% N1 Q: M8 I) s


这个公式也不该陌生吧，就是DFT的公式，见下篇：终于到来的DFT

" K% r+ y/ c9 R! g7 u里面给出DFT的公式：
. X6 l7 I  M* J* e" O6 p( x. ^
: q2 i( g3 ?1 A" U
, X. E5 l' o& E1 l6 {
5 h1 x' K6 E9 L, ^5 l对比一下，你会发现其实二者是一样的。
. o8 N0 v5 e1 I6 q; t
. M9 a+ y& }& n- f- M5 `, zMATLAB®中的fft函数使用快速傅里叶变换算法来计算数据的傅里叶变换。 考虑一个正弦信号x，它是时间t的函数，频率分量为15 Hz和20 Hz。 使用以10秒为单位以1/50秒为增量采样的时间向量。
, @7 t* U0 Z6 l% q

8 }/ O  o4 U6 Q  Jt = 0:1/50:10-1/50;                     
$ x- J9 z' y) Y. _# O
* e2 a5 T2 M; G' L  C3 vx = sin(2*pi*15*t) + sin(2*pi*20*t);

plot(t,x)
# U9 ^0 g+ F$ S


计算信号的傅立叶变换，并创建对应于频率空间中信号采样的向量f。

- V6 \0 c9 T$ w5 _9 [( Ny = fft(x);     
7 M6 c" `# k( F3 s
" d0 a1 J8 n0 r) hf = (0:length(y)-1)*50/length(y);

, W% m* `' L5 B3 r5 ^" d/ U当您将信号幅度绘制为频率的函数时，幅度峰值对应于15 Hz和20 Hz的信号频率分量。

plot(f,abs(y))

, v$ ]: r/ u, e1 `) V! q& }3 J* Jtitle('Magnitude')

/ U" ]. D( I2 Q' b/ p  K1 \. {

变换还会生成尖峰的镜像副本，这些副本对应于信号的负频率。 为了更好地可视化此周期性，您可以使用fftshift函数，该函数在变换上执行零中心的循环移位。

n = length(x);                        
0 ~4 b  J9 G; O! R& C7 }- j
7 J$ G$ C" ^$ O6 T. kfshift = (-n/2:n/2-1)*(50/n);

yshift = fftshift(y);

plot(fshift,abs(yshift))


) g+ }7 n1 j( f, U5 l+ ?* O
* H0 j9 I! O+ `
Noisy Signals

在科学应用中，信号经常被随机噪声破坏，掩盖了它们的频率成分。 傅里叶变换可以处理随机噪声并显示频率。

6 C8 J1 c4 }: Q4 kFor example, create a new signal, xnoise, by injecting Gaussian noise into the original signal, x.

" J3 z1 x1 a8 D7 a" N5 v: Bxnoise = x + 2.5*gallery('normaldata',size(t),4);

+ ?/ ^9 a; J& R
+ X9 ?$ q0 G) t6 \/ N# g画出这个信号的图像：
; X8 f7 w& c* s9 V) ?
plot(t,xnoise);


) V3 p" x% C3 X# \
( u+ G; h3 ^7 |: h) A作为频率函数的信号功率是信号处理中使用的常用度量。 功率是信号傅立叶变换的平方幅度，通过频率样本的数量进行归一化。

6 _) v5 X+ Q3 Q' B: R计算并绘制以零频率为中心的噪声信号的功率谱。 尽管存在噪音，但由于功率峰值，您仍然可以确定信号的频率。
) E& s' v8 \0 O& ~: f$ }: r4 J. i, J
. f7 ~% O* z) k5 ]ynoise = fft(xnoise);

0 ^  m  p1 u9 S5 ?! r! Eynoiseshift = fftshift(ynoise);   

" b- {) Q' D3 p4 S0 Y" ~power = abs(ynoiseshift).^2/n;
. u8 U4 g9 o* L$ {
plot(fshift,power)

title('Power')


; c8 S* r7 ^6 N0 w2 s
Computational Efficiency

直接使用傅立叶变换公式计算浮点运算顺序中的每个需求元素。 快速傅里叶变换算法仅需要按计算的操作顺序。 当处理具有数百万个数据点的数据时，这种计算效率是一个很大的优势。 当功率为2时，快速傅里叶变换算法的许多专门实现甚至更有效。
: o/ I& v# B& U5 p; s考虑从加利福尼亚海岸外的水下麦克风收集的音频数据。 这些数据可以在康奈尔大学生物声学研究计划维护的图书馆中找到。 在bluewhale.au中加载和格式化数据的子集，其中包含太平洋蓝鲸发声。 您可以使用命令声音（x，fs）来收听整个音频文件。     
% x0 }5 P8 a5 [; u
whaleFile = 'bluewhale.au';

4 G  V8 f' @: N$ {( N& c; A1 s7 y# b[x,fs] = audioread(whaleFile);

whaleMoan = x(2.45e4:3.10e4);
4 g( @1 i! ~3 Y' N
t = 10*(0 :1/fs: (length(whaleMoan)-1)/fs);

: @; |$ C: C% p4 g: A. ^
* ~; V& z" \1 Z" ~0 u) R8 q
plot(t,whaleMoan)

0 x5 o, q, G, I. cxlabel('Time (seconds)')
9 |, C7 R' b6 ]% b: p
ylabel('Amplitude')
9 D1 `8 H7 ]/ ~0 K$ |6 D- ]3 c
xlim([0 t(end)])

/ }  Q; B/ m) G' W3 n5 w7 d! |

指定新的信号长度，该信号长度是2的下一个幂，大于原始长度。 然后，使用fft使用新的信号长度计算傅立叶变换。 fft使用零自动填充数据以增加样本大小。 这种填充可以使变换计算明显更快，特别是对于具有大质数因子的样本大小。
/ b, z! p% B- j: L
m = length(whaleMoan);

7 {0 I9 y! v9 T. U  {% v. en = pow2(nextpow2(m));
9 [0 g3 ~' ^9 ]* d+ Z) {6 m
. z2 W: `3 S$ x  w- [3 \y = fft(whaleMoan,n);        
6 S  i7 j; `8 [- i- l
1 g7 I' c8 i3 t' U- a绘制信号的功率谱。 该图表明，呻吟由大约17 Hz的基频和一系列谐波组成，其中强调了二次谐波。

f = (0:n-1)*(fs/n)/10; % frequency vector
  Q& J6 [5 E! }- W" D9 i
power = abs(y).^2/n;   % power spectrum      

- f6 u, i% U# Iplot(f(1:floor(n/2)),power(1:floor(n/2)))
1 h* M5 m; h/ t6 P1 a! T, m# H' g
xlabel('Frequency')
8 {8 x* y& O0 s! E# i
ylabel('Power')
4 y. Z) B8 R# S/ D/ J


/ E4 O, G' j: H
, O- ~7 l, s  _

yin123

MATLAB之Fourier Transforms ( fft )