65.9K
CodeProject 正在变化。 阅读更多。
Home

DSPLib - 适用于 .NET 4 的 FFT / DFT 傅里叶变换库

starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.93/5 (68投票s)

2016年6月19日

MIT

26分钟阅读

viewsIcon

209031

downloadIcon

18030

DSPLib 是一个完整的 DSP 库,是使用 .NET 4 执行 FFT 的端到端解决方案

DSPLib Splash

引言

市场确实需要一个即用型傅里叶变换库,用户可以直接使用它来执行快速傅里叶变换 (FFT) 或离散傅里叶变换 (DFT),并获得经典的频谱与频率图。

您在商业软件包、开源库、教科书和网络上找到的绝大多数代码都根本不适合这项任务,需要数小时的进一步调整才能获得经典且正确缩放的频谱图。

此处介绍的库是一个实用、有组织且完整的 .NET 4+ 开源 DSP 相关例程库,根据非常宽松的 MIT 许可发布。

DSPLib 的作用

DSPLib 有几个主要部分,但其基本目标是允许对时间序列输入数组执行实际的傅里叶变换,从而生成可用的经典频谱输出,而无需用户进行任何进一步的调整。

基本的傅里叶变换 (FT) 有两种基本类型:最一般的形式可以从任意长度的输入数据生成频谱输出。这种类型的变换称为离散傅里叶变换 (DFT)。代码简单粗暴。

  • 优点是:输入数据可以是任意长度。
  • 缺点是:由于它是一种通用方法,计算量大,大型输入数据集的计算时间可能非常长。

一种更具体的 FT 类型称为快速傅里叶变换 (FFT)。

  • 优点是:它的计算速度比 DFT 快得多。
  • 缺点是:输入数据长度限制为2的幂次。代码理解起来更复杂。

例如:我的 i7 电脑上,一个 8192 点的 FFT 耗时不到 1 毫秒。相同长度的 DFT 耗时 360 毫秒。因此,您可以看到为什么 FFT 在实时应用中比暴力 DFT 更受欢迎。

DSPLib 实现了两种傅里叶变换。

注意:在本文的其余部分中,当讨论既适用于“FFT”又适用于“DFT”时(因为它们对等效输入数据产生相同的结果),广义傅里叶变换将被称为“FT”。

所有 FT 都可以接受实数或复数输入,并自然产生复数结果。迄今为止生成的所有 FT 都实现了它们自己的复数类型,这自然导致了库之间的不兼容性。.NET 4.0 最终(在 System.Numerics 命名空间中)包含一个复数结构和许多用于复数运算的数学方法。DSPLib 集成了 .NET 4 复数类型。

为了实现真正的速度改进并自动在具有多核和/或线程的处理器上缩放速度,DSPLib 还实现了 .NET 4 中内置的任务并行扩展。这导致执行时间真正改进,并自动随可用处理器核心/线程数量而缩放。例如,在 2 核/4 线程的 i7 处理器上,使用任务并行扩展将 10,000 点 DFT 的执行时间缩短了 3 倍以上。

对较小 DFT 的正弦和余弦乘法项进行智能缓存也使性能提高了约 3 倍。

即使在低端 i7 处理器上,这两个易于实现的功能也将原始 DFT 速度提高了约 9 倍。

频谱的实部和虚部

所有 FT 实现自然会产生一个与输入数组长度相同的输出数据数组。然而,输出不仅包含复数,还包含频谱本身的实部和虚部——有时根据人们如何看待它而称为负频率。频谱的虚部只是实部的镜像,不包含关于频谱的任何额外信息。

DSPLib 只返回频谱输出数据数组的实部(输入数据数组长度的一半,不包括 FS/2 数据点)。您会发现所有其他库都将整个实部和虚部作为返回值提供。这让许多人感到困惑,并且对于生成通常的真实频谱数据表示没有用处。

注意:如果出于某种数学原因需要复杂输出频谱的虚部和实部,那没问题,您知道为什么需要它。DSPLib 旨在生成经过缩放的经典频谱图,而不是一般的 FT 数学。

傅里叶变换缩放

您会发现大多数 FT 实现都没有针对经典频谱显示进行正确缩放。也就是说,变换信号的幅度随变换点的数量 N 而变化。这让大多数用户感到不安。大多数用户期望获得与 FT 中使用的点数无关的正确缩放输出。

例如,如果用户输入一个 1 Vrms 幅度的 10 Hz 信号,那么他们希望在频谱图中看到 10 Hz 处的 1 Vrms 峰值。DSPLib 是迄今为止唯一一个正确缩放 FT 输出的库,无论您如何更改变换中的点数,它们都是准确的。

