[FIR window function filter] C language implementation and comparison with Matlab window+fir1

Reference documentation

Main reference sources Digital signal processing Computer-based methods 4th edition (Chinese version) and Matlab window function [Matlab Help Center] windowing method two parts

Category

Filter type classification

  • Lowpass
  • Qualcomm
  • Bandpass
  • band resistance
    Different filter types correspond to differentideal unit impulse functions hd(n)

Window function classification

  • rectangular window
  • triangular window
  • Hanning window
  • Hamming window
  • blackman window
  • flat top window
  • Gaussian window
  • kaiser window
  • And some other uncommon windows, such as Parzen, Bohman, Chebyshev windows, etc.

C language implementation

Ideal unit impulse response function

It is implemented based on the introduction of the ideal unit impulse response function (as shown below) in Digital Signal Processing Computer-Based Methods 4th Edition (Chinese Version).
Ideal unit shock response function

  • Lowpass
// wc is the cut-off frequency, N is the window length
double *IdealLp(double wc, int N) //low pass
{<!-- -->
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));
    for (int n = 0; n < N; n + + )
    {<!-- -->
        if(n==M)
        {<!-- -->
            hd[n] = wc / PI;
        }
        else
        {<!-- -->
            hd[n] = sin(wc * (n - M)) / (PI * (n - M));
        }
    }

    return hd;
}
  • Qualcomm
// wc is the cut-off frequency, N is the window length
double *IdealHp(double wc, int N) //Qualcomm
{<!-- -->
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));
    for (int n = 0; n < N; n + + )
    {<!-- -->
        if(n==M)
        {<!-- -->
            hd[n] = 1.0 - wc / PI;
        }
        else
        {<!-- -->
            hd[n] = -sin(wc * (n - M)) / (PI * (n - M));
        }
    }

    return hd;
}
  • Bandpass
// wcl is the low-frequency cut-off frequency, wch is the high-frequency cut-off frequency, and N is the window length
double *IdealBp(double wcl, double wch, int N) // bandpass
{<!-- -->
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));

    hd = IdealLp(wcl, N);
    hd = subtract(IdealLp(wch, N), hd, N);

    return hd;
}

The bandpass ideal unit impulse response function can be obtained by subtracting the lowpass ideal unit impulse response function at the low frequency and high frequency cutoff frequencies.

Among them, subtract is the subtraction function in the encapsulated double array operation library, and IdealLp is the low-pass ideal unit impulse response function.

  • band resistance
// wcl is the low-frequency cut-off frequency, wch is the high-frequency cut-off frequency, and N is the window length
double *IdealBs(double wcl, double wch, int N) // band stop
{<!-- -->
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));
    for (int n = 0; n < N; n + + )
    {<!-- -->
        if(n==M)
        {<!-- -->
            hd[n] = 1.0 - (wch - wcl) / PI;
        }
        else
        {<!-- -->
            hd[n] = (sin(wcl * (n - M)) - sin(wch * (n - M))) / (PI * (n - M));
        }
    }

    return hd;
}

At this point, the ideal unit impulse response functions of low pass, high pass, band pass and band stop have all been realized.

Window function

Rectangular window
//Array of all 1’s
double *RectangleWindow(int n)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = 1.0;
    }

    return result;
}
Triangle window
double *TriangWindow(int n)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    if (n % 2 == 0) // even
    {<!-- -->
        for (size_t i = 1; i <= n; i + + )
        {<!-- -->
            if (i >= 1 & amp; & amp; i <= n / 2)
            {<!-- -->
                result[i - 1] = (2.0 * i - 1.0) / n;
            }
            else
            {<!-- -->
                result[i - 1] = 2.0 - (2.0 * i - 1.0) / n;
            }
        }
    }
    else if (n % 2 == 1) // odd
    {<!-- -->
        for (size_t i = 1; i <= n; i + + )
        {<!-- -->
            if (i >= 1 & amp; & amp; i <= (n + 1) / 2)
            {<!-- -->
                result[i - 1] = 2.0 * i / (n + 1.0);
            }
            else
            {<!-- -->
                result[i - 1] = 2.0 - 2.0 * i / (n + 1.0);
            }
        }
    }

    return result;
}
Hanning Window
double *HannWindow(int n, double *w)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = 0.5 - 0.5 * cos(w[i]);
    }

    return result;
}

The parameter w is calculated in advance as a public parameter

 double *w = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
            w[i] = 2 * PI * i / (n - 1);
    }
