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);
结果
正弦信号频率为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 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);
结果:频谱泄漏情况大为改善,采样点继续增多时,频谱泄漏会进一步减小。
正弦信号频率为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获取百度网盘下载链接。
Matlab如何进行利用离散傅里叶变换DFT (快速傅里叶变换FFT)进行频谱分析相关推荐
- 【从FT到DFT和FFT】(三)从离散傅里叶变换到快速傅里叶变换
文章目录 推荐阅读 前言 从离散傅里叶变换到快速傅里叶变换 单位根 对DFT进行分治得到FFT 计算前半截 计算后半截 快速傅里叶逆变换(IFFT) 推荐阅读 前置阅读 [从FT到DFT和FFT](一 ...
- 图像算法四:【图像增强--频率域】傅里叶变换、快速傅里叶变换、频域滤波、频域低通滤波、频域高通滤波
频率域滤波与空间域滤波殊途同归,空间域图像增强与频率域图像增强是两种截然不同的技术,实际上在相当程度上说它们是在不同的领域做相同的事情,只是有些滤波更适合在空间域完成,而有些则更适合在频率域中完成. ...
- 傅里叶变换与快速傅里叶变换
傅里叶变换与快速傅里叶变换 作为电子信息专业的学生老说,这个不知道,或者理解不清楚,是十分不应该的,作为一个学渣,有时候确实是理解不清楚的 1.首先离散傅里叶变换目的: 简单点说: 就是将一个信号从时 ...
- 卷积、傅里叶级数、傅里叶变换、快速傅里叶变换、pytorch中的fft,rfft
卷积: 连续形式: 离散形式: '卷' : 翻转 和 滑动 '积' : 积分 翻转:g(t) - > g(-t) 滑动:g(-t) - > g(n-t) 平移n个单位 举个例 ...
- 傅里叶变换 一维快速傅里叶变换(快速的一维离散傅里叶变换、分治法)
一.介绍 1.一维离散傅里叶变换DFT. DFT:(Discrete Fourier Transform)离散傅里叶变换是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域 ...
- MATLAB之傅里叶变换,快速傅里叶变换FFT
文章目录 傅里叶变换及傅里叶逆变换定义 窗函数/矩形脉冲信号的傅里叶变换 基于MATLAB的快速傅里叶变换FFT 傅里叶变换及傅里叶逆变换定义 能从时域的非周期连续信号转化到频域非周期连续信号. 窗函 ...
- matlab图片快速傅里叶变换,图像傅里叶变换(快速傅里叶变换FFT) | 学步园
#include "Image_FFT.h" /* 中心化,根据傅里叶性质的平移性质 */ void FFT_Shift(double * src,int size_w,int s ...
- 快速傅里叶变换python_快速傅里叶变换及python代码实现
一.前言 我想认真写好快速傅里叶变换(Fast Fourier Transform,FFT),所以这篇文章会由浅到细,由窄到宽的讲解,但是傅里叶变换对于寻常人并不是很容易理解的,所以对于基础不牢的人我 ...
- 图像傅里叶变换(快速傅里叶变换FFT)
学习DIP第7天,图像傅里叶变换 转载请标明出处:http://blog.csdn.net/tonyshengtan,欢迎大家转载,发现博客被某些论坛转载后,图像无法正常显示,无法正常表达本人观点,对 ...
最新文章
- WebForm 【上传图片--添加水印】
- 【高级绘图】MATLAB应用实战系列(八十)-圣诞前夜,想表白女神?教你如何用MATLAB绘制圣诞树动态图(附MATLAB代码)
- 07.MyBatis中的关联查询
- 数据仓库之电商数仓-- 2、业务数据采集平台
- windows installer没有正确安装_电脑还可以这样禁止软件自动安装,后悔知道得太晚...
- GossipView:圆圈布局的自定义view
- 使用PowerShell将字符串拆分为数组
- windbg查询内存泄笔记
- 部署exchange2010三合一:之二:先决条件
- sql函数 StringSplit(SELECT * from Split('黄色,蓝色,黑色',','))
- android ConstraintLayout布局 详解
- 英语听力采用计算机化考试,高考英语听力机考特点与应对建议
- Complier Validation via Equivalence Modulo Inputs
- 中职一年级计算机英语课件,职高一年级英语期中试题
- 莫圣宏:4.30黄金开启跌势,黑色星期五黄金操作建议!
- jenkins windows 20008 R2 msi 工作目录迁移
- Mybatis分页实现
- PRODUCT_CHARACTERISTICS 详解
- 内网ip映射到外网(路由器是tp-link)
- word中MathType公式不能 二次编辑解决方案