Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析

1. 定义

信号在频域能够呈现出时域不易发现的性质和规律，傅里叶变换是将信号从时域变换到频域，便于在频域对信号的特性进行分析。离散傅里叶变换 (DFT)，是傅里叶变换在时域和频域上的离散呈现形式，通俗的说就是将经过采样的有限长度时域离散采样序列变换为等长度的频域离散采样序列，通过对变换得到的频域采样序列进行适当的换算和处理，可以得到信号的频谱（频率-幅值曲线和频率-相位曲线）。
离散傅里叶变换 (DFT)的定义为：

$X(k)=\mathrm{D}\mathrm{F}\mathrm{T}[x(n)]=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk}$, $0\le k\le N-1$

式中，$x(n)$为时域离散采样序列（通常为实数序列），$N$为时域离散采样序列$x(n)$的长度，$X(k)$为频域离散采样序列（通常为复数序列）。

快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的一种快速算法，FFT的计算结果与DFT完全相同，但FFT相对于DFT减小了计算量、节约计算资源消耗，能够适应在线计算，因此实际DFT都是通过FFT算法来求得结果。

2. 变换和处理

Matlab软件自带fft函数实现快速傅里变换算法，但是光使用fft并不能直接得到信号的频谱，还需要解决以下问题，下图为DFT变换后的X(k)复数序列幅值、相位图。

a) 幅值变换：$X(k)$序列的幅值大小与参与变换的时域序列x(n)长度N有关，变换后的幅值$|X(k)|$需要乘以$2/N$得到真实幅值；
b) 有效频率区域：$X(k)$序列由两部分共轭复数序列组成（复数共轭表示幅值相等、相位相反），相当于只有一半的复数序列是独立有效的，这部分复数序列对应$0-f_s/2$的频率区域（$f_s$为时域离散采样序列$x(n)$的采样频率）。
c) 直流信号的处理：直流信号幅值（对应频率0Hz）为两部分共轭复数序列在频率0Hz处的加和，其真实幅值再乘以$2/N$后还需要再除以2得到真实的直流信号幅值。

3. 函数

初学的朋友若不理解上述变换和处理技巧，很难得到正确的频谱图。为此作者在fft函数的基础上，使用Matlab开发了函数DFT.m，通过函数来实现上述幅值变换、有效频率区域和直流信号的处理，能够直接分析出给定离散信号x(n)的幅值谱和相位谱，函数简单、易用、通用性好。

function [f,X_m,X_phi] = DFT(xn,ts,N,drawflag)
% [f,X_m,X_phi] = DFT(xn,ts,N,drawflag) 离散序列的快速傅里叶变换，时域转换为频域
% 输入  xn为离散序列 为向量  
%       ts为序列的采样时间/s
%       N为FFT变换的点数，默认为xn的长度  
%       drawflag为绘图标识位，取0时不绘图，其余非0值时绘图，默认为绘图
% 输出 f为频率向量
%      X_m为幅值向量
%      X_phi为相位向量，单位为°
% 注意计算出来的0频分量(直流分量应该除以2)  直流分量的符号应结合相位图来确定
% By ZFS@wust  2020
% 获取更多Matlab/Simulink原创资料和程序，清关注微信公众号：Matlab Fans

4. 实例演示

下面结合实例进行演示和分析。

例1：单频正弦信号（整数周期采样）

%% Eg 1 单频正弦信号
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 2;         % 频率
w = 2*pi*f;    % 角频率
phi = pi/3;    % 初始相位 
x = A*cos(w*t+phi);   % 时域信号
figure
plot(t,x)
xlabel('时间/s')
ylabel('时域信号x(t)')
% DFT变换将时域转换到频域,并绘制频谱图
[f,X_m,X_phi] = DFT(x,ts);

结果

正弦信号频率为2Hz，频谱分析频率为1.98Hz
正弦信号幅值为1.5，频谱分析幅值为1.495
正弦信号相位为60°，频谱分析相位为63.32°

例2：单频正弦信号（非整数周期采样）

%% Eg 2 单频正弦信号(非整数周期采样)
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 1.5;         % 频率
w = 2*pi*f;    % 角频率
phi = pi/3;    % 初始相位 
x = A*cos(w*t+phi);   % 时域信号
figure
plot(t,x)
xlabel('时间/s')
ylabel('时域信号x(t)')
% DFT变换将时域转换到频域,并绘制频谱图
[f,X_m,X_phi] = DFT(x,ts);

结果

正弦信号频率为1.5Hz，频谱分析频率为0.99Hz、1.98Hz
正弦信号幅值为1.5，频谱分析幅值为1.034、0.923
正弦信号相位为60°，频谱分析相位为160.93°、-33.76°

总结：DFT变换后频率序列的最小单位刻度为$f_s/N$（此例为1Hz），非整数周期采样时关心信号的频率（此例为1.5Hz）不是频率分辨率$f_s/N$的正整数倍，那这个频率成分信号会由前后两个正整数倍的频率成分信号（此例为1Hz和2Hz）的线性组合来替代，这就是频谱泄漏现象，非周期采样时某频率成分信号向两侧频率分辨率正整数倍的频点泄漏。实际频谱分析时并不清楚所关心的频率点精确值，避免此问题的一个解决方法是，取更多的点参加DFT，即时域序列$x(n)$长度$N$值取长一些，让频率分辨率$f_s/N$很小，以减小频谱泄漏现象。