注意:所有 DFT 和 FFT 实施在 DC (Bin 0) 和采样频率一半 (FS/2) 处也存在缩放差异。这些点没有虚部,只有实部。因此,这些点的所需缩放也不同。DSPLib 对 DC 点应用适当的比例。FS/2 点被认为是虚部的,并且按照惯例不包含在频谱输出的最终实部中。

添加数据窗函数后的缩放

您可能知道,所有 FT 都假设您输入的有限时间数据块实际上在数据块之前和之后都无限延伸。如果您使用偶数个周期变换纯正弦波,则频谱显示将是完美的峰值,位于输出频谱的一个 bin 中。

在现实生活中,通常不可能生成完全偶数的输入序列。任何小数周期都会导致输出在整个结果频谱中模糊不清。小数周期的幅度也会显得低于其实际值。

这就是使用窗函数的作用。窗函数应用于输入数据,它将时间序列的开始和结束逐渐减小到零或接近零。因此,使用窗函数可以欺骗 DFT,使其不会注意到输入数据中的不连续性。

窗函数通常根据四个标准进行选择

  1. 它们将不希望的旁瓣推向多低的幅度
  2. 它们减少幅度不确定性的程度
  3. 它们将数据模糊到多少个 FT 频点中
  4. 历史原因,即“在这种情况下,每个人都使用窗口:xyz 等。”

选择最佳窗口涉及许多权衡。

DSPLib 的窗函数基于马克斯·普朗克研究所工作人员撰写的一篇优秀文章 [1],包括所有常见的窗函数类型,例如:Welch、Hamming 和 Hann(或俗称 Hanning),以及其他 27 种非常有用但可能鲜为人知的窗函数。还包括 17 种非常独特的平顶窗函数,即使输入频率在 FT 输出频点之间,它们也能提供非常低的幅度误差。

由于对输入数据进行加窗会物理地改变它,因此它也可能改变最终的 FT 频谱输出幅度。DSPLib 包含两个函数,可以对窗增益(或损耗)进行补偿,以便可以找到并显示正确缩放的 FT 输出。其中一个比例因子与离散信号的峰值有关。当测量噪声时,必须使用另一种类型的比例因子。

生成一组窗函数系数后,将其输入到适当的比例因子函数中,并计算出一个标量比例因子。通过乘以结果频谱来应用此比例因子,以校正幅度。

有关窗函数的需求和选择的更多信息,请参阅参考资料 [2]。

零填充和缩放

零填充是 FFT 和 DFT 中一个非常有用的技巧。一种用途是将输入数据长度向上取整到下一个2的幂,以便可以使用更快的 FFT 进行变换。例如,如果您有 1000 个输入数据点,您可以对数据进行额外的 24 个点零填充,这样就可以执行 1024 点 FFT,而不是较慢的 1000 点 DFT(我的电脑上速度提高了近 9 倍)。

零填充的另一个用途是使显示看起来更像预期的频谱显示。零填充通过在显示屏上插入点来扩大任何峰值,这使得图表在显示屏或绘图上看起来更好。零填充还对计算点之间的空间进行插值,当信号不直接位于 FT 的 bin 中心时,可以减少幅度误差。

最后,零填充可以非常精细地观察窗函数添加的旁瓣和边带抑制。

即使信号位于 bin 中心,零填充也会影响最终的频谱输出幅度。此处介绍的 FT 例程接受一个包含所需零填充数量的初始化参数,以便可以自动考虑正确的重新缩放因子。由于 FT 例程了解出于缩放原因所需的零填充,因此它们会向输入数据添加请求的零数量,这样用户就不必自己添加了。

参考 [2] 中有更多关于零填充及其更多用途的信息。

测试数据生成

在程序开发和调试过程中,并不总是能获得真实数据进行分析。DSPLib 通过包含两个校准过的正弦波发生器和一个校准过的高斯随机噪声发生器来帮助生成包含噪声的逼真测试信号。

通用频谱输出操作方法

如前所述,任何 FT 的原始输出都是一个复数数组。除了 .NET 4 复数类型运算符之外,DSPLib 还包含许多易于使用的转换例程,用于将复数结果转换为线性或对数幅度格式的标量数组。

此外,还包含一些对标量数组进行操作的数学例程。这些例程允许在两个数组之间或数组与常数之间进行加、减、乘和除运算。

示例

