信号处理快速傅里叶变换(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.mzph.cn/diannao/55415.shtml

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

相关文章

webpack信息泄露

先看看webpack中文网给出的解释 webpack 是一个模块打包器。它的主要目标是将 JavaScript 文件打包在一起,打包后的文件用于在浏览器中使用,但它也能够胜任转换、打包或包裹任何资源。 如果未正确配置&#xff0c;会生成一个.map文件&#xff0c;它包含了原始JavaScript代码的映…

课设实验-数据结构-线性表-手机销售

题目&#xff1a; 代码&#xff1a; #include<stdio.h> #include<string.h> #define MaxSize 10 //定义顺序表最大长度 //定义手机结构体类型 typedef struct {char PMod[10];//手机型号int PPri;//价格int PNum;//库存量 }PhoType; //手机类型 //记录手机的顺序…

【HTTP(3)】(状态码,https)

【认识状态码】 状态码最重要的目的&#xff0c;就是反馈给浏览器:这次请求是否成功&#xff0c;若失败&#xff0c;则出现失败原因 常见状态码: 200:OK&#xff0c;表示成功 404:Not Found&#xff0c;浏览器访问的资源在服务器上没有找到 403:Forbidden&#xff0c;访问被…

springboot系列--web相关知识探索三

一、前言 web相关知识探索二中研究了请求是如何映射到具体接口&#xff08;方法&#xff09;中的&#xff0c;本次文章主要研究请求中所带的参数是如何映射到接口参数中的&#xff0c;也即请求参数如何与接口参数绑定。主要有四种、分别是注解方式、Servlet API方式、复杂参数、…

【案例】距离限制模型透明

开发平台&#xff1a;Unity 2023 开发工具&#xff1a;Unity ShaderGraph   一、效果展示 二、路线图 三、案例分析 核心思路&#xff1a;计算算式&#xff1a;透明值 实际距离 / 最大距离 &#xff08;实际距离 ≤ 最大距离&#xff09;   3.1 说明 | 改变 Alpha 值 在 …

stm32f103调试,程序与定时器同步设置

在调试定时器相关代码时&#xff0c;注意到定时器的中断位总是置1&#xff0c;怀疑代码有问题&#xff0c;经过增大定时器的中断时间&#xff0c;发现定时器与代码调试并不同步&#xff0c;这一点对于调试涉及定时器的代码是非常不利的&#xff0c;这里给出keil调试stm32使定时…

自用Proteus(8.15)常用元器件图示和功能介绍(持续更新...)

文章目录 一、 前言二、新建工程&#xff08;以51单片机流水灯为例&#xff09;2.1 打开软件2.2 建立新工程2.3 创建原理图2.4 不创建PCB布版设计2.5 创建成功2.6 添加元器件2.7 原理图放置完成2.8 编写程序&#xff0c;进行仿真2.9 仿真 三、常用元器件图示和功能介绍3.1 元件…

【回眸】Tessy 单元测试软件使用指南(四)常见报错及解决方案与批量初始化的经验

前言 分析时Tessy的报错 1.fatal error: Tricore/Compilers/Compilers.h: No such file or directory 2.error: #error "Compiler unsupported" 3.warning: invalid suffix on literal;C11 requires a space between literal and string macro 4.error: unknown…

螺蛳壳里做道场:老破机搭建的私人数据中心---Centos下Docker学习01(环境准备)

1 准备工作 由于创建数据中心需要安装很多服务器&#xff0c;这些服务器要耗费很所物理物理计算资源、存储资源、网络资源和软件资源&#xff0c;作为穷学生只有几百块的n手笔记本&#xff0c;不可能买十几台服务器来搭建数据中心&#xff0c;也不愿意跑实验室&#xff0c;想躺…

文件上传之%00截断(00截断)以及pikachu靶场

pikachu的文件上传和upload-lab的文件上传 目录 mime type类型 getimagesize 第12关%00截断&#xff0c; 第13关0x00截断 差不多了&#xff0c;今天先学文件上传白名单&#xff0c;在网上看了资料&#xff0c;差不多看懂了&#xff0c;但是还有几个地方需要实验一下&#…

SpringBoot整合异步任务执行

同步任务&#xff1a; 同步任务是在单线程中按顺序执行&#xff0c;每次只有一个任务在执行&#xff0c;不会引发线程安全和数据一致性等 并发问题 同步任务需要等待任务执行完成后才能执行下一个任务&#xff0c;无法同时处理多个任务&#xff0c;响应慢&#xff0c;影响…

VirtualBox+Vagrant快速搭建Centos7系统【最新详细教程】

VirtualBoxVagrant快速搭建Centos7系统 &#x1f4d6;1.安装VirtualBox✅下载VirtualBox✅安装 &#x1f4d6;2.安装Vagrant✅下载Vagrant✅安装 &#x1f4d6;3.搭建Centos7系✅初始化Vagrantfile文件生成✅启动Vagrantfile文件✅解决 vagrant up下载太慢的问题✅配置网络ip地…

咸鱼sign逆向分析与爬虫实现

目标&#xff1a;&#x1f41f;的搜索商品接口 这个站异步有点多&#xff0c;好在代码没什么混淆。加密的sign值我们可以通过搜索找到位置 sign值通过k赋值&#xff0c;k则是字符串拼接后传入i函数加密 除了开头的aff…&#xff0c;后面的都是明文没什么好说的&#xff0c;我…

SysML案例-电磁轨道炮

DDD领域驱动设计批评文集>> 《软件方法》强化自测题集>> 《软件方法》各章合集>> 图片示例摘自intercax.com&#xff0c;作者是Intercax公司总裁Dirk Zwemer博士。

C题(六) 1到 100 的所有整数中出现多少个数字9

场景&#xff1a;编写程序数一下 1到 100 的所有整数中出现多少个数字9 控制循环的变量不可以随意改动&#xff01;&#xff01;&#xff01; 控制循环的变量不可以随意改动&#xff01;&#xff01;&#xff01; 控制循环的变量不可以随意改动&#xff01;&#xff01;&#x…

看480p、720p、1080p、2k、4k、视频一般需要多大带宽呢?

看视频都喜欢看高清&#xff0c;那么一般来说看电影不卡顿需要多大带宽呢&#xff1f; 以4K为例&#xff0c;这里引用一位网友的回答&#xff1a;“视频分辨率4092*2160&#xff0c;每个像素用红蓝绿三个256色(8bit)的数据表示&#xff0c;视频帧数为60fps&#xff0c;那么一秒…

数据结构--二叉树的顺序实现(堆实现)

引言 在计算机科学中&#xff0c;二叉树是一种重要的数据结构&#xff0c;广泛应用于各种算法和程序设计中。本文将探讨二叉树的顺序实现&#xff0c;特别是堆的实现方式。 一、树 1.1树的概念与结构 树是⼀种⾮线性的数据结构&#xff0c;它是由 n(n>0) 个有限结点组成…

C#串口温度读取

背景&#xff1a;每天学点&#xff0c;坚持 要安装好虚拟串口和modbus poll&#xff0c;方便调试&#xff08;相关资源在文末&#xff0c;也可以私信找我要&#xff09; 传感器部分使用的是达林科技的DL11B-MC-D1&#xff0c;当时42软妹币买的&#xff08;官网上面有这个传感…

若依--文件上传前端

前端 ry的前端文件上传单独写了一个FileUpload.Vue文件。在main.js中进行了全局的注册&#xff0c;可以在页面中直接使用文件上传的组件。全局导入 在main.js中 import 组件名称 from /components/FileUpLoadapp.compoent(组件名称) //全局挂载组件在项目中使用 组件命令 中…

828华为云征文|华为云 Flexus X 实例之家庭娱乐中心搭建

话接上文《828华为云征文&#xff5c;华为云Flexus X实例初体验》&#xff0c;这次我们利用手头的 Flexus X 实例来搭建家庭影音中心和密码管理环境。 前置环境 为了方便小白用户甚至运维人员&#xff0c;我觉得现阶段的宝塔面板 和 1Panel 都是不错的选择。我这里以宝塔为例…