信号处理快速傅里叶变换(FFT)的学习

FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。采样频率要大于信号频率的两倍。为了方便进行FFT运算,通常N取2的整数次方。

一、DFT

1.1 原理

定义(来自百科):离散傅里叶变换(Discrete Fourier Transform,缩写为DFT),是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域采样。在形式上,变换两端(时域和频域上)的序列是有限长的,而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT,也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT。

DFT的公式是:设x(n)为M点的有限长序列,即在0≤n≤M-1内有值,则定义x(n)的N点(N≥M。当N>M时,补N-M个零值点),离散傅里叶变化定义为

其中,为旋转因子(为方便编辑后续记作W(N, nk),特此说明),其计算公式为,具有以下性质:

①周期性:W(N, nk) = W(N, (n+rN)k) = W(N, n(k+rN)),其中r问整数。

②共轭对称性:(W(N, nk))* = W(N, -nk)。

③可约性:W(N, nk) = W(mN, mnk),W(N, nk) = W(N/m, nk/m),其中m为整数,N/m为整数。

④特殊值:W(N, N/2) = -1;  W(N, (k+N/2)) = -W(N, k);  

                  W(N, (N-k)n) = W(N, (N-n)k) = W(N, -nk)。

所以,在计算旋转因子的过程中可以适当的使用特殊值来提高运算的效率。

在计算DFT时,如果数据的点数够计算的点数,则截取计算点数长的数据进行DFT运算,否则将数据点个数补0至计算点个数,然后进行计算,举例如下(举例只是为了说明问题,没有按照编程语言的书写格式)。

比如:数据点数组为 dataArray = {1, 2, 3, 4,5},可以看出数据长度为5。

如果要求做4点DFT运算,则只需截取前4个数作为运算数组进行运算即可,即为{1, 2, 3, 4};

如果要求做8点DFT运算,则需在原数组后补三个0,使长度为8后再进行计算,计算数组为{1, 2, 3, 4,5,0,0,0}。

1.2 流程图

二、FFT

2.1 原理

设序列x(n)点数为N=2^L,L为整数。将N=2^L,先按照n的奇偶分为两组,

其中r = 0, 1, ..., N/2-1,x(2r) = x1(r),x(2r-1) = x2(r);

则可将DFT化为

由上式可以看出一个N点的DFT可以分为两个N/2点的DFT,按照上式右组合成N点DFT。但是这里的x1(r)、x2(r)以及X1(k)、X2(k)都是N/2点的序列,即r,k满足r,k=0, 1, ..., N/2-1。而X(k)却有N点,上式计算的只是X(k)的前半项数的结果,因此还需要计算后半项的值。

将①②③代入(式1),就可将X(k)表达为前后两部分:

前半部分,X(k),当k=0, 1, ..., N/2-1

后半部分,X(k),当k=N/2, ..., N-1

2.2 流程图

①FFT运算总流程图,包括

②根据(1)中推导的内容,FFT的流程图可化为

蝶形运算中又三层循环

第一层(最外层),完成M次迭代过程,即算出A0(k), A1(k), ..., Am(k),其中k=0, 1, ..., N;A2(k)为蝶形运算第2级的结果,如A0(k)=x(k), Am(k)=X(k);

第二层(中间层),完成旋转因子的变化,步进为2^L;

第三层(最里层),完成相同旋转因子的蝶形运算。

三、DFT与FFT

1.DFT的运算量

复数乘法次数为N*N,复数加法次数为N(N-1),若N>>1,则这两者都近似为N*N,它随N增大为急速增大。改进途径:①利用旋转因子性质减小计算量;②由于运算量和N*N成正比,因而可将N点DFT分解成小点数的DFT,以减少运算量(点数越小,计算量越小);③改为用FFT计算,复数乘法次数为(N/2)*log(2,N)。

2.FFT可分为按照时间抽选的基-2算法(库利-图基算法DIT-FFT)和按频率抽选的基-2算法(桑德-图基算法DIF-FFT)。

3.FFT计算原理

FFT的计算要求点数必须为2的整数次幂,如果点数不够用0补齐。例如计算{2,3,5,8,4}的8点FFT,需要补3个0后进行计算,如果计算该数组的5点FFT,则先计算8点FFT后截取前5个值即可(不提倡)。

