Merge pull request #12 from mike919192/octave

Add octave script for generating impulse response

test_data/WriteImpulse.m

@@ -0,0 +1,97 @@
+
+%this requires signal 1.4.4 because there is a bug in the previous version where zp2sos function return wrong result!
+pkg load signal
+
+clear all
+close all
+fs = 39e3;
+fs2 = fs / 2;
+n = 1000;
+impulse = [1,zeros(1, n-1)];
+step = fs/n;
+x = 0:step:fs - step;
+f0 = 200;
+Q = 1.4;
+N = 8;
+
+%HP filter
+[z, p, k] = butter ( N, f0 / fs2, "high" );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/HPimpulse.csv", [2, fs, f0, Q, n, filtered]);
+
+%LP filter
+[z, p, k] = butter ( N, f0 / fs2 );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/LPimpulse.csv", [1, fs, f0, Q, n, filtered]);
+
+%BP filter
+[f1, f2] = findIIRCutoffFreq(fs, f0, Q);
+[z, p, k] = butter( N/2, [f1 / fs2, f2 / fs2] );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/BPimpulse.csv", [3, fs, f0, Q, n, filtered]);
+
+clear all
+close all
+fs = 39e3;
+fs2 = fs / 2;
+n = 1000;
+impulse = [1,zeros(1, n-1)];
+step = fs/n;
+x = 0:step:fs - step;
+f0 = 2000;
+Q = 0.8;
+N = 8;
+
+%HP filter
+[z, p, k] = butter ( N, f0 / fs2, "high" );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/HPimpulse2.csv", [2, fs, f0, Q, n, filtered]);
+
+%LP filter
+[z, p, k] = butter ( N, f0 / fs2 );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/LPimpulse2.csv", [1, fs, f0, Q, n, filtered]);
+
+%BP filter
+[f1, f2] = findIIRCutoffFreq(fs, f0, Q);
+[z, p, k] = butter( N/2, [f1 / fs2, f2 / fs2] );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/BPimpulse2.csv", [3, fs, f0, Q, n, filtered]);
+
+clear all
+close all
+fs = 39e3;
+fs2 = fs / 2;
+n = 1000;
+impulse = [1,zeros(1, n-1)];
+step = fs/n;
+x = 0:step:fs - step;
+f0 = 15000;
+Q = 2.0;
+N = 8;
+
+%HP filter
+[z, p, k] = butter ( N, f0 / fs2, "high" );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/HPimpulse3.csv", [2, fs, f0, Q, n, filtered]);
+
+%LP filter
+[z, p, k] = butter ( N, f0 / fs2 );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/LPimpulse3.csv", [1, fs, f0, Q, n, filtered]);
+
+%BP filter
+[f1, f2] = findIIRCutoffFreq(fs, f0, Q);
+[z, p, k] = butter( N/2, [f1 / fs2, f2 / fs2] );
+sos = zp2sos(z, p, k);
+filtered = sosfilt(sos, impulse);
+csvwrite("impulse_response/BPimpulse3.csv", [3, fs, f0, Q, n, filtered]);
+

test_data/findIIRCutoffFreq.m

@@ -0,0 +1,37 @@
+
+%finds approximate cutoff frequency using eq 8-12 pg 90 of DSP filters book
+%first find f1 and then uses the relationship with Q and f0 to find f2
+function [f1, f2] = findIIRCutoffFreq(fs, f0, Q)
+
+  theta0 = 2 * pi * f0 / fs;
+  beginTheta = 0;
+  endTheta = theta0;
+  if (theta0 > 0.01)
+    step = 0.01;
+  elseif (theta0 > 1e-4)
+    step = 1e-4;
+  else
+    step = 1e-6;
+  endif
+
+  while(step > 1e-16)
+    x = beginTheta : step : endTheta;
+    result = (sin(x) .* tan(theta0 ./ (2 .* Q))) ./ sqrt((sin(x) .* tan(theta0 ./ (2 .* Q))).^2 + (cos(x) - cos(theta0)).^2) - 1/sqrt(2);
+    %indicate zero crossing with sign change
+    signResult = sign(result);
+
+    %technically possible we could have a step land right on zero so check for it
+    if (isempty(find(signResult == 0)) == false)
+      beginTheta = (find(signResult == 0, 1) - 1) * step + beginTheta;
+      break;
+    endif
+
+    endTheta = (find(signResult == 1, 1) - 1) * step + beginTheta;
+    beginTheta = endTheta - step;
+    step = step * 0.01;
+  endwhile
+
+  f1 = beginTheta * fs / (2 * pi);
+  f2 = f0 / Q + f1;
+
+endfunction