Hamming window
double *HammingWindow(int n, double *w)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = 0.54 - 0.46 * cos(w[i]);
    }

    return result;
}
Brackman window
double *HammingWindow(int n, double *w)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = 0.54 - 0.46 * cos(w[i]);
    }

    return result;
}
Flat top window
double *FlattopWindow(int n, double *w)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    double fta0 = 0.215578995;
    double fta1 = 0.41663158;
    double fta2 = 0.277263158;
    double fta3 = 0.083578947;
    double fta4 = 0.006947368;
    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = fta0 -
                    fta1 * cos(w[i]) +
                    fta2 * cos(2 * w[i]) -
                    fta3 * cos(3 * w[i]) +
                    fta4 * cos(4 * w[i]);
    }

    return result;
}
Gaussian window
double *GaussianWindow(int n)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    double alpha = 2.5;
    double width = (n - 1) / 2.0;

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        double length = i - width;
        result[i] = width == 0 ? 0 : exp(-0.5 * pow(alpha * lenth / width, 2));
    }

    return result;
}

Among them, alpha can be added to the parameters. Here, the fixed value 2.5 is selected. The alpha selection principle is:

  • alpha < 1: If you need a narrower main lobe and a wider side lobe, choose a smaller alpha. This may be useful for applications requiring better frequency resolution, but will result in slower sidelobe roll-off.
  • 1 <= alpha <= 2: This is a balancing option and usually the most commonly used. It provides a relatively wide main lobe and moderately fast side lobe descent.
  • alpha > 2: If you need a very wide main lobe and a very fast side lobe drop, you can choose a larger alpha. This may be useful for certain applications, but may reduce frequency resolution.
Kaiser Window
double *KaiserWindow(int n)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    double beta = 7;
    // β satisfies the conditions, α>50, β=0.1102 (α-8.7);
    // 50>=α>=21,β=0.5842(α?21)^0.4 + 0.07886(α?21);
    // α<21,β=0;

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        // double x = beta * sqrt(1 - pow(1 - 2.0 * i / (n - 1), 2)); // old
        double x = beta * sqrt(1 - pow((i-(n-1)/2.0)/((n-1)/2.0), 2)); //new
        result[i] = I0function(x) / I0function(beta);
    }

    return result;
}

Among them, beta can be passed in as a parameter. It is generally calculated based on the filter performance index. Here, the fixed value 7* is selected;
The method I0function is the zero-order modified Bessel function of the first kind, and its specific implementation is:

double I0function(double x)
{<!-- -->
    double i0 = 0;
    for (size_t i = 0; i < 25; i + + )
    {<!-- -->
        i0 + = pow(x, 2 * i) / pow(4, i) / pow(Fact(i), 2);
    }

    return i0;
}

double Fact(int n) // factorial
{<!-- -->
    if(n==0)
    {<!-- -->
        return 1.0;
    }

    double y = (double)n;
    for (double m = n - 1; m > 0; m--)
    {<!-- -->
        y *= m;
    }

    return y;
}
Other window functions

Bartley Window

double *BartlettWindow(int n)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        if (i <= (n - 1) / 2)
        {<!-- -->
            result[i] = 2.0 * i / (n - 1.0);
        }
        else
        {<!-- -->
            result[i] = 2.0 - 2.0 * i / (n - 1.0);
        }
    }

    return result;
}
Bartley-Hanning Window
double *BartlettHanningWindow(int n)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = 0.62 - 0.48 * fabs((1.0 * i / (n - 1.0)) - 0.5) + 0.38 * cos(2.0 * PI * ((1.0 * i / (n - 1.0)) - 0.5)) ;
    }

    return result;
}
Blackman-Harris window
//type=1 means symmetry, =2 means periodicity
double *BlackmanHarrisWindow(int n, int type)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));
    double a0 = 0.35875;
    double a1 = 0.48829;
    double a2 = 0.14128;
    double a3 = 0.01168;

    double m = type == 0 ? (n - 1.0) : n;
    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = a0 -
                    a1 * cos(2.0 * PI * i / m) +
                    a2 * cos(2.0 * 2.0 * PI * i / m) -
                    a3 * cos(3.0 * 2.0 * PI * i / m);
    }

    return result;
}

Does it look familiar? It’s similar to the formula for Flat-roofed windows, but the coefficients are different.

Parzen
double *ParzenWindow(int n)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));

    int index = 0;
    for (int i = -(n - 1) / 2; i <= (n - 1) / 2; i + + )
    {<!-- -->
        if (abs(i) >= 0 & amp; & amp; abs(i) <= (n - 1) / 4)
        {<!-- -->
            result[index + + ] = 1 -
                              6.0 * pow(fabs(i) / (n / 2.0), 2) +
                              6.0 * pow(fabs(i) / (n / 2.0), 3);
        }
        else
        {<!-- -->
            result[index + + ] = 2 * pow((1.0 - fabs(i) / (n / 2.0)), 3);
        }
    }

    return result;
}
Nutall’s Blackman-Harris Window
double *NuttallWindow(int n, int type)
{<!-- -->
    double *result = (double *)malloc(n * sizeof(double));
    double a0 = 0.3635819;
    double a1 = 0.4891775;
    double a2 = 0.1365995;
    double a3 = 0.0106411;

    double m = type == 0 ? (n - 1.0) : n;
    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        result[i] = a0 -
                    a1 * cos(2.0 * PI * i / m) +
                    a2 * cos(2.0 * 2.0 * PI * i / m) -
                    a3 * cos(3.0 * 2.0 * PI * i / m);
    }

    return result;
}

