Calculate fir filter with c

1 files changed, 238 insertions, 0 deletions

diff --git a/md/writeup/calculate_fir_coefficients_with_c.md b/md/writeup/calculate_fir_coefficients_with_c.md
new file mode 100644
index 0000000..c8c9552
--- /dev/null
+++ b/md/writeup/calculate_fir_coefficients_with_c.md
@@ -0,0 +1,238 @@
+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. There is few more but those for later.
+
+## 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
+
+Impulse response of the 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);
+}
+```
+
+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]$ agains 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)
+
+## 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.
+
+### 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://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