关于 DSPLib 的一般方面就谈到这里。下面的示例将展示它在实践中应用起来是多么容易。

示例 1

void example1()
{
  // Generate a test signal,
  //  1 Vrms at 20,000 Hz
  //  Sampling Rate = 100,000 Hz
  //  DFT Length is 1000 Points
  double amplitude = 1.0;
  double frequency = 20000;
  UInt32 length = 1000;
  double samplingRate = 100000;

  double[] inputSignal = 
           DSP.Generate.ToneSampling(amplitude, frequency, samplingRate, length);

  // Instantiate a new DFT
  DFT dft = new DFT();

  // Initialize the DFT
  // You only need to do this once or if you change any of the DFT parameters.
  dft.Initialize(length);

  // Call the DFT and get the scaled spectrum back
  Complex[] cSpectrum = dft.Execute(inputSignal);

  // Convert the complex spectrum to magnitude
  double[] lmSpectrum = DSP.ConvertComplex.ToMagnitude(cSpectrum);

  // Note: At this point, lmSpectrum is a 501 byte array that 
  // contains a properly scaled Spectrum from 0 - 50,000 Hz (1/2 the Sampling Frequency)

  // For plotting on an XY Scatter plot, generate the X Axis frequency Span
  double[] freqSpan = dft.FrequencySpan(samplingRate);

  // At this point a XY Scatter plot can be generated from,
  // X axis => freqSpan
  // Y axis => lmSpectrum

  // In this example, the maximum value of 1 Vrms is located at bin 200 (20,000 Hz)
}

示例 1 – 使用 DFT 测量信号的基本设置

在上面的示例中,生成了一个测试信号,实例化并初始化了一个 DFT,然后执行了 DFT。对结果进行缩放,幅值结果在数组 lmSpectrum 中可用。

注意:要在您自己的项目中使用此代码,请参阅下面的“将 DSPLib 添加到您的项目”部分。

Plot of the results from Example 1

示例 1 – 示例 1 代码片段的输出图

那 FFT 呢?

使用 FFT 的过程与示例 1 中的 DFT 完全相同,只是:输入数据长度加上任何零填充必须是 2 的精确幂。

可以替换为实例化 FFT 的实际代码

DSPLib.FFT fft = new DSPLib.FFT();
fft.Initialize(inputDataArray.Length);
Complex[] cSpectrum = fft.Execute(inputDataArray);

示例 2

void example2()
{
  // Same Input Signal as Example 1 - Except a fractional cycle for frequency.
  double amplitude = 1.0; double frequency = 20000.5;
  UInt32 length = 1000; UInt32 zeroPadding = 9000; // NOTE: Zero Padding
  double samplingRate = 100000;

  double[] inputSignal = DSPLib.DSP.Generate.ToneSampling
                         (amplitude, frequency, samplingRate, length);

  // Apply window to the Input Data & calculate Scale Factor
  double[] wCoefs = DSP.Window.Coefficients(DSP.Window.Type.Hamming, length);
  double[] wInputData = DSP.Math.Multiply(inputSignal, wCoefs);
  double wScaleFactor = DSP.Window.ScaleFactor.Signal(wCoefs);

  // Instantiate & Initialize a new DFT
  DSPLib.DFT dft = new DSPLib.DFT();
  dft.Initialize(length, zeroPadding); // NOTE: Zero Padding

  // Call the DFT and get the scaled spectrum back
  Complex[] cSpectrum = dft.Execute(wInputData);

  // Convert the complex spectrum to note: Magnitude Format
  double[] lmSpectrum = DSPLib.DSP.ConvertComplex.ToMagnitude(cSpectrum);

  // Properly scale the spectrum for the added window
  lmSpectrum = DSP.Math.Multiply(lmSpectrum, wScaleFactor);

  // For plotting on an XY Scatter plot generate the X Axis frequency Span
  double[] freqSpan = dft.FrequencySpan(samplingRate);

  // At this point a XY Scatter plot can be generated from,
  // X axis => freqSpan
  // Y axis => lmSpectrum
}

示例 2 – 使用 DFT 和加窗测量信号的基本设置

在上面的示例中,生成了一个测试信号,但这次有一个小数周期 (20000.5 Hz),实例化并初始化了一个 DFT,并进行了零填充以显示精细的窗口细节。最后,执行 DFT。对结果进行缩放,幅值结果在数组 lmSpectrum 中可用。

Plot of the results from Example 2

示例 2 A – 示例 2 代码片段的输出图

零填充显示精细的频谱细节。

