diff options
Diffstat (limited to 'md')
| -rw-r--r-- | md/writeup/calculate_fir_coefficients_with_c.md | 238 |
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 |