Unity3D C# 信号FFT分析之5MHz超声波信号处理

Unity3D C# 信号FFT分析之5MHz超声波信号处理

2020年9月13日 0 作者 老王

上一篇《Unity3D C# 信号FFT分析之自定义波形UI控件》我们为了将采集到的波形展示出来,特意写了个脚本来绘制波形图。
这一篇我们将用一个实际采集到的5MHz的超声波探头信号来进行fft分析,并将频谱图绘制出来。效果如下。
5MHz探头频谱图
闲话少说,直接进入主题。咱们先来看看理论基础,然后再看代码怎么写。

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。如果要提高频率分辨率,则必须增加采样点数,也即采样时间。频率分辨率和采样时间是倒数关系。
下面这幅图更能够清晰地表示这种对应关系:
T与Δf的关系
变换之后的频谱的宽度(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的运算量。如下图所示:
03.ΔT与完整带宽的关系
可见,更高的频谱分辨率要求有更长的采样时间,更宽的频谱分布需要提高对于原始信号的采样率,当然我们希望频谱更宽,分辨率更精确,那么示波器的长存储就是必要的!它能提供您在高采样率下采集更长时间信号的能力!

这一点很重要,是我们的理论基础,只有弄明白了代码才写得出来。

1.3 FFT补零的作用

补零的作用有2个:
①在做fft分析时,如果点数是2的整数次方,这样可以提高运算效率(为什么可以提高,得看fft算法的代码具体实现了,这里我也没去详细了解,就不误人子弟啦)。所以如果如果数据点数不是以2为基数的整数次方,可以在原始数据开头或末尾补零,即将数据补到以2为基数的整数次方。由于Demo中的波形数据点数是2500点,不是2的整数次方所以需要补零,demo中是前后各补了一部分,一共补到了16384=214个点。
②补零后可以让频谱图看起来更加的平滑。注意,仅仅是看起来更加光滑,并不能降低可提取的最小分辨率。
这部分可看这篇文章《快速傅里叶变换(FFT)中为什么要“补零”》,也说得特别清楚。

1.4 分贝dB

大家在说声音大小的时候一般会用到分贝这个单位,那么分贝到底是什么呢?

分贝dB定义为两个数值的对数比率,这两个数值分别是测量值和参考值(也称为基准值)。
两种定义情况。
一种为功率之比:
dB定义功率执笔
一种为幅值之比:
幅值之比
下标为0的数值均为幅值和功率的参考值。功率量的例子如:声功率(W),声强(W/m^2),电功率,电强等。幅值量的例子如:声压(Pa),电压(V),加速度(m/s^2),温度等。但有一点要注意对于场量的幅值应该是RMS值,如声压场。
因为分贝值完全依赖于测量值与参考值之比,因此,计算时选择合适的参考值尤为关键。当测量结果相互比较时,这一点非常重要,选择的参考值不同,计算结果肯定不一样。常见信号的dB参考值如下表所示
常见信号的dB参考值
注:没有特殊要求时,参考值通常为1。

具体见这篇文章《什么是分贝dB?》。
有几个dB值比较常用,需要重点关注。
6dB表示幅值的2倍;12dB表示幅值的4倍;20dB表示幅值的10倍;40dB表示幅值的100倍。
分贝与幅值比对应关系

2 代码实现

2.1 Unity中使用MathNet.Numerics库

打开vs,进入Nuget包管理器。
打开NuGet包管理器
然后输入Math进行搜索,并选中MathNet.Numerics,右侧选中点击安装。
选中安装
将弹出提示,点击确定即可。
NuGet安装提示
然后打开项目目录,找到Packages文件夹中的MathNet.Numerics插件文件夹。
库文件夹
里面有好几个版本,我们选择net461那个版本。
MathNet.Numerics.dll net461版本
将其拷贝出来,然后在Unity的Assets文件夹下新建一个Plugins文件夹,将拷贝的dll复制过去。
Plugins文件夹下的MathNet.Numerics插件
然后我们就可以愉快的使用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-2
    int 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分析后是左右完全对称的,为什么是这样的?见这篇文章《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 参考文章