Plot of the results from Example 2

示例 2 B – 示例 2 代码片段的输出图

X 轴放大后显示了由于应用了窗函数和添加了零填充而产生的精细频谱细节。另请注意,即使在这种情况下 DFT 长度为 1000 点 + 9000 零填充 = 总共 10,000 点,幅度仍然正确。幅度校正在 DSPLib 中很容易处理:零填充在 FT 本身中处理,应用窗函数产生的幅度误差通过窗函数比例因子函数进行校正。

DSPLib 基本使用步骤回顾

DSPLib 的设计旨在保持 DFT 和 FFT 使用的一致性。此外,可以按照下面概述的相同基本步骤分析信号和噪声类信号。

基本使用步骤 - 信号幅度分析

  1. 实例化并初始化 FT,计算信号的窗函数比例因子
  2. 获取或生成输入数据
  3. 对输入数据应用窗函数
  4. 对加窗的输入数据执行 FT
  5. 转换为幅值格式
  6. 如果需要:逐个频点平均,然后对所有平均重复步骤 2-6
  7. 对数据应用窗函数比例因子

结果:信号的正确缩放频谱

基本使用步骤 - 噪声信号分析

  1. 实例化并初始化 FT,计算噪声的窗函数比例因子
  2. 获取或生成输入数据
  3. 对输入数据应用窗函数
  4. 对加窗的输入数据执行 FT
  5. 转换为平方幅值格式
  6. 如果需要:逐个频点平均平方幅值数据,然后对所有平均重复步骤 2-6
  7. 转换为幅值格式
  8. 对数据应用窗函数比例因子。

结果:噪声信号的正确缩放频谱。

注意:信号和噪声必须区别对待的原因在参考资料 [2] 中有所阐述。

测试应用和示例说明

代码下载包含一个示例测试应用程序,其中包含几个测试函数代码片段,还包含一个基于 GUI 的应用程序,允许对用户可定义的测试信号执行 DFT 和 FFT。这对于理解 DSPLib 的使用非常有帮助。

使用 DSPLib 的示例代码片段包含在“DSPLib_Test”解决方案中的一个名为“DSPLib_ExampleCode.cs”的文件中,其中包括以下示例:

  • 示例 1 – 一个基本的 DFT 示例,它生成数据并产生频谱输出
  • 示例 2 – 与示例 1 相同,但还包括添加窗函数和比例因子
  • 示例 3 – 与示例 2 相同,但还包括添加零填充
  • 示例 4 – 与示例 2 相同,但最终结果为对数幅度 dBV 单位
  • 示例 5 – 带窗函数的 FFT 示例
  • 示例 6 – 带噪声信号和窗函数的 DFT,演示了平均和测量噪声的正确方法
  • 示例 7 - 带窗函数和零填充的 FFT
  • 示例 8 - 带信号和底噪的 DFT(模拟真实系统)
  • 示例 9 - DFT 相位提取测试
  • 示例 10 - 带窗函数、零填充和自动峰值检测的 FFT 相位提取测试
  • 示例 11 - DFT 显示任务并行 DFT 也可以在后台工作线程上运行
  • 示例 12 - 复数输入数据 (I&Q) 的 DFT,还演示了如何对 IQ 数据进行加窗以及如何应用比例因子

完整库文档

下面将介绍 DSPLib 的每个部分和功能的完整库描述。该库还实现了 XML 注释,可在 Visual Studio 环境中启用 IntelliSense 帮助。

傅里叶变换

void DSPLib.DFT()

void DSPLib.FFT()

实例化 DFT 或 FFT 类。

注意:这些变换被实例化为对象,因为在实时应用中,我经常发现同时运行多个具有不同长度的 FT 很方便。通过使用独立对象,无需浪费时间在执行变换之间重新定义 FT 长度参数。

void DSPLib.DFT.Initialize
  (UInt32 inputDataLength, UInt32 zeroPaddingLength = 0, bool forceNoCache = false)

void DSPLib.FFT.Initialize(UInt32 inputDataLength, UInt32 zeroPaddingLength = 0)
  • 在使用 DFT 或 FFT 类之前初始化其内部变量
  • 必须在实例化 DFT 或 FFT 对象后调用
  • 只需调用一次,或者当 FT 定义改变时
  • InputDataLength - 要进行变换的时间序列的长度
  • (Optional) zeroPaddingLength - 在变换之前要添加到输入数据数组的零的可选数量。FT“Execute”函数将向输入数据数组添加请求的零数量,这样用户就不必自己添加了。
  • DFT 特定参数和属性