四、蝴蝶运算

4.1 基本流程

蝶形运算,符号表示如下图所示:

所以,FFT的蝶形运算图表示为(8点,2^3 = 8,运算级数最大为L=3)

在蝶形运算中变化规律由W(N, p)推导,其中N为FFT计算点数,J为下角标的值

        L = 1时,W(N, p) = W(N, J) = W(2^L, J),其中J = 0;

        L = 2时,W(N, p) = W(N, J) = W(2^L, J),其中J = 0, 1;

        L = 3时,W(N, p) = W(N, J) = W(2^L, J),其中J = 0, 1, 2, 3;

所以,W(N, p) = W(2^L, J),其中J = 0, 1, ..., 2^(L-1)-1

又因为2^L = 2^M*2^(L-M) = N*2^(L-M),这里N为2的整数次幂,即N=2^M,

W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))

所以,p = J*2^(M-L),此处J = 0, 1, ..., 2^(L-1)-1,当J遍历结束但计算点数不够N时,J=J+2^L,后继续遍历,直到计算点数为N时不再循环。

4.2 举例:N=8点的FFT计算

当L=2时,J = 0, 1两个值,因此p = J*2^(M-L) = 0, 2两个值,即旋转因子有两个值W(8, 0)和W(8, 2),计算中两行之间的距离B = 2^(L-1)=2^(2-1)=2,

代入J=0, 1可求得X(0)、X(0+2)和X(1)、X(1+2),即可求出第2级蝶形运算的X(0), X(1), X(2), X(3),也就是求出一半,此时J加步进2^L=4,即J=J+2^L=4, 5,

再代入J=4, 5可求出X(4), X(4+2)和X(5), X(5+2),即可求出第2级蝶形运算的X(4), X(5), X(6), X(7),已经全部求出,J循环结束。

4.3 二进制倒序说明

前面数的排列顺序是进行二进制倒序后的排序。二进制倒序是指将某数转化为二进制表示,将最高位看做最低位、次高位看做次低位...以此类推,计算后的纸进行排序。例如3的二进制表示为011b(3=0*2^2+1*2+1),二进制倒序值为6(011b,最高位看做最低位...即6=0+1*2+1*2^2)

即0到7的二进制排序是:

       0,                4,                 2,                6,                 1,                5,                3,                7

000→000    001→100    010→010    011→110    100→001    101→101    110→011    111→111

五、从振幅或快速傅立叶变换到dB

从振幅或快速傅立叶变换到dB是一个涉及信号处理和单位转换的问题。

函数f(t)为一元连续函数,其傅里叶变换定义为:

其中,i为虚数单位。欧拉公式:

5.1 相关名词

1.振幅(Amplitude):振幅是指信号的最大偏离值,表示信号的强度或能量大小。在信号处理中,振幅常用于描述波形的高低或振动的幅度。

2.dB(分贝):分贝是一种用于表示信号强度或功率比例的对数单位。在信号处理中,分贝常用于描述信号的增益或衰减程度。分贝的计算公式为:dB = 10 * log10(P1/P0),其中P1为信号的功率,P0为参考功率。当保持输入信号的幅度不变,改变频率使输出信号降至最大值的0.707倍,即用频响特性来表述即为-3dB点处即为截止频率,它是用来说明频率特性指标的一个特殊频率。

利用系统函数的模来表示电路的放大倍数,由于20lgA(ω)=-3dB,解得

A(ω)=10^-0.15=0.707945784≈1/√2,又因为A(ω)=|H(jω)|,则|H(jω)|^2=1/2

在高频端和低频端各有一个截止频率,分别称为上截止频率和下截止频率。两个截止频率之间的频率范围称为通频带。

5.2 物理意义理解

假设采样频率为Fs,信号频率F,采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。

1.幅值:这个点的模值,就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢?假设原始信号的峰值为A,那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。而第一个点就是直流分量,它的模值就是直流分量的N倍。

2.相位:而每个点的相位呢,就是在该频率下的信号的相位。

