C#实现和验证基于频率采样法的FIR滤波器（参考Matlab fir2函数实现）

关注

原创 2023-11-13 09:04:23 641 阅读 AI 写同款· GEO 优化 ›

从WFIR到WIIR：基于频率弯折的音频滤波器设计实战与MATLAB实现本文深入探讨了基于频率弯折技术的WFIR和WIIR音频滤波器设计，通过MATLAB实现详细展示了从理论到实战的全过程。文章重点解析了频率规整技术如何模拟人耳非线性听觉特性，并提供了优化计算效率和性能的实用技巧，适用于HRTF建模、语音增强等多种音频处理场景。阅读详情

具体代码

理论知识看各个课本，都有介绍，这里省略。
参考的是Matlab中fir2函数，参考【Matlab帮助中心】fir2函数
在这里插入图片描述

public static double[] fir2(int nn, double[] freqs, double[] mags)
{
    int npt = 512;
    int lap = 0;
    double[] wind;
    nn++;
    if (nn < 1024)
    {
        npt = 512;
    }
    else
    {
        npt = (int)(2 * Math.Ceiling(Math.Log(nn) / Math.Log(2)));
    }

    wind = hamming(nn);

    lap = (int)Math.Floor(npt / 25.0);

    freqs[0] = 0;
    freqs[mags.Length - 1] = 1;
    npt++;

    double[] H = new double[npt];
    int nint = mags.Length - 1;
    double[] df = diff(freqs);
    double nb = 1.0;
    double ne;
    H[0] = mags[0];
    for (int i = 0; i < nint; i++)
    {
        if (df[i] == 0)
        {
            nb = Math.Ceiling(nb - lap / 2.0);
            ne = nb + lap;
        }
        else
        {
            ne = (int)Math.Floor(freqs[i + 1] * npt);
        }

        for (int j = (int)nb; j <= ne; j++)
        {
            double inc = 0;
            if (nb == ne)
            {
                inc = 0;
            }
            else
            {
                inc = (j - nb) / (ne - nb);
            }

            if (j <= H.Length)
            {
                H[j - 1] = inc * mags[i + 1] + (1.0 - inc) * mags[i];
            }
        }
        nb = ne + 1;
    }

    double dt = 0.5 * (nn - 1);
    Complex[] rad = new Complex[npt];
    for (int i = 0; i < npt; i++)
    {
        rad[i] = -dt * Complex.Sqrt(-1.0) * Math.PI * i / (npt - 1.0);
    }

    Complex[] Hc = new Complex[npt];

    for (int i = 0; i < npt; i++)
    {
        Hc[i] = H[i] * Complex.Exp(rad[i]);
    }

    Complex[] newH = new Complex[(npt - 1) * 2];


    newH[0] = Hc[0];
    newH[npt - 1] = Hc[npt - 1];
    for (int i = 0; i < (npt - 1) * 2; i++)
    {
        if (i < npt)
        {
            newH[i] = Hc[i];
        }
        else
        {
            newH[i] = Complex.Conjugate(Hc[(npt-1)*2-i]);//取共轭复数
        }
    }

    var temp = FFTHelper.iFFT(newH, newH.Length); //做快速傅里叶逆变换

    var ht = temp.Select(x => x.Real).ToArray(); //取实部数据

    double[] result = new double[nn];
    Array.Copy(ht, result, nn);  //取前nn个数据
    
    for (int i = 0; i < nn; i++)
    {
        result[i] *= wind[i];
    }

    return result;
}

其中出现的 快速傅里叶逆变换（IFFT） 方法FFTHelper.iFFT参考快速傅里叶变换（FFT）的C#实现及详细注释；

方法diff为差分函数，具体代码如下：

double[] diff(double[] data)
{
    double[] tmp = new double[data.Length - 1];

    for (int i = 0, length = data.Length; i < length - 1; i++)
    {
        tmp[i] = data[i + 1] - data[i];
    }

    return tmp;
}

方法hamming为汉明窗函数，参考【Matlab帮助中心】汉明窗，具体代码：

double[] hamming(int nn)
{
    double[] result = new double[nn];

    for (Int32 i = 0; i < nn; i++)
    {
        result[i] = 2 * Math.PI * i / (nn - 1);
    }

    return result.Select(o => 0.54 - 0.46 * Math.Cos(o)).ToArray();
}

数据验证

int order = 42;
double[]  freqs = { 0, 1 / 20.0, 2 / 20.0, 3 / 20.0, 4 / 20.0, 5 / 20.0, 6 / 20.0, 7 / 20.0, 8 / 20.0, 9 / 20.0, 10 / 20.0, 11 / 20.0, 12 / 20.0, 13 / 20.0, 14 / 20.0, 15 / 20.0, 16 / 20.0, 17 / 20.0, 18 / 20.0, 19 / 20.0, 1 };
double[] mags = { 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };

C#结果

0.0006145517341182119
0.000880580535935807
0.0008096195245461328
0.0001259640160925836
-0.0012933844793908147
-0.0030050738259155526
-0.0038286109689454305
-0.002307170038675055
0.0022055253858808514
0.008417746435538696
0.012858839141055355
0.01118769172861607
0.0008358319385345649
-0.016268669759934275
-0.03265882126190443
-0.03721408591907582
-0.019766562072078633
0.023447476881272418
0.08625793019711207
0.15316735202088866
0.20442751665090597
0.2236328125
0.20442751665090603
0.1531673520208887
0.0862579301971121
0.023447476881272435
-0.019766562072078615
-0.03721408591907584
-0.032658821261904454
-0.016268669759934285
0.0008358319385345515
0.011187691728616065
0.012858839141055357
0.008417746435538699
0.0022055253858808523
-0.0023071700386750596
-0.0038286109689454327
-0.0030050738259155517
-0.0012933844793908192
0.00012596401609258273
0.0008096195245461329
0.0008805805359358077
0.0006145517341182123

Matlab结果

0.00061455173411821193123899531585153
0.00088058053593580685753511305691177
0.00080961952454613279536732584062975
0.00012596401609258324146879692850831
-0.0012933844793908163932349975056013
-0.0030050738259155529950861218679847
-0.0038286109689454304871869538828832
-0.0023071700386750560937687559714959
 0.0022055253858808505570965330377931
 0.0084177464355386956429416756009232
  0.012858839141055353541553962770649
  0.011187691728616073333357761043771
0.00083583193853456266234763782918549
 -0.016268669759934274721135949448581
 -0.032658821261904440180767750234736
 -0.037214085919075837183278565589717
 -0.019766562072078629197502408487708
  0.023447476881272404092282357623844
  0.086257930197112073211762606206321
   0.15316735202088868628145235106786
   0.20442751665090599932916859415855
                         0.2236328125
   0.20442751665090599932916859415855
   0.15316735202088871403702796669677
  0.086257930197112100967338221835234
    0.0234474768812724422561988291136
 -0.019766562072078611850267648719637
  -0.03721408591907582330549075777526
 -0.032658821261904440180767750234736
 -0.016268669759934288598923757263037
0.00083583193853455182032591297414115
  0.011187691728616066394463857136543
  0.012858839141055355276277438747456
 0.0084177464355386991123886275545374
  0.002205525385880854460224353985609
-0.0023071700386750543590452799946888
-0.0038286109689454322219104298596903
 -0.003005073825915552561405252873783
-0.0012933844793908155258732595171978
 0.0001259640160925829704182538071322
0.00080961952454613279536732584062975
0.00088058053593580729121598205111354
0.00061455173411821225649964706150286