%% Eg 2 单频正弦信号(非整数周期采样)
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 1.5;         % 频率
w = 2*pi*f;    % 角频率
phi = pi/3;    % 初始相位 
x = A*cos(w*t+phi);   % 时域信号
figure
plot(t,x)
xlabel('时间/s')
ylabel('时域信号x(t)')
% DFT变换将时域转换到频域,并绘制频谱图
[f,X_m,X_phi] = DFT(x,ts);

结果：频谱泄漏情况大为改善，采样点继续增多时，频谱泄漏会进一步减小。

正弦信号频率为1.5Hz，频谱分析频率主要成分为1.46Hz、
正弦信号幅值为1.5，频谱分析频率主要成分对应幅值为1.41
正弦信号相位为60°，频谱分析频率主要成分对应相位为89.5°

例3：含有直流分量的单频正弦信号

%% Eg 3 含有直流分量的单频正弦信号
ts = 0.01;
t = 0:ts:1;
A = 1.5;       % 幅值  
f = 5;         % 频率
w = 2*pi*f;    % 角频率
phi = pi/6;    % 初始相位 
x = 0.5 + A*cos(w*t+phi);   % 时域信号,带有直流偏移0.5
figure
plot(t,x)
xlabel('时间/s')
ylabel('时域信号x(t)')
% DFT变换将时域转换到频域,并绘制频谱图
[f,X_m,X_phi] = DFT(x,ts);

结果

正弦信号频率为5Hz，频谱分析频率为4.95Hz
正弦信号幅值为1.5，频谱分析幅值为1.498
正弦信号相位为30°，频谱分析相位为38.66°
正弦信号直流分量0.5，频谱分析直流分量为0.51

例4：正弦复合信号

%% Eg 4 正弦复合信号
ts = 0.01;
t = 0:ts:2;
A = [1.5 1 0.5 0.2];    % 幅值  
f = [3 6 9 15];         % 频率
w = 2*pi*f;             % 角频率
phi = (1:4)*pi/4;       % 初始相位 
x = -0.5 + A(1)*cos(w(1)*t+phi(1)) + A(2)*cos(w(2)*t+phi(2)) + A(3)*cos(w(3)*t+phi(3)) + A(4)*cos(w(4)*t+phi(4));     % 时域信号
figure
plot(t,x)
xlabel('时间/s')
ylabel('时域信号x(t)')
% DFT变换将时域转换到频域,并绘制频谱图
[f,X_m,X_phi] = DFT(x,ts);

结果：

正弦信号频率为3、6、9、15Hz，频谱分析频率为2.985、5.97、8.96、14.93Hz
正弦信号幅值为1.5、1、0.5、0.2，频谱分析幅值为1.499、0.989、0.485、0.192
正弦信号相位为45°、90°、135°、180°，频谱分析相位为50.6°、101.5°、152.9°、210°
正弦信号直流分量-0.5，频谱分析直流分量为-0.497

注意：频率为0Hz时对应的直流信号的幅值的正负号，是通过零频相位来确定的，相位为0°表示幅值为正，相位为180°表示幅值为负。

例5：含有随机干扰的正弦信号

%% Eg 5 含有随机干扰的正弦信号
ts = 0.01;
t = 0:ts:2;
A = [1 0.5];    % 幅值  
f = [3 10];         % 频率
w = 2*pi*f;             % 角频率
phi = (1:2)*pi/3;       % 初始相位 
x =  A(1)*cos(w(1)*t+phi(1)) + A(2)*cos(w(2)*t+phi(2)) + 0.8*(rand(size(t))-0.5);     % 时域信号
figure
plot(t,x)
xlabel('时间/s')
ylabel('时域信号x(t)')
% DFT变换将时域转换到频域,并绘制频谱图
[f,X_m,X_phi] = DFT(x,ts);

结果：

正弦信号频率为3、10Hz，频谱分析频率为2.985、9.95Hz
正弦信号幅值为1、0.5，频谱分析幅值为0.978、0.456
正弦信号相位为60°、135°，频谱分析相位为65.1°、139.8°

例6：实际案例

load data
ts = 0.001;
x = Jsd;
t = [0:length(x)-1]*ts;
figure
plot(t,x)
xlabel('时间/s')
ylabel('时域信号x(t)')
% DFT变换将时域转换到频域,并绘制频谱图
[f,X_m,X_phi] = DFT(x,ts);

结果：

频谱分析主要频率成分为18.996、37.992Hz
频谱分析主要频率成分对应幅值为1.741、1.117

该项目为作者在强振环境下测得加速度信号，加速度是机械结构周期运动激励产生，需要通过频谱分析获取机械结构周期运动的频率。由于噪声幅度远大于有效信号幅度，信号的信噪比很低，从时域上很难辨别机械结构周期运动的频率。但经过DFT后，从频域上可以看出信号的主要频率成分为19Hz和其倍频38Hz，可以判断机械结构周期运动的频率为19Hz，38Hz为结构响应的非线性特性所产生的倍频。

5. 拓展

工程上我们还会遇到这样的问题：获取了信号的频谱，希望从信号的频谱来恢复时域信号。这就涉及离散傅里叶逆变换问题，下一篇详谈。

6. 联系作者

有Matlab/Simulink方面的技术问题，欢迎发送邮件至944077462@qq.com讨论。
更多Matlab/Simulink原创资料，欢迎关注微信公众号：Matlab Fans

源程序下载:
1. csdn资源: Matlab如何进行离散傅里叶变换DFT(快速傅里叶变换FFT)进行频谱分析
2. 扫码关注微信公众号Matlab Fans，回复BK07获取百度网盘下载链接。