Получение частоты звука с микрофона

#c_sharp #аудио #фурье #fft #naudio

Есть задача получать частоту звука с микрофона для дальнейших преобразований на C#.
Подобное уже делал на Python с numpy, но тут как то не клеится...

private void write(byte[] angles, int byte_len)
{
    Complex[] _fftBuffer = new Complex[byte_len];
    var _m = (int)Math.Log(byte_len, 2.0);

    for (var n=0; n < byte_len; n++)
    {
        var r = angles[n];
        var i = 0;
        _fftBuffer[n].X = (float)(r * FastFourierTransform.HammingWindow(n, byte_len));
        _fftBuffer[n].Y = i;

    }
    FastFourierTransform.FFT(true, _m, _fftBuffer);
    float[] fft_x = new float[_fftBuffer.Length];
    for (var i=0; i<_fftBuffer.Length; i++)
    {
        fft_x[i] = Math.Abs(_fftBuffer[i].X);
    }
    int i_peak = fft_x.ToList().IndexOf(fft_x.Max());
    for (var i = 0; i < _fftBuffer.Length; i++)
    {
        fft_x[i] = (float)Math.Log(Math.Abs(_fftBuffer[i].X));
    }
    ///var i_peak = fft_x.Max();
    var i_interp = parabolic(fft_x, i_peak);

    float freq = byte_len * i_interp / angles.Length;
    Console.WriteLine("Debug stop");

    ///this.port.Write(this.first_command, 0, 8);
    ///this.port.Write(this.get_comand(0.0f, 0.0f), 0, 32);


}
private float parabolic(float[] f, int peak)
{
    var xv = 0.5f * (f[peak-1] - f[peak+1])/(f[peak-1]-2*f[peak] + f[peak+1]) + peak;

    return xv;
}


Это что на С# сделал. Сделано так, ибо так же было на Python.
В конечном итоге работает, но неправильно - частота получается одна и та же (+- пара
герц), но по идее там должны быть абсолютно другие частоты - передаем данные с помощью
звука.

Помогите разобраться, что в коде может быть не так?

UPD

Вот что получилось в итоге. Спасибо товарищу @MSDN.WhiteKnight - натолкнул на правильные
мысли.
Плюс использовался проект 
вот отсюда

выкладываю только код, который несет смысл по вытаскиванию частот звука с микрофона(МОНО).
Можно переделать на стерео - не особо сложно будет

public partial class MainWindow : Window
{
    static double Fs = 48000; // Частота дискретизации !В данной программе ТОЛЬКО
целые числа
    static double T = 1.0 / Fs; // Шаг дискретизации
    static int N; //Длина сигнала (точек)
    static double Fn = Fs / 2;// Частота Найквиста
    WaveIn waveIn;

    public MainWindow()
    {
        InitializeComponent();
    }

    void waveIn_DataAvailable(object sender, WaveInEventArgs e)
    {

        //данные из буфера распределяем в массив чтобы в нем они были в формате ?PCM?
        byte[] buffer = e.Buffer;
        N = buffer.Length;
        int bytesRecorded = e.BytesRecorded;
        Complex[] sig = new Complex[bytesRecorded / 2];
        for (int i = 0, j = 0; i < e.BytesRecorded; i += 2, j++)
        {
            short sample = (short)((buffer[i + 1] << 8) | buffer[i + 0]);
            sig[j] = sample / 32768f;
        }

        Fourier.Forward(sig, FourierOptions.Matlab);
        // обнуляем спектр на небольших частотах (там постоянная составляющая и вообще
много помех)
        for (int i = 0; i < 35 * sig.Length / Fn; i++)
        {
            sig[i] = 0;
        }

