# 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}$$

 1 2 3 4 for (i=0;i

Sinc function implementation

 1 2 3 4 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

 1 2 3 for (i=0;i

Flip window, against middle as curve is symmetrical against center

 1 2 3 for (i=0;i

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

 1 2 3 for (i=0;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

 1 2 3 4 for (i=0;i

### Highpass

 1 2 3 4 5 6 7 8 9 for (i=0;i

### Bandpass

 1 2 3 4 5 6 7 8 for (i=0;i

### Bandstop

 1 2 3 4 5 6 7 8 for (i=0;i

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

 1 2 3 4 //rectangular window for (i=0;i

### Hamming

 1 2 3 4 5 for (i=0;i

### Hanning

 1 2 3 4 5 for (i=0;i

### Blackman

 1 2 3 4 5 6 7 for (i=0;i

## 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.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 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

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

run program

 1 ./simple_fir