| Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析 | 您所在的位置：网站首页 › matlab如何调用函数得到复指数序列 › Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析 |
Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析
|
文章目录
1. 定义2. 变换和处理3. 函数4. 实例演示例1：单频正弦信号（整数周期采样）例2：单频正弦信号（非整数周期采样）例3：含有直流分量的单频正弦信号例4：正弦复合信号例5：含有随机干扰的正弦信号例6：实际案例
5. 拓展6. 联系作者
1. 定义
信号在频域能够呈现出时域不易发现的性质和规律，傅里叶变换是将信号从时域变换到频域，便于在频域对信号的特性进行分析。离散傅里叶变换 (DFT)，是傅里叶变换在时域和频域上的离散呈现形式，通俗的说就是将经过采样的有限长度时域离散采样序列变换为等长度的频域离散采样序列，通过对变换得到的频域采样序列进行适当的换算和处理，可以得到信号的频谱（频率-幅值曲线和频率-相位曲线）。 离散傅里叶变换 (DFT)的定义为： X ( k ) = D F T [ x ( n ) ] = ∑ n = 0 N − 1 x ( n ) e − j 2 π N n k X(k) = {\rm{DFT}}[x(n)] = \sum\limits_{n = 0}^{N - 1} {x(n){e^{ - j\frac{{2\pi }}{N}nk}}} X(k)=DFT[x(n)]=n=0∑N−1x(n)e−jN2πnk, 0 ≤ k ≤ N − 1 0 \le k \le N-1 0≤k≤N−1 式中， x ( n ) x(n) x(n)为时域离散采样序列（通常为实数序列）， N N N为时域离散采样序列 x ( n ) x(n) x(n)的长度， X ( k ) X(k) X(k)为频域离散采样序列（通常为复数序列）。 快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的一种快速算法，FFT的计算结果与DFT完全相同，但FFT相对于DFT减小了计算量、节约计算资源消耗，能够适应在线计算，因此实际DFT都是通过FFT算法来求得结果。 2. 变换和处理Matlab软件自带fft函数实现快速傅里变换算法，但是光使用fft并不能直接得到信号的频谱，还需要解决以下问题，下图为DFT变换后的X(k)复数序列幅值、相位图。
a) 幅值变换： X ( k ) X(k) X(k)序列的幅值大小与参与变换的时域序列x(n)长度N有关，变换后的幅值 ∣ X ( k ) ∣ |X(k)| ∣X(k)∣需要乘以 2 / N 2/N 2/N得到真实幅值； b) 有效频率区域： X ( k ) X(k) X(k)序列由两部分共轭复数序列组成（复数共轭表示幅值相等、相位相反），相当于只有一半的复数序列是独立有效的，这部分复数序列对应 0 − f s / 2 0 - f_s/2 0−fs/2的频率区域（ f s f_s fs为时域离散采样序列 x ( n ) x(n) x(n)的采样频率）。 c) 直流信号的处理：直流信号幅值（对应频率0Hz）为两部分共轭复数序列在频率0Hz处的加和，其真实幅值再乘以 2 / N 2/N 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);结果 结果 总结：DFT变换后频率序列的最小单位刻度为 f s / N f_s/N fs/N（此例为1Hz），非整数周期采样时关心信号的频率（此例为1.5Hz）不是频率分辨率 f s / N f_s/N fs/N的正整数倍，那这个频率成分信号会由前后两个正整数倍的频率成分信号（此例为1Hz和2Hz）的线性组合来替代，这就是频谱泄漏现象，非周期采样时某频率成分信号向两侧频率分辨率正整数倍的频点泄漏。实际频谱分析时并不清楚所关心的频率点精确值，避免此问题的一个解决方法是，取更多的点参加DFT，即时域序列 x ( n ) x(n) x(n)长度 N N N值取长一些，让频率分辨率 f s / N f_s/N fs/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);结果：频谱泄漏情况大为改善，采样点继续增多时，频谱泄漏会进一步减小。 结果 结果： 注意：频率为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);结果： 结果： 该项目为作者在强振环境下测得加速度信号，加速度是机械结构周期运动激励产生，需要通过频谱分析获取机械结构周期运动的频率。由于噪声幅度远大于有效信号幅度，信号的信噪比很低，从时域上很难辨别机械结构周期运动的频率。但经过DFT后，从频域上可以看出信号的主要频率成分为19Hz和其倍频38Hz，可以判断机械结构周期运动的频率为19Hz，38Hz为结构响应的非线性特性所产生的倍频。 5. 拓展工程上我们还会遇到这样的问题：获取了信号的频谱，希望从信号的频谱来恢复时域信号。这就涉及离散傅里叶逆变换问题，下一篇详谈。 6. 联系作者有Matlab/Simulink方面的技术问题，欢迎发送邮件至[email protected]讨论。 更多Matlab/Simulink原创资料，欢迎关注微信公众号：Matlab Fans 源程序下载: 1. csdn资源: Matlab如何进行离散傅里叶变换DFT(快速傅里叶变换FFT)进行频谱分析 2. 扫码关注微信公众号Matlab Fans，回复BK07获取百度网盘下载链接。 |
【本文地址】