        write(sig);

    }
    //Окончание записи
    private void waveIn_RecordingStopped(object sender, EventArgs e)
    {
        waveIn.Dispose();
        waveIn = null;
    }

    private void start_button_Click(object sender, RoutedEventArgs e)
    {
        this.waveIn = new WaveIn();
        this.waveIn.DeviceNumber = 0;
        this.waveIn.DataAvailable += this.waveIn_DataAvailable;
        this.waveIn.RecordingStopped += this.waveIn_RecordingStopped;
        this.waveIn.WaveFormat = new WaveFormat((int)Fs, 1);
        this.waveIn.StartRecording();
        Start_button.IsEnabled = false;
        Stop_button.IsEnabled = true;
        this.log_box("старт записи");
    }


    private void stop_Button_Click(object sender, RoutedEventArgs e)
    {
        this.stop_recording();
    }

    private void stop_recording()
    {
        this.waveIn.StopRecording();
        Start_button.IsEnabled = true;
        Stop_button.IsEnabled = false;
        this.log_box("конец записи");
    }


    private void log_box(string message)
    {
        Log_Box.AppendText("\n" + message);
        Log_Box.ScrollToEnd();
    }

    private void write(Complex[] signal)
    {
        PointPairList list1 = new PointPairList();
        int max_index = 0;
        double freq = 0;
        double K = signal.Length / 2;
        for (int i = 0; i < K; i++)
        {
            list1.Add(i * Fn / K, Complex.Abs(signal[i]) / N * 2);
        }

        foreach (ZedGraph.PointPair i in list1)
        {
            if (i.Y > list1[max_index].Y)
            {
                max_index = list1.IndexOf(i);
            }
        }
        freq = list1[max_index].X;

        string s = freq.ToString();
        log_box(s);


    }

}


что использовалось...
NAudio - для получения потока звука с микрофона

MathNET - Фурье

ZedGraph - как переходник для работы с сигналом после преобразования Фурье - его
можно убрать, но в моем случае был удобен.
    

Ответы

Ответ 1

Ваш код будет работать нормально, только если на вход подать данные определенного
формата: моно, 1 байт на сэмпл, определенная частота дискретизации и т.п. Кроме того,
он не учитывает несколько деталей: из результата БПФ нужно отбросить первое значение
("постоянная составляющая") и вторую половину значений (которая не несет полезной информации);
количество сэмплов должно быть в степени 2. 

Лучше написать код, который может корректно обрабатывать разные форматы, для этого
возьмем за основу класс SampleAggregator из примера на Github:

using System;
using System.Collections.Generic;
using System.Text;
using System.Diagnostics;
using NAudio.Dsp;

namespace WindowsFormsTest1
{    
    public class SampleAggregator
    {
        // volume
        public event EventHandler MaximumCalculated;
        private float maxValue;
        private float minValue;
        public int NotificationCount { get; set; }
        public Complex[] FftBuffer { get { return this.fftBuffer; } }
        int count;

        // FFT
        public event EventHandler FftCalculated;
        public bool PerformFFT { get; set; }
        private Complex[] fftBuffer;
        private FftEventArgs fftArgs;
        private int fftPos;
        private int fftLength;
        private int m;

        public SampleAggregator(int fftLength = 1024)
        {
            if (!IsPowerOfTwo(fftLength))
            {
                throw new ArgumentException("FFT Length must be a power of two");
            }
            this.m = (int)Math.Log(fftLength, 2.0);
            this.fftLength = fftLength;
            this.fftBuffer = new Complex[fftLength];
            this.fftArgs = new FftEventArgs(fftBuffer);
        }

        bool IsPowerOfTwo(int x)
        {
            return (x & (x - 1)) == 0;
        }


        public void Reset()
        {
            count = 0;
            maxValue = minValue = 0;
        }

        public void Add(float value)
        {
            if (PerformFFT)
            {
                fftBuffer[fftPos].X = (float)(value * FastFourierTransform.HammingWindow(fftPos,
fftBuffer.Length));
                fftBuffer[fftPos].Y = 0;
                fftPos++;
                if (fftPos >= fftBuffer.Length)
                {
                    fftPos = 0;
                    // 1024 = 2^10
                    FastFourierTransform.FFT(true, m, fftBuffer);

                    if(FftCalculated != null) FftCalculated(this, fftArgs);
                }
            }

            maxValue = Math.Max(maxValue, value);
            minValue = Math.Min(minValue, value);
            count++;
            if (count >= NotificationCount && NotificationCount > 0)
            {
                if (MaximumCalculated != null)
                {
                    MaximumCalculated(this, new MaxSampleEventArgs(minValue, maxValue));
                }
                Reset();
            }
        }
    }

