在 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 的音频。

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 不太明智。

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 函数的计算公式有误，直接抛弃。

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条逢。

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

六、封装

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

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

<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

/// <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