title:Calculate fir coefficients with C
keywords:c,linux,fir,dsp,octave
# 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
<script>
MathJax = {
tex: {
inlineMath: [['$', '$'], ['\\(', '\\)']]
}
};
</script>
<script type="text/javascript" id="MathJax-script" async
src="/js/tex-chtml.js">
</script>
Low pass ideal impulse response
Ideal impulse response for low-pass filter
$$2f_c\frac{sin(n \omega_c)}{n \omega_c}$$
```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
```c
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
```c
for (i=0;i<filter_N1;i++) {
w[i] = 1.0;
}
```
Flip window, against middle as curve is symmetrical against center
```c
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]$
```c
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](/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
```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);
}
```
### Highpass
```c
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
```c
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
```c
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
```c
//rectangular window
for (i=0;i<filter_N1;i++) {
w[i] = 1.0;
}
```
### Hamming
```c
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
```c
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
```c
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.
```matlab
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](http://git.main.lv/cgit.cgi/code-snippets.git/tree/fir1)
to compile get and compile code run
```bash
git clone http://git.main.lv/cgit.cgi/code-snippets.git
cd code-snippets/fir1
make
```
run program
```bash
./simple_fir
```
## Links
[main.lv/writeup/dsp_lp_filter.md](/writeup/dsp_lp_filter.md)
[http://git.main.lv/cgit.cgi/code-snippets.git/tree/fir1](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)
