1 files changed, 85 insertions, 0 deletions
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");
+