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

```
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](/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)