3.分辨率:第一个点表示直流分量(即0Hz),而最后一个点N的再下一个点(实际上这个点是不存在的,这里是假设的第N+1个点,也可以看做是将第一个点分做两半分,另一半移到最后)则表示采样频率Fs,这中间被N-1个点平均分成N等份,每个点的频率依次增加。例如某点n所表示的频率为:Fn=(n-1)*Fs/N。由上面的公式可以看出,Fn所能分辨到频率为为Fs/N,如果采样频率Fs为1024Hz,采样点数为1024点,则可以分辨到1Hz。1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时间的信号并做FFT,则结果可以分析到0.5Hz。如果要提高频率分辨力,则必须增加采样点数,也即采样时间。频率分辨率和采样时间是倒数关系。

4.表达式:假设FFT之后某点n用复数a+bi表示,那么这个复数的模就是An=根号a*a+b*b,相位就是Pn=atan2(b,a)。根据以上的结果,就可以计算出n点(n≠1,且n<=N/2)对应的信号的表达式为:An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号,是直流分量,幅度即为A1/N。

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

5.3 示例

假设我们有一个信号,它含有2V的直流分量,频率为50Hz、相位为-30度、幅度为3V的交流信号,以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下(式中cos参数为弧度,所以-30度和90度要分别换算成弧度):

S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

我们以256Hz的采样率对这个信号进行采样,总共采样256点。按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个点之间的间距就是1Hz,第n个点的频率就是n-1。我们的信号有3个频率:0Hz、50Hz、75Hz,应该分别在第1个点、第50个点、第76个点上出现峰值,其它各点应该接近0。

幅值 相位(弧度)相位(角度)
直流分量

512/N

=512/256

=2

\\
50Hz

384/(N/2)

=384/(256/2)

=3

atan2(-192, 332.55)

=-0.5236

(-0.5236)/pi

=-30.0001

75Hz

192/(N/2)

=192/(256/2)

=1.5

atan2(192, 3.4315E-12)

=1.5708

180*1.5708/pi

=90.0002

FFT结果的物理意义

(图片横坐标为第几个点,即Hz;纵坐标为模)

5.4 波特图

1、伯德图的定义

伯德图由对数幅频特性和对数相频特性两条曲线构成。 伯德图一般是由二张图组合而成, 伯德图由两张图组成:①G(jω)的幅值(以分贝,dB表示)-频率(以对数标度)对数坐标图,其上画有对数幅频曲线;②G(jω)的相角-频率(以对数标度)对数坐标图,其上画有相频曲线。 

对数幅值的标准表达式为20 lg|G(jω)|,单位是分贝,相角的单位是度。由于增益用对数来表示(log(ab)=log(a)+log(b)),因此一传递函数乘以一常数,在伯德增益图只需将图形的纵向移动即可,二传递函数的相乘,在波德幅频图就变成图形的相加。幅频图纵轴0分贝以下具有正增益裕度、属稳定区,反之属不稳定区。

配合波德相频图可以估算一信号进入系统后,输出信号及原始信号的比例关系及相位。例如一个Asin(ωt) 的信号进入系统后振幅变原来的k倍,相位落后原信号Φ,则其输出信号则为(Ak)sin(ωt−Φ),其中的k和Φ都是频率的函数。相频图纵轴-180度以上具有正相位裕度、属稳定区,反之属不稳定区。

2、伯德图的画法

画伯德图时,分三个频段进行,先画幅频特性,顺序是中频段、低频段和高频段。将三个频段的频率特性(或称频率响应)合起来就是全频段的幅频特性,然后再根据幅频特性画出相应的相频特性来。 作伯德图时,首先写出频率特性,然后按常数因子K、积分和微分因子(jω)、一阶因子(1+jωT)和二阶因子[1+2ζ(jω/ωn)+(jω)/ω]1这样四种基本因子分别画出伯德图,再总加而成。

下图是常见简易传递函数的伯德图趋势:

五、代码实现(C语言)

FFT.H

FFT.C

六、代码实现(Matlab)

基于MATLAB/Simulink的2ASK数字带通传输系统建模与仿真

相关文章

十分简明易懂的FFT(快速傅里叶变换)

彻底搞懂快速傅里叶变换FFT核心思想--分而治之

