Calculate fir coefficients with C

Intro

FIR filter's is probably ones of fascinating parts when you want to do with DSP, its kind of too simple from first glance and they works. As before I have learned about FIR and IIR and was able to use them just from tools that take as input parameters and give as output coefficients. Now time came just to calculate FIR from ground up. Also C is perfect to do stuff without abstraction levels hiding details, and can be ported to other languages and platforms. Lets calculate filter coefficients with windows method. Theory is not explained, there is a lot of playlists on youtube and lecture notes online to check on that.

Implementation

Low pass ideal impulse response

Ideal impulse response for low-pass filter

$$2f_c\frac{sin(n \omega_c)}{n \omega_c}$$

for (i=0;i<filter_N;i++) {
        arg = (double)i-(double)(filter_N-1)/2.0;
        h[i] = omega_c * sinc(omega_c*arg*M_PI);
}

Sinc function implementation

double sinc(double x) {
    if (x>-1.0E-5 && x < 1.0E-5) return (1.0);
    return sin(x)/x;
}

Rectangular window for low pass filter

for (i=0;i<filter_N1;i++) {
        w[i] = 1.0;
}

Flip window, against middle as curve is symmetrical against center

for (i=0;i<filter_N1;i++) {
        w[filter_N-i-1] = w[i];
}

Convolution of ideal filter response $h[i]$ against window function $w[i]$

for (i=0;i<filter_N;i++) {
        h[i] = h[i]*w[i];
}

How to use calculated coefficients /writeup/dsp_lp_filter.md#toc-6

Ideal impulse response

Filter type Ideal impulse response
Lowpass $2f_c\frac{sin(n \omega_c)}{n \omega_c}$
Highpass $-2f_c\frac{sin(n \omega_c)}{n \omega_c}$
Bandpass $2f_c\frac{sin(n \omega_2)}{n \omega_2}-2f_c\frac{sin(n \omega_1)}{n \omega_1}$
Bandstop $2f_c\frac{sin(n \omega_1)}{n \omega_1}-2f_c\frac{sin(n \omega_2)}{n \omega_2}$

Lowpass

for (i=0;i<filter_N;i++) {
        arg = (double)i-(double)(filter_N-1)/2.0;
        h[i] = omega_c * sinc(omega_c*arg*M_PI);
}

Highpass

for (i=0;i<filter_N;i++) {
    arg = (double)i-(double)(filter_N-1)/2.0;
    if (arg == 0.0) {
        h[i] = 0.0;
    } else {

        h[i] = cos(omega_c*arg*M_PI)/M_PI/arg + cos(arg*M_PI);
    }
}

Bandpass

for (i=0;i<filter_N;i++) {
        arg = (double)i-(double)(filter_N-1)/2.0;
        if (arg==0.0) {
            h[i] = 0.0;
        } else {
            h[i] = (cos(omega_c1*arg*M_PI) - cos(omega_c2*arg*M_PI))/M_PI/arg;
        }
    }

Bandstop

    for (i=0;i<filter_N;i++) {
        arg = (double)i-(double)(filter_N-1)/2.0;
        if (arg==0.0) {
            h[i] = 0.0;
        } else {
            h[i] = sinc(arg*M_PI)  - omega_c2*sinc(omega_c2*arg*M_PI) - omega_c1*sinc(omega_c1*arg*M_PI);
        }
    }

Windowing methods

Here is most common windows for window methods that can give you better results then naive rectangular windows filter. There is better ways how to calculate filters, but that for laters.

Name Main lobe Stop band attenuation Window function
Rectangluar 13dB 21dB 1
Hanning 31dB 44dB $0.5+0.5cos(\frac{2\pi n}{N})$
Hamming 41dB 53dB $0.54 + 0.46 cos\frac{2\pi n}{N}$
Blackman 57dB 74dB $0.42+0.5cos(\frac{2\pi n}{N-1}) + 0.08 cos(\frac{4\pi n }{N-1})$

Rectangular

//rectangular window
for (i=0;i<filter_N1;i++) {
        w[i] = 1.0;
}

Hamming

for (i=0;i<filter_N1;i++) {
        arg = (double)i-(double)(filter_N-1)/2.0;
        w[i] = 0.54+0.46*cos(2*M_PI*arg/filter_N);

}

Hanning

for (i=0;i<filter_N1;i++) {
        arg = (double)i-(double)(filter_N-1)/2.0;
        w[i] = 0.5+(1-0.5)*cos(2*M_PI*arg/filter_N);

}

Blackman

for (i=0;i<filter_N1;i++) {
        arg = (double)i-(double)(filter_N-1)/2.0;
        w[i] = 0.42
        +0.50*cos(2*M_PI*arg/(filter_N-1))
        +0.08*cos(4*M_PI*arg/(filter_N-1));

}

Octave verification code

Try to run this code if there is complain that fir1 function is missing, then need to install signal and control packages for Octave. There is may need to install Fortran compiler if its not yet installed on system.

from Octave command line

pkg install -forge control
pkg install -forge signal

After FIR filters calculated then best approach it to check filter response diagram. Replace in octave source code values of fir_coef with newly calculated values, and check if filter response is as expected. If not then there is some kind of error.

pkg load signal
Fs = 10000;
t = 0:1/Fs:1-1/Fs;

x = cos(2*pi*1000*t) + cos(2*pi*1500*t) + cos(2*pi*2000*t) \
  + cos(2*pi*4000*t) + cos(2*pi*5000*t) + cos(2*pi*6000*t);

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

figure("name","Figure 1: Periodogram Using FFT")
plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
%filter
%sampling frequency
Fs = length(x);
Fnyq = Fs/2;
fc = 2000;
hc = fir1(41,fc/Fs)
f_filt = filter( hc, 1, x );


N = length(f_filt);
xdft = fft(f_filt);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

figure("name","Figure 2: Octave fir1 filtering");
plot( freq,10*log10(psdx) );
title("filtered fft");

figure("name","Figure 3: Octave fir1 filter response diagram");
freqz(hc);
title("coef freqz");

%draw my own filter results
fir_coef =[
-0.010354
-0.030296
-0.042441
-0.039618
-0.017884
0.021858
0.073577
0.127324
0.171679
0.196726
0.196726
0.171679
0.127324
0.073577
0.021858
-0.017884
-0.039618
-0.042441
-0.030296
-0.010354
];

my_hw = fir_coef;
figure("name","Figure 4: Verify filter using custom coefficients");
my_filt = filter(my_hw,1,x);
N = length(my_filt);
xdft = fft(my_filt);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;
plot( freq,10*log10(psdx) );
title("my filtered coeficients fft");

figure("name","Figure 5: Show custom filter response diagram");
freqz(my_hw);
title("my coef freqz");

Source code

Snippet code is located at http://git.main.lv/cgit.cgi/code-snippets.git/tree/fir1
to compile get and compile code run

git clone http://git.main.lv/cgit.cgi/code-snippets.git
cd code-snippets/fir1
make

run program

./simple_fir

Links

main.lv/writeup/dsp_lp_filter.md
http://git.main.lv/cgit.cgi/code-snippets.git/tree/fir1
https://en.wikipedia.org/wiki/Sinc_function
https://ccrma.stanford.edu/~jos/st/FFT_Simple_Sinusoid.html
https://ccrma.stanford.edu/~jos/st/Example_Applications_DFT.html
https://ccrma.stanford.edu/~jos/st/Use_Blackman_Window.html
http://www.iowahills.com/A7ExampleCodePage.html
Digital Signal Processing: A Practical Approach by (Emmanuel Ifeachor, Barrie Jervis)