Unity3D C# 信号FFT分析之5MHz超声波信号处理
目录
- 1 快速傅里叶分析
- 1.1 快速傅里叶(FFT)分析是什么
- 1.2 示波器采集到的实际波形如何进行傅里叶分析
- 1.3 FFT补零的作用
- 1.4 分贝dB
- 2 代码实现
- 2.1 Unity中使用MathNet.Numerics库
- 2.2 补零的FFT实现
- 2.3 将幅值换算为分贝值
- 2.4 提取需重点观察的频谱部分
- 2.5 思考问题
- 3 完整项目
- 4 参考文章
上一篇《Unity3D C# 信号FFT分析之自定义波形UI控件》我们为了将采集到的波形展示出来,特意写了个脚本来绘制波形图。
这一篇我们将用一个实际采集到的5MHz的超声波探头信号来进行fft分析,并将频谱图绘制出来。效果如下。
闲话少说,直接进入主题。咱们先来看看理论基础,然后再看代码怎么写。
1 快速傅里叶分析
1.1 快速傅里叶(FFT)分析是什么
学过信号与系统的同学应该知道,任何一个周期信号都可以分解为无数个正弦波的叠加,傅里叶分析就是把这些正弦波给求出来。
那么快速傅里叶又是什么呢?我们现实生活中实际采集到信号都是离散的,而且大多数不是周期信号,那能否对这样的信号进行傅里叶分析呢,当然可以的,对应的分析就是离散傅里叶分析。快速傅里叶就是离散傅里叶分析的一种实现。
其实博主之前一直有个问题没弄明白,我们为什么要进行傅里叶分析?
答案是换个角度看待问题,由于看待问题的角度变了,将会得出一些有意义的结论。一个杂乱无章的信号,比如说声音,从时域上看都是杂乱无章的,但是从频域上看就会发现男生的声音低频的能量较强,女生的声音高频的能量较强。
可参考李永乐老师的视频:傅立叶变换如何理解?美颜和变声都是什么原理?李永乐老师告诉你
1.2 示波器采集到的实际波形如何进行傅里叶分析
首先拿到一个波形的数据,我们必须先知道这个波形一共有多少个点,整个波形的采集时间是多长。比如Demo中提供的数据是2500个点,整个波形的时间是10μs。那么采集周期(即两个点的时间间隔)为10μs/2500 = 4ns,整个频谱带宽是1/4ns = 250MHz。因为我们的波形有2500个点,那么经过fft分析后也将得到一个2500个点的频谱数据,且频谱数据中任意相邻两个点的频率间隔是相同的,即Δf=1/10μs=0.1MHz(这也说明Demo中提供的波形,能分解出的正弦波最小正弦分量为0.1MHz)。
为什么是这样的对应关系以及傅里叶变换的完整解析,见《FFT后的物理意义》说得特别清楚。这里我就直接引用过来了。
我们以一个实际的信号为例来说明:
示波器采样得到的数字信号,就可以做FFT变换了。N个采样点,经过FFT之后,就可以得到N个点的FFT结果。为了方便进行FFT运算,通常N取2的整数次方。
假设采样频率为Fs,信号频率F,采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值,就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢?假设原始信号的峰值为A,那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。而第一个点就是直流分量,它的模值就是直流分量的N倍。而每个点的相位呢,就是在该频率下的信号的相位。第一个点表示直流分量(即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。如果要提高频率分辨率,则必须增加采样点数,也即采样时间。频率分辨率和采样时间是倒数关系。
下面这幅图更能够清晰地表示这种对应关系:
变换之后的频谱的宽度(Frequency Span)与原始信号也存在一定的对应关系。根据Nyquist采样定理,FFT之后的频谱宽度(Frequency Span)最大只能是原始信号采样率的1/2,如果原始信号采样率是4GS/s,那么FFT之后的频宽最多只能是2GHz。时域信号采样周期(Sample Period)的倒数,即采样率(Sample Rate)乘上一个固定的系数即是变换之后频谱的宽度,即 Frequency Span = K(1/ΔT),其中ΔT为采样周期,K值取决于我们在进行FFT之前是否对原始信号进行降采样(抽点),因为这样可以降低FFT的运算量。如下图所示:
可见,更高的频谱分辨率要求有更长的采样时间,更宽的频谱分布需要提高对于原始信号的采样率,当然我们希望频谱更宽,分辨率更精确,那么示波器的长存储就是必要的!它能提供您在高采样率下采集更长时间信号的能力!
这一点很重要,是我们的理论基础,只有弄明白了代码才写得出来。
1.3 FFT补零的作用
补零的作用有2个:
①在做fft分析时,如果点数是2的整数次方,这样可以提高运算效率(为什么可以提高,得看fft算法的代码具体实现了,这里我也没去详细了解,就不误人子弟啦)。所以如果如果数据点数不是以2为基数的整数次方,可以在原始数据开头或末尾补零,即将数据补到以2为基数的整数次方。由于Demo中的波形数据点数是2500点,不是2的整数次方所以需要补零,demo中是前后各补了一部分,一共补到了16384=214个点。
②补零后可以让频谱图看起来更加的平滑。注意,仅仅是看起来更加光滑,并不能降低可提取的最小分辨率。
这部分可看这篇文章《快速傅里叶变换(FFT)中为什么要“补零”》,也说得特别清楚。
1.4 分贝dB
大家在说声音大小的时候一般会用到分贝这个单位,那么分贝到底是什么呢?
分贝dB定义为两个数值的对数比率,这两个数值分别是测量值和参考值(也称为基准值)。
两种定义情况。
一种为功率之比:
一种为幅值之比:
下标为0的数值均为幅值和功率的参考值。功率量的例子如:声功率(W),声强(W/m2),电功率,电强等。幅值量的例子如:声压(Pa),电压(V),加速度(m/s2),温度等。但有一点要注意对于场量的幅值应该是RMS值,如声压场。
因为分贝值完全依赖于测量值与参考值之比,因此,计算时选择合适的参考值尤为关键。当测量结果相互比较时,这一点非常重要,选择的参考值不同,计算结果肯定不一样。常见信号的dB参考值如下表所示
注:没有特殊要求时,参考值通常为1。
具体见这篇文章《什么是分贝dB?》。
有几个dB值比较常用,需要重点关注。
6dB表示幅值的2倍;12dB表示幅值的4倍;20dB表示幅值的10倍;40dB表示幅值的100倍。
2 代码实现
2.1 Unity中使用MathNet.Numerics库
打开vs,进入Nuget包管理器。
然后输入Math进行搜索,并选中MathNet.Numerics,右侧选中点击安装。
将弹出提示,点击确定即可。
然后打开项目目录,找到Packages文件夹中的MathNet.Numerics插件文件夹。
里面有好几个版本,我们选择net461那个版本。
将其拷贝出来,然后在Unity的Assets文件夹下新建一个Plugins文件夹,将拷贝的dll复制过去。
然后我们就可以愉快的使用MathNet.Numerics进行fft分析啦。
2.2 补零的FFT实现
demo中是将信号前后补充相同个数的零到16384=214个点,至于为什么补零,见1.3节。
我们需要使用原有信号来构建一个Complex32(复数)的数组,然后才能进行fft分析。
补零的地方,实部和虚部全为0;采集到的信号部分,实部就用采集到的值,虚部为0。
复数数组构建完毕后使用Fourier.Forward来进行fft分析。
//fouierArr为刚构建的复数数组
Fourier.Forward(fouierArr, FourierOptions.Matlab);
然后将经过fft分析后的复数数组求模即可得到频率点的幅值数组。
public const int FftWithZeroAllPoints = 16384; // 补零后的波形总点数/// <summary>
/// 补零后的傅里叶变换.
/// </summary>
public static float[] FftAddZeros(float[] wave)
{int fftPoint = FftWithZeroAllPoints;// 补零,补零后的波形总点数为16384,分别在波形前面、后面各补一半// 补零是为了让频谱图看起来更加光滑,见这篇文章:https://blog.csdn.net/xingsongyu/article/details/103390317?depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-2&utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromBaidu-2int zeroPoint = (fftPoint - wave.Length) / 2;float[] real = new float[fftPoint]; // 实部float[] img = new float[fftPoint]; // 虚部for (int i = 0; i < fftPoint; i++){if (i >= zeroPoint && i < (zeroPoint + wave.Length)){real[i] = wave[i - zeroPoint];}else{real[i] = 0;}img[i] = 0;}Complex32[] fouierArr = new Complex32[fftPoint];for (int i = 0; i < fftPoint; i++){fouierArr[i] = new Complex32(real[i], img[i]);}Fourier.Forward(fouierArr, FourierOptions.Matlab);float[] result = new float[fftPoint]; // 傅里叶变换后的内容for (int i = 0; i < fftPoint; i++){result[i] = fouierArr[i].Magnitude;}return result;
}
2.3 将幅值换算为分贝值
幅值换算为分贝值的定义如下,见1.4节。
参考值我们选择为1,注意前面乘的是20。
// 计算分贝值
// 幅值转换为dB
// 参考值选择的1
float dbValue = (float)(20 * Math.Log(xx) / Math.Log(10));
2.4 提取需重点观察的频谱部分
上图是我们补零的信号后的完整频谱图。
大家会发现信号经过fft分析后是左右完全对称的,为什么是这样的?见这篇文章《FFT(快速傅里叶变换)中频率和实际频率的关系》它能给你答案。
我们只需关注红框部分,那才是我们真正的信号,其他部分都是补零的部分。
Demo中的给出的波形是0.5MHz~5MHz,那么我们只需要关注这部分的频谱图即可。那么问题来了,0.5MHz和5MHz分别对应到频谱数据中的哪个点呢?
1.2节中原理已经说得比较明白,结合代码相信大家一看就能明白。
//此波形的采集总时间是10微秒,共2500个点,则采集周期为4ns = 10μs/2500
//所以整个波形的带宽为250MHz = 2,5000,0000Hz = 1.0 / 4ns
//这个公式是怎么来的?见这篇文章 FFT后的物理意义:https://blog.csdn.net/iloveyoumj/article/details/53308142
float allBandWidth = 1f / (4 * Mathf.Pow(10, -9));
//我们只观察0.5MHz~5MHz这段数据
//为什么只观察这段数据?因为用仪器(这里用的阻抗仪)来产生这段波形的时候,使用的频率范围是0.5MHz~5MHz
//这个范围不影响分析结果的,只是提取出来方便我看一点
//0.5MHz和5MHz对应频谱上的点
int startFrequencyPointOnFFT = (int)(FftUtil.FftWithZeroAllPoints * 5 * Mathf.Pow(10, 5) / allBandWidth); //0.5MHz
int stopFrequencyPointOnFFT = (int)(FftUtil.FftWithZeroAllPoints * 5 * Mathf.Pow(10, 6) / allBandWidth); //5MHz
//对0~0.5MHz这部分进行高通滤波(其实就是缩小点倍数)
for (int i = 0; i < startFrequencyPointOnFFT; i++)
{fftAddZeros[i] *= 0.05f;
}
float[] watchRangeFftData = new float[stopFrequencyPointOnFFT];
Array.Copy(fftAddZeros, watchRangeFftData, stopFrequencyPointOnFFT);
代码中对0.5MHz之前的部分进行了高通滤波,其实就是将之前的部分乘了个很小的数,如0.05,将其过滤掉。
Demo中还画出了一个只观察[最高dB-44dB,最高dB]的频谱图,为什么提取这部分出来呢,因为这部分能量最强,我们需要重点观测。三个频谱图如下。
2.5 思考问题
大家有空可以想想怎么求出最高dB-6dB的起始频率、结束频率、中心频率和相对带宽。
最高dB-6dB起始频率和结束频率很好求,最高dB我们是能计算出来的,那么-6dB的值我们也知道了。然后根据这个值我们去0.5MHz~5MHz这部分的频谱数据中去找到最接近的那两个点(其实就是遍历就行了)。
-6dB中心频率 = 结束频率 - 起始频率 / 2
-ddB相对带宽 = (结束频率 - 起始频率) / 中心频率 * 100
最后计算出的结果如下:
-6dB起始频率 1.98MHz
-6dB结束频率 2.98MHz
-6dB中心频率 2.48MHz
-6dB相对带宽 40.56%
3 完整项目
项目上传到这儿了。
链接:https://pan.baidu.com/s/1qrkavvpUnw6Nkfh-2NK-FQ
提取码:k8ex
博主本文博客链接。
4 参考文章
- FFT后的物理意义
- 快速傅里叶变换(FFT)中为什么要“补零”
- 什么是分贝dB?
- FFT(快速傅里叶变换)中频率和实际频率的关系
- 傅立叶变换如何理解?美颜和变声都是什么原理?李永乐老师告诉你
Unity3D C# 信号FFT分析之5MHz超声波信号处理相关推荐
- 信号处理与分析-确定性信号的分析
目录 一.引言 二.确定性信号的定义 三.确定性信号的分类 四.确定性信号的分析方法 4.1 傅里叶变换 4.2 离散傅里叶变换 4.3 离散余弦变换 4.4 小波变换 五.确定性信号的处理方法 六. ...
- 在FFT分析在而立之年的展望与总结
这是Sri Welaratna 在 原文 : https://www.dataphysics.com/30_Years_of_FFT_Analyzers_by_Sri_Welaratna.pdf# ...
- FFT分析的加窗和重叠
FFT分析的加窗和重叠 这将通过实例说明加窗和重叠对频谱分析的影响.用10Hz的正弦波,以说明重叠窗函数在频谱分析过程中的不同之处. 重叠分析,就是连续分析的时域数据块通过指定的时间纪录百分比进行重叠 ...
- I/Q信号解调分析过程
现代的卫星通信.5G移动通信.无线局域网的通信带宽都非常宽,可能超过500MHz甚至达到2GHz左右.对于这么宽带的信号来说,传统频谱仪虽然可以看到信号频谱,但受限于实时分析带宽,已经很难再对信号的调 ...
- 飞控中加速度计数据fft分析
加速度计原始数据高频噪声很严重,使用时需要设计合适的滤波器将其滤掉. matlab对加速度计原始数据,30Hz低通滤波之后的数据进行fft分析如下: clc;source=simout1; filte ...
- 语音信号线性预测分析(MATLAB实战篇)
文章目录 前言 基本概念 基本参数的求解及其用途 1.线性预测系数(LPC) 2.线性预测系数LPC的频谱 3.线性预测系数的倒谱LPCC 4.线性预测误差e(n)及其自相关 5.预测误差滤波器A(z ...
- 数字示波器FFT分析
数字示波器的FFT分析功能 FFT是一项很强大的分析功能,在数字示波器中普遍存在,基于先进的FFT分析,设计人员可以准确了解信号中引入的干扰信号频点,信号功率谱,信号频率构成,滤波电路截频特性等.为了 ...
- matlab时域信号如何分析方法,信号时域采样频谱分析(matlab)
<信号时域采样频谱分析(matlab)>由会员分享,可在线阅读,更多相关<信号时域采样频谱分析(matlab)(12页珍藏版)>请在读根文库上搜索. 1.基于matlab的时域 ...
- 信号完整性分析2——时域与频域
信号完整性分析2--时域与频域 Date:2020/06/07 2.1 时域 定义: 时域就是真实世界,是唯一实际存在的域 Fclock=1TclockF_{clock}={\frac {1}{T_{ ...
最新文章
- SQL脚本--有关压缩数据库日志
- 【java读书笔记】——java的异常处理
- linux vlc流媒体服务器,vlc media server rtsp 流媒体服务器搭建成功经验分享
- CI中创建你自己的类库
- Webservice更新时出错。下载”。。。”时出错。请求失败,错误信息为:
- Qt中调用C语言函数库
- 计算机考研专业课——c语言
- 【转】推荐系统入门实践:世纪佳缘会员推荐(完整版)
- Smart-Link、Monitor-Link介绍与配置举例
- 微软软件实现技术授课系列内容之五:软件测试基础
- 【JAVA基础】java基础之-泛型详解
- 图片标签z-index设置不起作用
- 有利可图网_有利可图的项目手册-现在可用
- python作排产计划表_生产排程计划表
- springboot easyexcel不创建对象导入excel 通用版
- java基本类型char
- Oracle与SQL Server在企业应用中的比较(转)
- file上传代码 ios_iOS-实现文件上传下载
- 撮合系统的数据流转过程
- 个人知识管理系统 mysql_个人知识管理系统Version1.0开发记录(12)