diff options
| -rw-r--r-- | fir1/Makefile | 2 | ||||
| -rw-r--r-- | fir1/simple_fir.c | 97 | ||||
| -rw-r--r-- | fir1/verify_fir.m | 85 |
3 files changed, 184 insertions, 0 deletions
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"); + |