DFT 尝试在内存中分配一些常见的正弦和余弦乘法项,以加快 DFT 计算。对于小型 DFT,这不会产生太大影响。对于非常大的 DFT,无论如何都没有足够的内存来分配数组,所以那时也不会产生任何影响。

DFT 初始化函数包含一个可选参数:forceNoCache。当设置为 true 时,此参数将阻止 DFT 在 DFT 计算中使用预分配内存。这可能会节省堆栈空间内存,但也可能导致执行时间更长。默认值为:forceNoCache = false

通常,如果您在程序中使用多个 DFT,您应该首先初始化最大的 DFT,以确保它们有足够的程序内存来分配数组以加快 DFT 执行。

DFT 对象实现了一个名为 IsUsingCached 的布尔只读属性。可以查询此属性以查看 DFT 在当前初始化时是否正在使用缓存值。

Complex[] DSPLib.DFT.Execute(double[] inputData)
Complex[] DSPLib.DFT.Execute(Complex[] inputData)

Complex[] DSPLib.FFT.Execute(double[] inputData)
Complex[] DSPLib.FFT.Execute(Complex[] inputData)
  • 在提供的 inputData 数组上执行 FT
  • 返回频谱结果的 System.Numerics.Complex[] 数组
  • 只返回频谱的实部(输入数据长度的一半)。这对应于:DC 到采样频率 / 2。
  • 注意:在 V2.0 中,DFT 和 FFT Execute 方法添加了一个接受 Complex[] 数据输入类型的重载。
double[] DSPLib.DFT.FrequencySpan(double samplingRateHz)

double[] DSPLib.FFT.FrequencySpan(double samplingRateHz)

一个辅助函数,用于返回与 FT 输出数组中每个值对应的频率点数组。此函数在需要绘制幅度与频率图的应用程序中非常有用。此函数依赖于首先使用正确参数调用“Initialized”函数。

此数组仅在 FT 中的点数或零填充长度发生变化,或采样率发生变化时才需要重新生成。如果您再次调用 .Initialize(),您可能需要再次调用此函数。

傅里叶变换格式转换

double[] DSPLib.ConvertComplex.ToMagnitude(Complex[] cSpectrum)

double[] DSPLib.ConvertComplex.ToMagnitudeDBV(Complex[] cSpectrum)

double[] DSPLib.ConvertComplex.ToMagnitudeSquared(Complex[] cSpectrum)

double[] DSPLib.ConvertComplex.ToPhaseDegrees(Complex[] cSpectrum)

double[] DSPLib.ConvertComplex.ToPhaseRadians(Complex[] cSpectrum)

这组函数将 FT 输出的 Complex[] 数组转换为更方便最终用户使用的格式。提供多种输出格式

  • ToMagnitude: 表示每个频点处复数的幅值
  • ToMagnitudeDBV: 表示以 1 Vrms 为参考的对数10格式幅值。例如:1 Vrms = 0 dBV,0.1 Vrms = -20 dBV
  • ToMagnitudeSquared: 表示频谱的幅值平方值。幅值平方 (V2) 与功率成正比。由于噪声与功率成正比,因此在平均噪声时,幅值平方是正确的格式。
  • ToPhaseDegrees: 表示每个点的复数量的相角(以度为单位)
  • ToPhaseRadians: 表示每个点的复数量的相角(以弧度为单位)
double[] DSPLib.ConvertMagnitude.ToMagnitudeSquared(double[] spectrum)

double[] DSPLib.ConvertMagnitude.ToMagnitudeDBV(double[] spectrum)

这组函数将 FT 转换后的频谱输出的幅度 double[] 数组转换为更方便最终用户使用的格式。提供多种输出格式

  • ToMagnitude: 表示每个频点处复数的幅值。
  • ToMagnitudeDBV: 表示以 1 Vrms 为参考的对数10格式幅值。例如:1 Vrms = 0 dBV,0.1 Vrms = -20 dBV
  • ToMagnitudeSquared: 表示频谱的幅值平方值。幅值平方 (V2) 与功率成正比。由于噪声与功率成正比,因此在平均噪声时,幅值平方是正确的格式。
double[] DSPLib.ConvertMagnitudeSquared.ToMagnitude(double[] spectrum)

double[] DSPLib.ConvertMagnitudeSquared.ToMagnitudeDBV(double[] spectrum)

