目录
介绍
确实需要一个现成的傅里叶变换库,用户可以开箱即用,执行快速傅里叶变换(FFT)或离散傅里叶变换(DFT),并获得经典的频谱与频率图。
您将在商业软件包、开源库、教科书和Web上找到的绝大多数代码根本不适合这项任务,需要数小时的进一步调整才能获得经典且适当缩放的频谱图。
此处介绍的库是一个实用、有条理且完整的面向DSP的.NET 4+开源例程库,在非常非限制性的MIT许可证下发布。
DSPLib的作用
DSPLib 有几个主要部分,但其基本目标是允许在时间序列输入阵列上预形成真正的傅里叶变换,从而产生可用的经典频谱输出,而无需用户进行任何进一步的调整。
基本傅里叶变换(FT)有两种基本类型:最一般的形式可以从任何长度的输入数据中产生频谱输出。这种类型的变换称为离散傅里叶变换或DFT。代码简单且蛮力。
- 优点是: 输入数据可以是任意长度。
- 缺点是:由于它是一种通用方法,因此计算密集型,并且大型输入数据集可能需要很长时间才能计算。
一种更具体的FT类型称为快速傅里叶变换或FFT。
- 优点是:它的计算速度比DFT快得多。
- 缺点是:输入数据长度被限制为二的幂。代码更难理解。
例如:8192点FFT需要:在我的i7计算机上不到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处理器上,使用Task Parallel扩展可将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不注意到输入数据中的不连续性。
通常根据四个标准选择窗口:
- 它们推动不需要的旁瓣的振幅有多低
- 它们在多大程度上降低了振幅不确定性
- 他们将数据涂抹到多少个FT Bin中
- 历史原因,即“在这种情况下,每个人都使用Window:xyz”等。
选择最佳窗口涉及许多权衡。
DSPLib的窗口是基于马克斯·普朗克研究所[1]工作人员撰写的一篇优秀文章,包括所有常用的窗户类型,如:韦尔奇、汉明和汉恩(或通常所说的汉宁)和其他27个非常有用但可能不太知名的窗户。还包括17个非常独特的平顶窗口,即使输入频率在FT输出Bin之间,也能提供非常低的幅度误差。
由于对输入数据进行窗口化会物理上改变它,因此它也可以改变生成的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添加到您的项目”部分。
示例1——示例1代码片段的输出图
FF呢?
使用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.5Hz),DFT被实例化并使用零填充进行初始化,以显示精细的窗口细节。最后,执行DFT。对结果进行缩放,大小结果在数组中可用:lmSpectrum。
示例2 A——示例2代码片段的输出图
零填充可显示精细的光谱细节。
示例2 B——示例2代码片段的输出图
在X轴上放大后,由于应用了窗口并添加了零填充,因此显示了精细的光谱细节。另请注意,即使本例中的DFT长度为1000点+ 9000零填充=总共10,000点,振幅仍然正确。振幅校正很容易在DSPLib中进行:零填充在FT中处理,应用窗口的振幅误差通过窗口比例因子函数进行校正。
回顾使用DSPLib的基本使用步骤
DSPLib被设计为在DFT和FFT的使用上保持一致。 此外,信号和类似噪声的信号可以按照相同的基本步骤进行分析,如下所述。
基本使用步骤——信号幅度分析
- 实例化和初始化FT,计算信号的窗口比例因子
- 获取或生成输入数据
- 将窗口应用于输入数据
- 对窗口输入数据执行FT
- 转换为幅度格式
- 如果需要:逐个箱中平均,然后对所有平均值重复步骤2-6
- 将窗口比例因子应用于数据
结果:正确缩放的信号频谱
基本使用步骤——噪声信号分析
- 实例化和初始化FT,计算噪声的窗口比例因子
- 获取或生成输入数据
- 将窗口应用于输入数据
- 对窗口输入数据执行FT
- 转换为幅度平方格式
- 如果需要:逐个条柱对数据条柱进行平均量级平方,然后对所有平均值重复步骤2-6
- 转换为幅度格式
- 将窗口比例因子应用于数据。
结果:正确缩放噪声信号频谱。
注:参考文献[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也可以在Back Ground Worker线程上运行
- 示例12——复杂输入数据(I&Q)的DFT,还显示了如何窗口化IQ数据以及如何应用比例因子
完整的库文档
下面介绍了每个部分和DSPLib功能的完整库说明。该库还实现了在Visual Studio环境中启用IntelliSense帮助的XML注释。
傅里叶变换
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[]数组
- 仅返回频谱的实部(输入数据长度的一半)。这对应于:直流通过采样频率/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:表示以log10格式表示以1伏RMS为基准的幅度。 例如:1 Vrms = 0dBV,0.1Vrms = -20dBV
- ToMagnitudeSquared:表示频谱的幅度平方值。幅度平方(V2)与功率成正比。由于噪声与功率成正比,因此在平均噪声时,幅度平方是正确的格式。
- ToPhaseDegrees:表示每个点处复数的相位(以度为单位)
- ToPhaseRadians:表示每个点处复数的相位(以弧度为单位)
double[] DSPLib.ConvertMagnitude.ToMagnitudeSquared(double[] spectrum)
double[] DSPLib.ConvertMagnitude.ToMagnitudeDBV(double[] spectrum)
这组函数将FT转换后的频谱输出的幅度double[]数组转换为更方便的最终使用格式。有几种输出格式可供选择:
- ToMagnitude:表示每个条柱处复数的大小。
- ToMagnitudeDBV:表示以Log10格式表示以1伏RMS为基准的幅度。 例如:1 Vrms = 0dBV,0.1Vrms = -20dBV
- ToMagnitudeSquared:表示频谱的幅度平方值。幅度平方(V2)与功率成正比。由于噪声与功率成正比,因此在平均噪声时,幅度平方是正确的格式。
double[] DSPLib.ConvertMagnitudeSquared.ToMagnitude(double[] spectrum)
double[] DSPLib.ConvertMagnitudeSquared.ToMagnitudeDBV(double[] spectrum)
将FT转换后的频谱输出的MagnitudeSquareddouble[]阵列转换为更方便的最终使用格式的功能。有几种输出格式可供选择:
- ToMagnitude:表示每个点处复数的大小
- ToMagnitudeDBV:表示以Log10格式表示以1伏RMS为基准的幅度。 例如:1 Vrms = 0dBV,0.1Vrms = -20dBV
频谱数据分析功能
double DSPLib.Analyze.FindMaxAmplitude(double[] inData)
double DSPLib.Analyze.FindMaxPosition(double[] inData)
double DSPLib.Analyze.FindMaxFrequency(double[] inData, double[] fSpan)
给定一个输入数组,这些函数会找到指定的最大数量。在阵列中找到的最大振幅、找到最大振幅的箱位置或找到最大振幅的频率。这些函数设计用于频谱输出,并直接接受任何格式的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)
给定输入数组,这些函数确定数组的平均值(平均值)或RMS值。由于DC或Fs/2(最后一个数据点)存在泄漏到频谱中,因此可选参数允许从结果中排除第一个和最后一个指定的bin。两端的默认值为10个bins,足以排除随附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值指定窗口这与所需的点数一起传递到函数系数中。结果是一个双精度数组,可以通过与输入数据相乘来应用,也可以由其中一个比例因子例程用于确定窗口增益或损失。窗口名称和参数在附录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函数计算以箱为单位提供的窗口系数的归一化等效噪声带宽。
信号生成
这些例程可用于生成用于仿真目的的测试信号。
DSPLib.Generate.LinSpace(double startVal, double stopVal, UInt32 points)
LinSpace生成一个数据数组,其线性间隔值介于startVal、StopVal之间,点数为:点。(类似于同名的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)
ToneSampling和ToneCycles生成具有指定点数的正弦波音信号的真实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(Volts RMS / Root-Hertz)。
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将从输入数据数组中删除平均值或平均值。这对于从输入数据数组中删除DC值很有用。
Sqrt和Square且是输入数据数组的平方根或平方根的通用例程。
注意:在调试构建模式下,如果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.Numerics和System.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()说明)。
同样,在x32或x64模式下,将编译器选项设置为 Release(或Optimize)而不是Debug可立即将速度提高3倍。
关于优化的进一步说明:由于该库广泛使用了.NET 4“并行任务线程”过程,并且这些过程是“核心感知”的,因此在运行代码的处理器中添加任何核心/线程都将立即提高执行速度,而无需任何用户干预或重新编译。举个例子,我在三个不同的英特尔I7处理器上运行了具有不同DFT大小的相同DFT代码,功率从15瓦到45瓦,从2核到6核。可以看出,从最慢的处理器到最快的处理器,计算速度提高了4倍以上,而编译后的代码完全没有变化。
运行同一可执行文件的三个不同处理器的速度比较。DFT大小从2000点更改为60000点(X轴),并记录以毫秒为单位的计算时间(Y轴)。蓝色跟踪是2核/4线程处理器,橙色跟踪是4核/8线程处理器,黄色跟踪是6核/12线程处理器。速度提升会自动发生,无需重新编译或用户干预。
引用
- [1] G. Heinzel、A. Rudiger和R. Schilling,“通过离散傅里叶变换(DFT)进行光谱和光谱密度估计,包括窗口函数和一些新的平顶窗口的综合列表”,Max-Planck-Institut fur Gravitationsphysik,2002年2月15日
- [2] Steve Hageman,“傅里叶变换的实践仪器工程师指南”,2012年6月EDN.com日发布。
- [3] DSPLib中使用的内核FFT源自Gerald T. Beauregard在2010年撰写的一篇博客文章。它根据其发布的MIT许可证的规定在此处使用。
附录A——包含的窗口
以下是DSPLib中实现的窗口的表格性能参数。
有关这些窗口的详细信息,请参阅参考[1]。
名字 | 峰值边带电平 (dB) | 平坦度(dB) | 3dB带宽(bins) | NENBW(bins) |
Rectangular or None | 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 or 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带宽(bins) | NENBW(bins) |
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#有许多将ENUM转换为string的方法,反之亦然。有关详细信息,请参阅C#联机帮助。
注:表1和表2中的表格值来自参考文献[1]。
https://www.codeproject.com/Articles/1107480/DSPLib-FFT-DFT-Fourier-Transform-Library-for-NET-6