数字信号处理公式变程序(一)——DFT、FFT

实验笔记【3】| 利用matlab进行fft分析

信号处理之快速傅里叶变换(FFT)-通俗易懂

深入浅出的理解频谱泄露

快速傅里叶变换(FFT)和小波分析在信号处理上的应用

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/435911.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

leetcode每日一题day19(24.9.29)——买票需要的时间

思路&#xff1a;在最开始的情况下每人需要买的票数减一是能保持相对位置不变的&#xff0c; 如果再想减一就有可能 有某些人只买一张票&#xff0c;而离开了队伍&#xff0c; 所有容易想到对于某个人如果比当前的人买的多就按当前的人数量算 因为在一次次减一的情况下&#xf…

从零开始手写STL库:Stack

从零开始手写STL库–Stack的实现 Gihub链接&#xff1a;miniSTL 文章目录 从零开始手写STL库–Stack的实现一、stack是什么&#xff1f;二、stack要包含什么函数总结 一、stack是什么&#xff1f; 栈是一种后进先出&#xff08;LIFO&#xff0c;Last In First Out&#xff09…

计算机网络--TCP、UDP抓包分析实验

计算机网络实验 目录 实验目的 实验环境 实验原理 1、UDP协议 2、TCP协议 实验具体步骤 实验目的 1、掌握使用wireshark工具对UDP协议进行抓包分析的方法&#xff0c;掌握UDP协议的报文格式&#xff0c;掌握UDP协议校验和的计算方法&#xff0c;理解UDP协议的优缺点&am…

【数据库文档】数据库设计说明书(Word原件参考)

一、 总述 &#xff08;一&#xff09; 编写目的 二、 外部设计 &#xff08;一&#xff09; 环境说明 &#xff08;二&#xff09; 指导 三、 物理实现 &#xff08;一&#xff09; 物理结构 &#xff08;二&#xff09; 安全设计 四、 表设计结构 &#xff08;一&#xff09;…

SpringAOP学习

面向切面编程&#xff0c;指导开发者如何组织程序结构 增强原始设计的功能 oop:面向对象编程 1.导入aop相关坐标&#xff0c;创建 <!--spring依赖--><dependencies><dependency><groupId>org.springframework</groupId><artifactId>spri…

Python 读取与处理出入库 Excel 数据实战案例(HTML 网页展示)

有如下数据&#xff0c;需要对数据合并处理&#xff0c;输出到数据库。 数据样例&#xff1a;&#x1f447; excel内容&#xff1a; 出入库统计表河北库.xlsx: 出入库统计表天津库.xlsx: 01实现过程 1、创建test.py文件&#xff0c;然后将下面代码复制到里面&#xff0c;最后…

广西容县霞烟鸡,品牌兴农,助力乡村振兴!

在两广与港澳地区,流传着一句深入人心的饮食谚语——“无鸡不成宴”,它不仅是一种习俗的体现,更是对餐桌礼仪与待客之道的深刻诠释。每逢家宴欢聚、祭祖庆典或盛宴宾客,一只精心烹制的鸡总是不可或缺的主角,其缺席往往被视为对宾客的不敬。在这片美食文化的沃土上,广西容县的霞…

个性化大语言模型:PPlug——让AI更懂你

在当今数字化转型的时代&#xff0c;大型语言模型&#xff08;LLMs&#xff09;已经成为了不可或缺的工具&#xff0c;它们在自然语言理解、生成和推理方面展现了非凡的能力。然而&#xff0c;这些模型普遍采用的是“一刀切”的方式&#xff0c;即对于相同的输入给予所有用户相…

uniapp监听滚动实现顶部透明度变化

效果如图&#xff1a; 实现思路&#xff1a; 1、使用onPageScroll监听页面滚动&#xff0c;改变导航条的透明度&#xff1b; 2、关于顶部图片的高度&#xff1a; 如果是小程序&#xff1a;使用getMenuButtonBoundingClientRect获取胶囊顶部距离和胶囊高度&#xff1b; 如果…

YOLOv8改进 - 注意力篇 - 引入SCAM注意力机制