将 FT 转换后的频谱输出的 MagnitudeSquared double[] 数组转换为更方便最终用户使用的格式的函数。提供多种输出格式

  • ToMagnitude: 表示每个点复数的幅值
  • ToMagnitudeDBV: 表示以 1 Vrms 为参考的对数10格式幅值。例如:1 Vrms = 0 dBV,0.1 Vrms = -20 dBV

频谱数据分析函数

double DSPLib.Analyze.FindMaxAmplitude(double[] inData)

double DSPLib.Analyze.FindMaxPosition(double[] inData)

double DSPLib.Analyze.FindMaxFrequency(double[] inData, double[] fSpan)

给定一个输入数组,这些函数找到指定的最大量。可以是数组中找到的最大幅度、找到最大幅度的 bin 位置或找到最大幅度的频率。这些函数设计用于频谱输出,并直接接受 ConvertTo 函数的任何格式输出。

double DSPLib.Analyze.FindMean(double[] inData, UInt32 startBin=10, UInt32 stopBin=10)

double DSPLib.Analyze.FindRms(double[] inData, UInt32 startBin=10, UInt32 stopBin=10)

给定一个输入数组,这些函数确定数组的平均值 (Average) 或 RMS 值。由于存在从 DC 或 Fs/2(最后一个数据点)到频谱的泄漏,可选参数允许从结果中排除第一个和最后一个指定的 bin。默认情况下,从两端排除 10 个 bin 足以排除 DSPLib 提供的任何窗口的所有泄漏。

double[] DSPLib.Analyze.UnwrapPhaseDegrees(double[] inPhaseDeg)

double[] DSPLib.Analyze.UnwrapPhaseRadians(double[] inPhaseRad)

计算出的相位总是在 +/-PI 区间处存在不连续性。这些函数尝试解卷输入相位数组以消除跳跃型不连续性。这对于干净信号效果很好,但随着信噪比的下降,很难从噪声中确定真正的不连续性。这些函数是一个很好的起点,但对于真实世界的信号和低信噪比,最佳性能是在使用足够信号平均时。其行为是根据等效的 Octave / Matlab 功能建模的。

窗函数系数生成

enum DSPLib.Window.Type.{WindowName}

double[] DSPLib.Window.Coefficients(DSPLib.Window.Type windowName, UInt32 points)

窗函数使用类型为“DSPLib.Window.Type”的 enum 值指定。此值连同所需的点数一起传递给函数 Coefficients。结果是一个 double 数组,可以通过与输入数据相乘来应用,也可以由其中一个比例因子例程使用来确定窗增益或损耗。窗函数名称和参数在附录 A 中指定。

窗函数分析与缩放

double DSPLib.Window.Signal(double[] winCoeffs)

double DSPLib.Window.Noise(double[] winCoeffs, double samplingFrequencyHz)

double DSPLib.Window.NENBW(double[] winCoeffs)

给定任何一组窗函数系数,这些例程将确定所提供窗函数数据数组的增益或损耗。窗函数系数可以来自任何窗函数,它们不限于 DSPLib 提供的窗函数。

Noise 函数还考虑了以 Hz 为单位的 bin 宽度,以校正多赫兹或小数赫兹的 bin 宽度。有关更多信息,请参阅示例 6。

函数 NENBW 计算所提供窗函数系数在 bin 中的归一化等效噪声带宽。

信号生成

这些例程可用于生成测试信号以进行仿真。

DSPLib.Generate.LinSpace(double startVal, double stopVal, UInt32 points)

LinSpace 生成一个数据数组,其中包含 startValStopVal 之间线性间隔的值,点数为:points。(类似于 Octave / Matlab 中同名函数)。

double[] DSPLIb.Generate.ToneSampling(double amplitudeVrms, double frequencyHz, 
  double samplingFrequencyHz, UInt32 points, double dcV = 0.0, double phaseDeg = 0)

double[] DSPLIb.Generate.ToneCycles(double amplitudeVrms, double cycles, 
         UInt32 points, double dcV = 0.0, double phaseDeg = 0)

ToneSamplingToneCycles 生成一个实数 double[] 数组,其中包含指定点数的正弦波音信号。

  • ToneSampling 根据典型的采样属性生成信号。
  • ToneCycles 根据点值生成音调。

可选参数为

  • dcV:添加到信号的直流电压,可以是正值或负值
  • PhaseDeg:添加到生成信号的相角(以度为单位)。可以使用不同的相角生成多个信号,然后使用包含的数组数学函数将它们相加,以进行进一步分析。
double[] DSPLIb.Generate.NoisePsd(double amplitudePsd, 
                                  double samplingFrequencyHz, UInt32 points)

