Pārlūkot izejas kodu

change the way the filter is designed, get a spline to interpolate a smooth frequency response curve, plot some graphics to see what is going on

Fernando Pablo Lopez-Lezcano 10 gadi atpakaļ
vecāks
revīzija
fdaeb3574a
1 mainītis faili ar 23 papildinājumiem un 3 dzēšanām
  1. 23 3
      sphear_fir_design.m

+ 23 - 3
sphear_fir_design.m

@@ -51,6 +51,7 @@ function [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, XYZ)
 	if (F(d) <= Fl)
 	  %% below the limit frequency we use the A2B coefficients
 	  vlast = A2BF(d, c, b);
+	  dlast = d;
 	  if (exist("W") == 1)
 	    wlast = W(d);
 	  end
@@ -74,10 +75,29 @@ function [FIR] = sphear_fir_design(F, A2BF, n, Fs, Fl, W, XYZ)
 	  end
 	end
       end
-      %% full response
-      r = [ A2BF(1,c,b) A2BF(:,c,b)' A2BF(end,c,b)]';
+      
+      %% interpolate last point at nyquist from the previous two points
+      r = [ A2BF(1,c,b) A2BF(:,c,b)' (A2BF(end,c,b) + (A2BF(end,c,b) - A2BF(end-1,c,b)))]';
+
+      %% get a spline fitted over the points to smooth the response
+      %% endpoints have null first derivative
+      fi = (0:200/(Fs/2):1)';
+      ri = spline(f, [0 r' 0], fi);
+
+      if (b == 1)
+	semilogx(f, r,'*g');
+      else
+	semilogx(f, r,'*b');
+      end
+      hold on;
+      if (b == 1)
+	semilogx(fi,ri,'g');
+      else
+	semilogx(fi,ri,'r');
+      end
+      
       %% design the filter
-      FIR(b, c, :) = fir2(n, f, r);
+      FIR(b, c, :) = fir2(n, fi, ri);
     end
   end