一、本文介绍 作为入门性篇章&#xff0c;这里介绍了SCAM注意力在YOLOv8中的使用。包含SCAM原理分析&#xff0c;SCAM的代码、SCAM的使用方法、以及添加以后的yaml文件及运行记录。 二、SCAM原理分析 SCAM官方论文地址&#xff1a;SCAM文章 SCAM官方代码地址&#xff1a;SC…

Google Protocol Buffers快速入门指南

声明&#xff1a;未经作者允许&#xff0c;禁止转载。 概念 Portocol Buffer是谷歌提出来的一种序列化结构数据的机制&#xff0c;它的可扩展性特别强&#xff0c;支持C、C#、Java、Go和Python等主流编程语言。使用Portocol Buffer时&#xff0c;仅需要定义好数据的结构化方式…

Mysql梳理10——使用SQL99实现7中JOIN操作

10 使用SQL99实现7中JOIN操作 10.1 使用SQL99实现7中JOIN操作 本案例的数据库文件分享&#xff1a; 通过百度网盘分享的文件&#xff1a;atguigudb.sql 链接&#xff1a;https://pan.baidu.com/s/1iEAJIl0ne3Y07kHd8diMag?pwd2233 提取码&#xff1a;2233 # 正中图 SEL…

每日OJ题_牛客_添加逗号_模拟_C++_Java

目录 牛客_添加逗号_模拟 题目解析 C代码1 C代码2 Java代码 牛客_添加逗号_模拟 添加逗号_牛客题霸_牛客网 题目解析 读取输入&#xff1a;读取一行字符串。分割字符串&#xff1a;使用空格将字符串分割成单词数组。拼接字符串&#xff1a;将单词数组中的每个单词用逗号…

Oracle控制文件全部丢失如何使用RMAN智能恢复?

1.手动删除所有控制文件模拟故障产生 2.此时启动数据库发现控制文件丢失 3.登录rman 4.列出故障 list failure; 5.让RMAN列举恢复建议 advise failure; 6.使用RMAN智能修复 repair failure;

Java常用三类定时器快速入手指南

文章目录 Java常用三类定时器快速入手指南一、序言二&#xff0c;Timer相关1、概念2、Timer类3、TimerTask类4、ScheduleExecutorService接口 三&#xff0c;Scheduled相关1、配置1.1 SpringMVC配置1.2 SpringBoot配置&#xff08;1&#xff09;单线程&#xff08;2&#xff09…

cpp,git,unity学习

c#中的? 1. 空值类型&#xff08;Nullable Types&#xff09; ? 可以用于值类型&#xff08;例如 int、bool 等&#xff09;&#xff0c;使它们可以接受 null。通常&#xff0c;值类型不能为 null&#xff0c;但是通过 ? 可以表示它们是可空的。 int? number null; // …

桥接(桥梁)模式

简介 桥接模式&#xff08;Bridge Pattern&#xff09;又叫作桥梁模式、接口&#xff08;Interface&#xff09;模式或柄体&#xff08;Handle and Body&#xff09;模式&#xff0c;指将抽象部分与具体实现部分分离&#xff0c;使它们都可以独立地变化&#xff0c;属于结构型…

【AAOS】CarService -- Android汽车服务

概述 Android Automative OS理解为Android OS + Android Automative Service,而CarService就是提供汽车相关功能的最主要模块。 CarService与Android OS的关系:CarService运行于独立的进程中,其作为原有Android服务的补充,在汽车设备上运行。CarService在整体车载通信中起…

JAVA线程基础二——锁的概述之乐观锁与悲观锁

乐观锁与悲观锁 乐观锁和悲观锁是在数据库中引入的名词,但是在并发包锁里面也引入了类似的思想&#xff0c;所以这里还是有必要讲解下。 悲观锁指对数据被外界修改持保守态度&#xff0c;认为数据很容易就会被其他线程修改&#xff0c;所以在数据被处理前先对数据进行加锁&…

python 如何引用变量

在字符串中引入变量有三种方法&#xff1a; 1、 连字符 name zhangsan print(my name is name) 结果为 my name is zhangsan 2、% 字符 name zhangsan age 25 price 4500.225 print(my name is %s%(name)) print(i am %d%(age) years old) print(my price is %f%(pric…