double[] DSPLIb.Generate.NoiseRms(double amplitudeRms, double samplingFrequencyHz, 
                                  UInt32 points, double dcV)

噪声例程根据应用要求生成校准过的噪声信号。可以指定功率谱密度 (PSD),也可以使用噪声的 RMS 值。

PSD 的单位是 Vrms / rt-Hz (均方根电压 / 根赫兹)。

NoiseRMS 有一个可选值 dcV,用于向生成的噪声信号添加直流电压。默认值 = 0.0 伏特。注意:直流值不包含在 RMS 电压计算中。

数组数学辅助函数

注意:所有函数均定义为静态。

double[] DSPLib.Math.Add(double[]a, double[]b)

double[] DSPLib.Math.Subtract(double[]a, double[]b)

double[] DSPLib.Math.Multiply(double[]a, double[]b)

double[] DSPLib.Math.Divide(double[]a, double[]b)

double[] DSPLib.Math.Add(double[]a, double b)

double[] DSPLib.Math.Subtract(double[]a, double b)

double[] DSPLib.Math.Multiply(double[]a, double b)

double[] DSPLib.Math.Divide(double[]a, double b)

double[] DSPLib.Math.RemoveMean(double[] a)

double[] DSPLib.Math.Sqrt(double[] a)

double[] DSPLib.Math.Square(double[] a)

double[] DSPLib.Math.Multiply(Complex[] a, double[] b)

DSPLib 定义了一组有用的数学函数,允许 double[] 数组与另一个 double[] 数组或一个常量值进行加、减、乘或除运算。

RemoveMean 将从输入数据数组中移除平均值。这对于从输入数据数组中移除直流值非常有用。

SqrtSquare 是用于对输入数据数组进行平方或平方根运算的通用例程。

注意:在 Debug 构建模式下,如果 double[] a 数组的长度与 double[] b 数组的长度不匹配,则会断言警告。

将 DSPLib 添加到您的项目

DSPLib 需要 .NET 4 及以上版本。

  • 从上面的链接下载 DSPLib_source.zip。此 zip 文件只包含 DSPLib.cs 源文件。
  • DSPLib.cs 源文件作为现有文件添加到您的项目。
  • DSPlib 需要引用 System.Numerics 库中的复数数据类型,因此 System.Numerics 也必须添加到您的项目中作为引用。
  • 为了方便起见,将 "using DSPLib;" 添加到任何需要访问 DSPLib 命名空间的文件顶部。

可以使用上面的链接下载示例测试项目。测试应用程序(和解决方案)是用 VS2022 编写的。它应该可以直接编译。如果您有任何引用问题,请将 System.NumericsSystem.Windows.Forms.DataVisulization(用于 .NET 图表控件)添加到 DSPLib_Test 项目中。

这些示例位于 DSPLib_Test.zip 文件中,作为一个单独的文件,名为 DSPLib_ExampleCode.cs。此文件包含许多如上所述的示例使用场景,并将回答许多常见的用法问题。

优化

现在是 64 位世界,大多数 Windows 计算机都运行某种形式的 64 位操作系统。根据我的测试:将 CPU 类型设置为 x64 进行编译不会加快任何基本的 DSPLib 操作,但会增加带缓存的 DFT 的最大大小。缓存是 DFT.Initialize() 例程预先计算并将正弦和余弦项存储在数组中以用于 DFT 执行的地方,从而使计算时间减少 3 倍。缓存发生在 x32 或 x64 位模式下,但在 x64 位模式下运行时,可用于缓存的内存可以翻倍(请参阅上面 DFT.Initialize() 的描述)。

同样,将编译器选项设置为 Release(或 Optimize)而不是 Debug,无论在 x32 还是 x64 模式下,都能立即提供 3 倍的速度提升。

关于优化的进一步说明:由于该库广泛使用了 .NET 4 的“并行任务线程”过程,并且这些过程是“核心感知”的,因此处理器中运行代码的任何核心/线程的增加都将立即提高执行速度,而无需用户干预或重新编译。例如(见下图),我在三款不同的 Intel I7 处理器上,使用从 15 到 45 瓦、从 2 核到 6 核的各种 DFT 大小运行了相同的 DFT 代码。可以看出,从最慢的处理器到最快的处理器,计算速度提高了 4 倍以上,而编译代码完全没有改变。

dftspeed

