Alby's blog

世上没有巧合,只有巧合的假象。

0%

在 WPF 中绘制功率谱(Power Spectrum)

一、概述

本文简述在 WPF 中绘制功率谱(Power Spectrum)。会将功率谱密度(PSD,Power spectral density)单位转换为分贝(dB, Decibel), 实现类似于 Adobe Audition 软件中频谱分析的功能。为简单起见及本人需要,FFT 数量使用 2048,窗口函数使用汉宁窗(Hanning)

二、准备音频文件

为了方便观察,生成一条采样率为 44.1K、单声道、位深 32 bit Float 的,时长为 20 秒的,频率为 40Hz、音量为 0 dB 的正弦波形的,封装格式为 wav 的音频。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
var signalGenerator = new SignalGenerator()
{
Gain = 1.0,
Frequency = 40,
Type = SignalGeneratorType.Sin,
};
var sine20Seconds = signalGenerator.Take(TimeSpan.FromSeconds(20));
using (var wo = new WaveFileWriter(@"C:\your\wav\file.wav", signalGenerator.WaveFormat))
{
var buffer = new float[signalGenerator.WaveFormat.SampleRate];
var count = 0;
while((count = sine20Seconds.Read(buffer, 0, buffer.Length)) > 0)
{
wo.WriteSamples(buffer, 0, count);
}
}


图1:波形图

三、读取音频文件

读取文件使用 NAudo,一是因为其对常见音频格式支持也算友好,二是在 API 方面也比 Bass(ManagedBass) 等其他库更 C# 些,三是在 .Net 下并且在当前场景下用 FFmpeg 不太明智。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
using NAudio.Wave;

/// <summary>
/// 从单声道 Wav 文件中读取指定数量的采样
/// <para>Note: 如果对多声道音频文件进行操作,可进行改写。</para>
/// </summary>
/// <param name="wavFile"></param>
/// <param name="numberOfSamples"></param>
/// <returns></returns>
private static float[] ReadMonoWavFile(string wavFile, int numberOfSamples = 2048)
{
using (var wfr = new WaveFileReader(wavFile))
{
var data = new float[numberOfSamples];
for (var i = 0; i < numberOfSamples; i++)
{
data[i] = wfr.ReadNextSampleFrame()[0];
}
return data;
}
}

四、计算 Power(纵坐标)

计算 FFT 使用 Math.Net。做信号分析一般使用 MATLAB 和 Python 用得多一些,在信号处理方面 .Net 中工具支持相对较弱。NAudio 及其他库也可以计算 FFT,但 Math.Net 比 NAudio 等更完善一些。

备注:开源项目 FftSharp 看其源码发现 Window 函数的计算公式有误,直接抛弃。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
using System.Numerics;
using MathNet.Numerics.IntegralTransforms;

const int FFTLength = 2048;

// 0. 准备数据
// Note: 读取 2048 个采样
var samples = new float[FFTLength];

// 1.使用窗口函数创建窗口
var window = MathNet.Numerics.Window.Hann(FFTLength);

// 2.将 float 采样转换为虚数为0的复数,并能应用窗口(点积)
// Note: 因为 Math.Net 没有类似于 numpy.fft.rfft 的方法,故需先转换成复数。fft 变量作为 Fourier.Forward 的输入和输出。
var fft = new Complex[numberOfSamples];
for (int i = 0; i < numberOfSamples; i++)
{
fft[i] = new Complex(samples[i] * window[i], 0.0);
}

// 3. FFT
Fourier.Forward(fft, FourierOptions.Matlab);

// 4. 获取音量
// Note: referenceValue: 参考值。用于 dB FS 拉伸。int16 是 32768, float 是 1。
var dbFS = new double[realFFTCount];
}
for (var i = 0; i < realFFTCount; i++)
{
// Scale the magnitude of FFT by window and factor of 2, because we are using half of FFT spectrum.
var abs = Complex.Abs(fft[i]);
var magnitude = abs * 2 / window.Sum();
// Convert to dBFS
dbFS[i] = 20 * Math.Log10(magnitude / 1.0/*referenceValue*/);
}

五、计算频率(横坐标)

根据生成的 FFT 的数量生成一个集合。如果是 FFT 长度是 2048,则生成 1024 个区间共 1025 个点;每个点的值表示频率明确标识是哪个 Hz —— 想象5根手指4条逢。

1
2
var realFFTCount = FFTLength / 2 + 1;
var frequencies = Enumerable.Range(0, realFFTCount).Select(m => (double)m / FFTLength / SampleRate).ToArray();

