From 068d08bf44aafc1b36fec7302130f1a5b2fde08c Mon Sep 17 00:00:00 2001
From: ZoRo <dos21h@gmail.com>
Date: Fri, 9 Jul 2021 09:01:22 +0100
Subject: FIR filter C snippet

---
 fir1/Makefile     |  2 ++
 fir1/simple_fir.c | 97 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 fir1/verify_fir.m | 85 ++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 184 insertions(+)
 create mode 100644 fir1/Makefile
 create mode 100644 fir1/simple_fir.c
 create mode 100644 fir1/verify_fir.m

diff --git a/fir1/Makefile b/fir1/Makefile
new file mode 100644
index 0000000..70c416d
--- /dev/null
+++ b/fir1/Makefile
@@ -0,0 +1,2 @@
+make:
+	gcc simple_fir.c -o simple_fir -lm
\ No newline at end of file
diff --git a/fir1/simple_fir.c b/fir1/simple_fir.c
new file mode 100644
index 0000000..8f48569
--- /dev/null
+++ b/fir1/simple_fir.c
@@ -0,0 +1,97 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <math.h>
+
+double sinc(double x) {
+    if (x>-1.0E-5 && x < 1.0E-5) return (1.0);
+    return sin(x)/x;
+}
+
+int main() {
+
+	const int buf_N=200;
+	//int N1 = N/2;
+    
+    int i=0;
+    double h[buf_N],w[buf_N];
+    int filter_N=20;
+    int even = filter_N ? 0: 1;
+    int filter_N1=filter_N/2;
+    int filter_FS_hz=10000;
+    int filter_FC_hz=2000;
+    double omega_c=1.0*filter_FC_hz/filter_FS_hz;
+    double arg;
+    double dN;
+
+    //nullify all coefficients
+    for (i=0;i<filter_N;i++) {
+        h[i] = 0.0;
+        w[i] = 0.0;
+    }
+
+    //sinc the
+    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);
+    }
+
+
+    if (even) filter_N1+=1;
+
+    #if 1
+    //rectangular window
+    
+    for (i=0;i<filter_N1;i++) {
+        w[i] = 1.0;
+    }
+    #endif
+
+    #if 0
+    //hamming window
+    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);
+
+    }
+    #endif
+
+    #if 0
+    //hanning window
+    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);
+
+    }
+    #endif
+
+    #if 0
+    //blackman window
+    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));
+
+    }
+    #endif
+
+    //flip window
+    for (i=0;i<filter_N1;i++) {
+        w[filter_N-i-1] = w[i];
+    }
+
+    //convolution
+    for (i=0;i<filter_N;i++) {
+        h[i] = h[i]*w[i];
+    }
+
+    
+    //print coefficients
+    printf("h=[\n");
+    for (i=0;i<filter_N;i++) {
+        printf("%f\n",h[i]);
+    }
+    printf("];\n");
+
+}
diff --git a/fir1/verify_fir.m b/fir1/verify_fir.m
new file mode 100644
index 0000000..ebdcdf3
--- /dev/null
+++ b/fir1/verify_fir.m
@@ -0,0 +1,85 @@
+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.000000
+-0.000310
+-0.001913
+-0.004532
+-0.004058
+0.008358
+0.041693
+0.095752
+0.155110
+0.194532
+0.194532
+0.155110
+0.095752
+0.041693
+0.008358
+-0.004058
+-0.004532
+-0.001913
+-0.000310
+0.000000
+];
+ 
+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");
+
-- 
cgit v1.2.3