三款不同处理器运行相同可执行文件的速度比较。DFT 大小从 2000 点更改为 60000 点(X 轴),并记录了计算时间(毫秒,Y 轴)。蓝色曲线是 2 核/4 线程处理器,橙色曲线是 4 核/8 线程处理器,黄色曲线是 6 核/12 线程处理器。速度提升是自动发生的,无需重新编译或用户干预。

参考文献

附录 A – 包含的窗函数

下面是 DSPLib 中实现的窗函数的性能参数表。

有关这些窗函数的更多信息,请参阅参考资料 [1]。

名称 峰值旁瓣电平 (dB) 平坦度 (dB) 3dB 带宽 (bin) NENBW (bin)
矩形或无 13.3 3.9224 0.8845 1
Welch 21.3 2.2248 1.1535 1.2
Bartlett 26.5 1.8242 1.2736 1.3333
Hanning 或 Hann 31.5 1.4236 1.4382 1.5
Hamming 42.7 1.7514 1.3008 1.3628
Nuttall3 46.7 0.863 1.8496 1.9444
Nuttall4 60.9 0.6184 2.1884 2.31
Nuttall3A 64.2 1.0453 1.6828 1.7721
Nuttall3B 71.5 1.1352 1.6162 1.7037
Nuttall4A 82.6 0.7321 2.0123 2.1253
BH92 92 0.8256 1.8962 2.0044
Nuttall4B 93.3 0.8118 1.9122 2.0212
表 1 - DSPLib 中实现的常用窗函数
名称 峰值旁瓣电平 (dB) 平坦度 (dB) 3dB 带宽 (bin) NENBW (bin)
SFT3F 31.7 0.0082 3.1502 3.1681
SFT3M 44.2 0.0115 2.9183 2.9452
FTNI 44.4 0.0169 2.9355 2.9656
SFT4F 44.7 0.0041 3.7618 3.797
SFT5F 57.3 0.0025 4.291 4.3412
SFT4M 66.5 0.0067 3.3451 3.3868
FTHP 70.4 0.0096 3.3846 3.4279
HFT70 70.4 0.0065 3.372 3.4129
FTSRS 76.6 0.0156 3.7274 3.7702
SFT5M 89.9 0.0039 3.834 3.8852
HFT90D 90.2 0.0039 3.832 3.8832
HFT95 95 0.0044 3.759 3.8112
HFT116D 116.8 0.0028 4.1579 4.2186
HFT144D 144.1 0.0021 4.4697 4.5386
HFT169D 169.5 0.0017 4.7588 4.8347
HFT196D 196.2 0.0013 5.0308 5.1134
HFT223D 223 0.0011 5.3 5.3888
HFT248D 248.4 0.0009 5.5567 5.6512
表 2 - DSPLib 中实现的平顶窗函数

窗函数以具有相同名称后缀的 ENUM 进行引用。例如,要调用“None”窗函数,ENUM 将被称为 DSPLib.DSP.Window.Type.None

类型为 HFT248D 的窗函数将被称为 DSPLib.DSP.Window.Type.HFT248D

C# 提供了多种方法用于在 ENUMstring 之间进行转换。请查阅 C# 在线帮助以获取更多信息。

注意:表 1 和表 2 中的表格值来自参考文献 [1]。

历史

  • 2016年6月18日
    • 首次发布
  • 2016年7月16日
    • 添加了更多示例
    • 添加了“PhaseUnwrap”方法
    • 修复了 DFT 相位符号
    • 修复了多次零填充执行的 FFT
    • 有关完整的更改列表,请参阅 DSPLib.cs 文件头
    • 下载文件已更新。
  • 2016年8月15日
    • 添加了优化说明
    • 添加了示例 11
    • 测试解决方案下载已更新
  • 2017年7月4日
    • 在 Log 计算中添加了 <= 零检查以捕获可能的异常
    • C# 源代码 V1.03
  • 2017年10月15日
    • V1.03 的工作方式略有改进
    • 相同的结果,不同的设计模式
  • 2018年2月9日
    • 将示例测试下载项目解决方案与最新库版本同步
    • 将测试项目更新到 VS2017
    • 无底层库或操作更改
  • 2020年9月19日
    • 添加了关于优化的更多说明
  • 2023年7月11日
    • 添加了 Complex[] 数据输入类型 FFT 和 DFT(参见示例 12)
    • 修复了 FT 返回长度与学术规范一致的问题。与此库的早期版本相比,这可能是一个突破性问题,因为返回的频谱点数可能与 1.x 版本中不同
    • 项目已更新至 .NET Framework 4.8
    • 项目已更新至 VS2022
© . All rights reserved.