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

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<filter_N1;i++) {
w[i] = 1.0;
}
``` |

Flip window, against middle as curve is symmetrical against center

1 2 3 | ```
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]$

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

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

1 2 3 4 5 6 7 8 9 | ```
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

1 2 3 4 5 6 7 8 | ```
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

1 2 3 4 5 6 7 8 | ```
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

1 2 3 4 | ```
//rectangular window
for (i=0;i<filter_N1;i++) {
w[i] = 1.0;
}
``` |

### Hamming

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

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

1 2 3 4 5 6 7 | ```
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.

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

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