六、封装

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
public class PowerSpectrumGenerator
{
private double[] window;
private double windowSum;
private Complex[] fft;
private double[] dbFS;

private int _sampleRate;
public int SampleRate
{
get
{
return _sampleRate;
}
set
{
if (value != _sampleRate)
{
_sampleRate = value;
ResetSampleRateOrFFTLength();
}
}
}

private int _fftLength;
public int FFTLength
{
get
{
return _fftLength;
}
set
{
if (value != _fftLength)
{
_fftLength = value;
ResetSampleRateOrFFTLength();
}
}
}

public double[] Frequencies { get; private set; }

public PowerSpectrumGenerator()
{

}

public PowerSpectrumGenerator(int sampleRate, int fftLength)
{
_sampleRate = sampleRate;
_fftLength = fftLength;
ResetSampleRateOrFFTLength();
}

/// <summary>
/// 计算 Power
/// </summary>
/// <param name="width"></param>
/// <returns></returns>
public double[] Compute(float[] samples)
{
if (samples == null || samples.Length != FFTLength)
{
throw new ArgumentException("采样数请保持和 fftLength 长度一致");
}

//referenceValue: 参考值。用于 dB FS 拉伸。int16 是 32768, float 是 1。
const double referenceValue = 1.0;

var numberOfSamples = samples.Length;
var realFFTCount = numberOfSamples / 2 + 1;

// Note: fft 变量可复用,不用每次构造。
if (fft == null || fft.Length < numberOfSamples)
{
fft = new Complex[numberOfSamples];
}

// 该 for 循环有两个作用:
// 1. 将采样集合和窗口进行点积(Dot product);
// 2. Math.Net 没有类似于 numpy.fft.rfft 方法,故构造虚数为 0 的复数。
for (int i = 0; i < numberOfSamples; i++)
{
fft[i] = new Complex(samples[i] * window[i], 0.0);
}

// fft 变量作为Fourier.Forward的输入和输出。
Fourier.Forward(fft, FourierOptions.Matlab);

if(dbFS == null || dbFS.Length < realFFTCount)
{
dbFS = new double[realFFTCount];
}
for (var i = 0; i < realFFTCount; i++)
{
// Scale the magnitude of FFT by window and factor of 2, because we are using half of FFT spectrum.
var abs = Complex.Abs(fft[i]);
var magnitude = abs * 2 / windowSum;

// Convert to dBFS
dbFS[i] = 20 * Math.Log10(magnitude / referenceValue);
}

return dbFS;
}

private void ResetSampleRateOrFFTLength()
{
var realFFTCount = _fftLength / 2 + 1;
Frequencies = Enumerable.Range(0, realFFTCount).Select(m => (double)m / ((double)_fftLength / _sampleRate)).ToArray();

window = MathNet.Numerics.Window.Hann(_fftLength);
windowSum = window.Sum();
}
}

七、绘图

绘图采用 OxyPlot 。如果要更好性能,建议用 DirectX 绘制。
在获取横纵两组数据后,就可以使用 OxPlot 绘制出来了。但获取的数据封装格式不符合 OxPlot,做些小调整即可。