Another window function similar to Flat-top window. Again, their coefficients are different.

Chebyshev window
double *ChebyshevWindow(int n, double ripple) // Chebyshev window function, ripple is the weakening coefficient, the default is 100dB
{<!-- -->

}

Leave it empty until it is finished, add it later

Filter coefficient

Through the above code, we can obtain the ideal unit shock response function coefficient hd(n) and window function coefficient w(n) of the corresponding length, and finally multiply them correspondingly, and then do Normalized.

double *ApplyWindow(double *hd, int windowType, int n)
{<!-- -->
    double *w = GetWindowCoefficient(windowType, n);
    double *h = multiply(hd, w, n);//corresponding to multiplication
    double sum_h = 0.0;

    for (size_t i = 0; i < n; i + + )
    {<!-- -->
        sum_h + = h[i];
    }

               
    h = multiply_num(h, 1.0 / sum_h, n);

    return h;
}

Test

  • matlab code implementation
function CreateFirByWindow()
disp('Hello, Test Beginning!')

N = input('Please enter order!\
'); % Order
windowType=input('Please enter window type!(Rectangle 0,Hann 1,Hamming 2,Blackman 3,Flattop 4,Kaiser 5,Gaussian 6,Triang 7,Bartlett 8,Bartlett-Hann 9,Blackman-Harris(Symmetric) 10 ,Blackman-Harris(Periodic) 11,Parzen 12,Nuttall(Symmetric) 13,Nuttall(Periodic) 14)\
'); %Window Type
filterResponseType=input('Please enter filter response type!(LowPass 0, HighPass 1,BandPass 2,BandStop 3)\
'); %Filter Response Type
wc1 = input('Please enter wc1!(0<wc1<1)\
');
% wc1 = 1.0e-323;
wc2=0.6;
if filterResponseType==2||filterResponseType==3
    wc2 = input('Please enter wc2!(0<wc2<1)\
');
else
    wc2=0.6;
end

win = CreateWindow(windowType,N);

type='low';
switch filterResponseType
    case 0
        type='low';
    case 1
        type='high';
    case 2
        type='bandpass';
    case 3
        type='stop';
end
%
if filterResponseType<2
    w=fir1(N, wc1, type, win, "scale");
else
    w=fir1(N,[wc1,wc2],type,win,"scale");
end
w=vpa(w,20);
disp(w.'); %Print the coefficient before normalization

function win = CreateWindow(windowType,order)

win=zeros(1,order + 1);

switch windowType
    case {<!-- -->0}
        win=rectwin(order + 1);
    case {<!-- -->1}
        win=hann(order + 1);
    case {<!-- -->2}
        win=hamming(order + 1);
    case {<!-- -->3}
        win=blackman(order + 1);
    case {<!-- -->4}
        win=flattopwin(order + 1);
    case {<!-- -->5}
        beta=7; %kaiser parameter defaults to 7
        win=kaiser(order + 1,beta);
    case {<!-- -->6}
        alpha=2.5;%gaussian parameter defaults to 2.5
        win=gausswin(order + 1,alpha);
    case {<!-- -->7}
        win=triang(order + 1);
    case {<!-- -->8}
        win=bartlett(order + 1);
    case {<!-- -->9}
        win=barthannwin(order + 1);
    case {<!-- -->10}
        win=blackmanharris(order + 1,"symmetric");
    case {<!-- -->11}
        win=blackmanharris(order + 1,"periodic");
    case {<!-- -->12}
        win=parzenwin(order + 1);
    case {<!-- -->13}
        win=nuttallwin(order + 1,"symmetric");
    case {<!-- -->14}
        win =nuttallwin(order + 1,"periodic");
    case {<!-- -->15}
        win=chebwin(order + 1,100);
end

Tested by different orders (single digits, tens, hundreds, thousands), parity (the order is divided into odd and even numbers) and limit values

Order and Parity test results:

There are some differences after 10 decimal places, but the overall results are basically the same

Limit valueTest results

Matlab:
Lower limit frequency: 1.0e-323, upper limit frequency 0.9999999999999999 (16 decimal places)
C:
The lower limit frequency is: 1.0e-309, the upper limit frequency is greater than 0.9999999999999999 (greater than 16 decimal places)

Basically passed, test completed