MATLAB：快速傅里叶变换（FFT）

快速傅里叶变换介绍

傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的余弦（或正弦）波信号的无限叠加。FFT 是离散傅立叶变换的快速算法，可以将一个信号变换到频域。那其在实际应用中，有哪些用途呢？

1. 有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征（频率，幅值，初相位）；

2. FFT 可以将一个信号的频谱提取出来，进行频谱分析，为后续滤波准备；

3. 通过对一个系统的输入信号和输出信号进行快速傅里叶变换后，两者进行对比，对系统可以有一个初步认识。

假设采样频率 $Fs$，信号频率 $F$，信号长度 $L$，采样点数 $N$。那么 FFT 之后结果就是一个为 $N$ 点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。

具体跟原始信号的幅度有什么关系呢？

1. 假设原始信号的峰值为 $A$，那么 FFT 的结果的每个点（除了第一个点直流分量之外）的模值就是 $A$ 的 $N/2$ 倍，而第一个点就是直流分量（即 0 Hz），它的模值是直流分量的 $N$ 倍；

2. 每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量，它的相位是该频率的初相位，MATLAB 以 $cos$ 为底的，若信号时正弦形式 $sin(t)$，则变成 $cos(t-\pi /2)$ 即可。

采样频率 $Fs$，被 $N-1$ 个点平均分成 $N$ 等份，每个点的频率依次增加。为了方便进行 FFT 运算，通常 $N$ 取大于信号长度 L 的 2 的整数次方。

例如某点 $n$ 所表示的频率为：$Fn=(n-1)\ast Fs/N$ 。由上面的公式可以看出，$Fn$ 所能分辨到频率为为 $Fs/N$。如果采样频率 $Fs$ 为 1024Hz，采样点数为 1024 点，则可以分辨到 1Hz。

1024Hz 的采样率采样 1024 点，刚好是 1 秒，也就是说，采样 1 秒时间的信号并做 FFT，则结果可以分析到 1Hz。如果采样2秒时间的信号，则 N 为2048，并做 FFT，则结果可以分析到 0.5Hz。

如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。

假设 FFT 之后某点 $n$ 用复数 $a+bi$ 表示，该复数的模就是$An=sqrt(a\ast a+b\ast b)$，相位就是 $Pn=atan2(b,a)$。根据以上的结果，就可以计算出 n 点（n≠1，且n<=N/2）对应的信号的表达式为：$A_n/(N/2)\ast cos(2\pi \ast Fn\ast t+Pn)$，即$2\ast An/N\ast cos(2\ast pi\ast Fn\ast t+Pn)$；对于 $n=1$ 点的信号，是直流分量，幅度即为 $A_1/N$。

由于 FFT 结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。

例子

假设我们有一个信号，它含有 5V 的直流分量，频率为 50Hz、相位为 -30 度、幅度为 7V 的交流信号以及一个频率为 90Hz、相位为 90度、幅度为 3V 的交流信号。数学表达式为：

$$x=5+7cos(2\pi \ast 15\ast t-30\pi /180)+3cos(2\pi \ast 40\ast t-90\pi /180)$$

我们以 128Hz 的采样率对这个信号进行采样，总共采样 256 点。按照我们上面的分析，$Fn=(n-1)\ast Fs/N$，我们可以知道，每两个点之间的间距就是 0.5Hz。我们的信号有3个频率：0Hz、15Hz、40Hz

出于编程方便，因为直流分量的幅值 $A_1/N$，其他点幅值为 $A_n/(N/2)$，故直流分量最后要除以 2 才是对的。

一般 FFT 所用数据点数 $N$ 与原含有信号数据点数 $L$ 相同，这样的频谱图具有较高的质量，可减小因补零或截断而产生的影响。

Matlab代码

clear
Fs = 128;       % 采样频率
T = 1/Fs;       % 采样时间
L = 256;        % 信号长度
t = (0:L-1)*T; % 时间
x = 5 + 7*cos(2*pi*15*t - 30*pi/180) + 3*cos(2*pi*40*t - 90*pi/180);   %cos为底原始信号
y = x + randn(size(t));     %添加噪声

figure;
plot(t,y)
title('加噪声的信号')
xlabel('时间(s)')
N = 2^nextpow2(L); %采样点数，采样点数越大，分辨的频率越精确，N>=L，超出的部分信号补为0
Y = fft(y,N)/N*2;   %除以N乘以2才是真实幅值，N越大，幅值精度越高
f = Fs/N*(0:1:N-1); %频率
A = abs(Y);     %幅值
P = angle(Y);   %相值
figure;
subplot(211);plot(f(1:N/2),A(1:N/2));   %函数fft返回值的数据结构具有对称性,因此我们只取前一半
title('幅值频谱')
xlabel('频率(Hz)')
ylabel('幅值')
subplot(212);plot(f(1:N/2),P(1:N/2));
title('相位谱频')
xlabel('频率(Hz)')
ylabel('相位')

参考：
https://blog.csdn.net/deng_sai/article/details/52473503