1、创建 ViewModel

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
public class MainViewModel : INotifyPropertyChanged
{
private LineSeries _lineSeries = new LineSeries { MarkerType = MarkerType.Circle, Color = OxyColor.Parse("#00FF00") };

public MainViewModel()
{
_mainModel = new PlotModel { Title = "Spectrum Analyzer", TitleColor = OxyColor.Parse("#E8E8E8"), Background = OxyColor.Parse("#000000") };
_mainModel.Axes.Add(new LinearAxis
{
Position = AxisPosition.Bottom,
AxislineStyle = LineStyle.Solid,
AxislineColor = OxyColor.Parse("#E8E8E8"),
MajorGridlineStyle = LineStyle.Solid,
MajorGridlineColor = OxyColor.Parse("#E8E8E8"),
MinorGridlineStyle = LineStyle.Dot,
MinorGridlineColor = OxyColor.Parse("#E8E8E8"),
TextColor = OxyColor.Parse("#E8E8E8"),
TickStyle = TickStyle.None,
Minimum = 0,
Maximum = 22050,
IsZoomEnabled = false,
Title = "Frequencies(Hz)",
TitleColor = OxyColor.Parse("#00CDCD"),
});
_mainModel.Axes.Add(new LinearAxis
{
Position = AxisPosition.Left,
AxislineStyle = LineStyle.Solid,
AxislineColor = OxyColor.Parse("#E8E8E8"),
MajorGridlineStyle = LineStyle.Solid,
MajorGridlineColor = OxyColor.Parse("#E8E8E8"),
MinorGridlineStyle = LineStyle.Dot,
MinorGridlineColor = OxyColor.Parse("#E8E8E8"),
TextColor = OxyColor.Parse("#E8E8E8"),
TickStyle = TickStyle.None,
Minimum = -200,
Maximum = 20,
AbsoluteMinimum = -214,
AbsoluteMaximum = 44,
IsZoomEnabled = true,
Title = "Amplitude(dB FS)",
TitleColor = OxyColor.Parse("#00CDCD")
});
_mainModel.Series.Add(new LineSeries
{
MarkerSize = 1,
Color = OxyColor.Parse("#00FF00")
});
_mainModel.Series.Add(_lineSeries);
}

private PlotModel _mainModel;
public PlotModel MainModel
{
get { return _mainModel; }
set { _mainModel = value; NotifyPropertyChanged("MainModel"); }
}

public void Update(double[] frequencies, double[] dbFS)
{
// X 坐标可能根据频率范围改变
var minFrequencie = frequencies.Min();
var maxFrequencie = frequencies.Max();
_mainModel.Axes[0].Minimum = minFrequencie;
_mainModel.Axes[0].Maximum = maxFrequencie;

// Y 坐标的最小值可能根据最小音量改变
var mindBFS = dbFS.Min();
_mainModel.Axes[1].Minimum = mindBFS;

_lineSeries.Points.Clear();
for(var i = 0; i <dbFS.Length; i++)
{
var x = frequencies[i];
var y = dbFS[i];
_lineSeries.Points.Add(new DataPoint(x, y));
}

NotifyPropertyChanged("MainModel");
}

#region NotifyPropertyChanged

public event PropertyChangedEventHandler PropertyChanged;

private void NotifyPropertyChanged(string info)
{
PropertyChanged?.Invoke(this, new PropertyChangedEventArgs(info));
}

#endregion

2、修改 MainWindow.xaml

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
<Window x:Class="Tubumu.Audio.SpectrumAnalyzer.MainWindow"
xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
xmlns:d="http://schemas.microsoft.com/expression/blend/2008"
xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006"
xmlns:local="clr-namespace:Tubumu.Audio.SpectrumAnalyzer"
xmlns:viewModels="clr-namespace:Tubumu.Audio.SpectrumAnalyzer.ViewModels"
xmlns:oxy="clr-namespace:OxyPlot.Wpf;assembly=OxyPlot.Wpf"
mc:Ignorable="d"
Title="Spectrum Analyzer" Height="600" Width="800" Loaded="Window_Loaded">
<Grid Background="Black">
<Grid.ColumnDefinitions>
<ColumnDefinition/>
<ColumnDefinition Width="20"/>
</Grid.ColumnDefinitions>
<oxy:PlotView Model="{Binding MainModel}" Grid.Column="0"/>
</Grid>
</Window>

3、修改 MainWindow.xaml.cs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
/// <summary>
/// MainWindow.xaml 的交互逻辑
/// </summary>
public partial class MainWindow : Window
{
private PowerSpectrumGenerator _powerSpectrumGenerator = new PowerSpectrumGenerator();
private MainViewModel _mainViewModel = new MainViewModel();

public MainWindow()
{
InitializeComponent();
DataContext = _mainViewModel;
}

private void Window_Loaded(object sender, RoutedEventArgs e)
{
var file = @"C:\your\wav\file.wav";

ReadWavFile(file);
}

/// <summary>
/// 从 Wav 文件中读取指定数量的采样
/// <para>Note: 如果对多声道音频文件进行操作,可进行改写。</para>
/// </summary>
/// <param name="wavFile"></param>
/// <param name="numberOfSamples"></param>
/// <returns></returns>
private float[] ReadWavFile(string wavFile, int numberOfSamples = 2048)
{
using (var wfr = new WaveFileReader(wavFile))
{
_powerSpectrumGenerator.SampleRate = wfr.WaveFormat.SampleRate;
_powerSpectrumGenerator.FFTLength = numberOfSamples;

var data = new float[numberOfSamples];
for (var i = 0; i < numberOfSamples; i++)
{
data[i] = wfr.ReadNextSampleFrame()[0];
}

var dbFS = _powerSpectrumGenerator.Compute(data);
_mainViewModel.Update(_powerSpectrumGenerator.Frequencies, dbFS);

return data;
}
}
}

八、效果图


图2:本文生成


图3:Audition生成

九、包依赖

  1. NAudo
  2. Math.Net
  3. System.Runtime.Numerics

参考资料

https://dsp.stackexchange.com/questions/32076/fft-to-spectrum-in-decibel
https://blog.csdn.net/weixin_42930928/article/details/89143908
https://www.cnblogs.com/luge/p/Wpf_Oxyplot_Sample.html