    public class MaxSampleEventArgs : EventArgs
    {
        [DebuggerStepThrough]
        public MaxSampleEventArgs(float minValue, float maxValue)
        {
            this.MaxSample = maxValue;
            this.MinSample = minValue;
        }
        public float MaxSample { get; private set; }
        public float MinSample { get; private set; }
    }

    public class FftEventArgs : EventArgs
    {
        [DebuggerStepThrough]
        public FftEventArgs(Complex[] result)
        {
            this.Result = result;
        }
        public Complex[] Result { get; private set; }
    }
}


Тогда для определения частоты порции из первых 1024 сэмплов Wav-файла можно использовать
вот такой код:

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Linq;
using System.Text;
using System.Windows.Forms;

using NAudio;
using NAudio.Wave;
using NAudio.Wave.SampleProviders;

namespace WindowsFormsTest1
{
    public partial class Form1 : Form
    {        

        private float parabolic(float[] f, int peak)
        {
            if (peak == 0) return f[0];
            var xv = 0.5f * (f[peak - 1] - f[peak + 1]) / (f[peak - 1] - 2 * f[peak]
+ f[peak + 1]) + peak;

            return xv;
        }       


        public Form1()
        {
            InitializeComponent();            
        }

        void PrintFrequency(float[] samples, int n_samples, WaveFormat fmt)
        {
            textBox1.Text = "";
            for (int i = 0; i < fmt.Channels; i++)
            {
                SampleAggregator aggregator = new SampleAggregator(n_samples);
                aggregator.PerformFFT = true;

                int j;
                float f;

                //выделяем данные одного канала
                for (j = 0; j < n_samples; j++)
                {
                    int index = (j * fmt.Channels) + i;
                    f = samples[index];
                    aggregator.Add(f);
                }

                float[] fft_x = new float[aggregator.FftBuffer.Length / 2]; //только
первая половина БПФ имеет смысл

                for (j = 0; j < fft_x.Length; j++)
                {
                   float real = aggregator.FftBuffer[j].X;
                   float imag = aggregator.FftBuffer[j].Y;
                   fft_x[j] = (float)Math.Sqrt(real * real + imag * imag); //получаем
амплитуду  
                }

                fft_x[0] = 0.0f;//избавляемся от постоянной составляющей

                int i_peak = fft_x.ToList().IndexOf(fft_x.Max());                

                for (j = 0; j < fft_x.Length; j++)
                {
                    fft_x[j] = (float)Math.Log(Math.Abs(aggregator.FftBuffer[j].X));
                }

                var i_interp = parabolic(fft_x, i_peak);

                float freq = fmt.SampleRate * i_interp / (float)n_samples;

                textBox1.Text += ("Channel " + i.ToString() + ": " + freq.ToString()
+ " Hz" + Environment.NewLine);
            }

        }


        private void button1_Click(object sender, EventArgs e)
        {

            WaveStream readerStream = new WaveFileReader("c:\\Test\\sound_01.wav");
            WaveStream pcmStream;
            WaveStream stream;

            //создаем поток в PCM-формате
            if (readerStream.WaveFormat.Encoding != WaveFormatEncoding.Pcm)
            {
                pcmStream = WaveFormatConversionStream.CreatePcmStream(readerStream);
                stream = new BlockAlignReductionStream(pcmStream);
            }
            else
            {
                pcmStream = readerStream;
                stream = readerStream;
            }

            float[] samples;
            const int N_SAMPLES = 1024; //количество сэмплов для спектрального анализа
            ISampleProvider prov;

            using(stream)
            using(readerStream)
            using (pcmStream)
            {
                prov = stream.ToSampleProvider();                
                samples = new float[N_SAMPLES * prov.WaveFormat.Channels];

                int res = prov.Read(samples, 0, N_SAMPLES * prov.WaveFormat.Channels);
                if (res < N_SAMPLES * prov.WaveFormat.Channels) throw new Exception("Not
enough data");

            }

            PrintFrequency(samples,N_SAMPLES,prov.WaveFormat);

        }  

    }

}