fix bug in calculation of power in frequency bands

calibration/src/capsule_polar_response.m

@@ -55,7 +55,7 @@ function [CAP] = capsule_polar_response(measurements, fmin, fmax, N, Npre)
     equalized = fft(capsule(range), N) ./ ref_fft;
     
     %% return total power of signal in frequency range
-    CAP(i,:) = sqrt(sum(equalized(imin:imax).*conj(equalized(imin:imax)))/size(equalized(imin:imax),1));
+    CAP(i,:) = sqrt(sum(equalized(imin:imax).*conj(equalized(imin:imax)))/(imax - imin));
     i = i + 1;
   end
 

calibration/src/sphear_af_power.m

@@ -57,7 +57,7 @@ function [R Fs] = sphear_af_power(measurements, fmin, fmax, Nfft)
       equalized = fft(capsule, Nfft);
       
       %% return total power of signal in frequency range
-      R(i, c) = sqrt(sum(equalized(imin:imax).*conj(equalized(imin:imax)))/size(equalized(imin:imax),1));
+      R(i, c) = sqrt(sum(equalized(imin:imax).*conj(equalized(imin:imax)))/(imax - imin));
 
     end
     i = i + 1;

calibration/src/sphear_bf_power.m

@@ -47,10 +47,10 @@ function [B] = sphear_bf_power(fmin, fmax, N, Fs, plotz)
       cfft = fft(bsignal(:,1) + 1.0 * bsignal(:,2), N);
     end
     %% return total power of signal in frequency range
-    B(m, :) = sqrt(sum(bfft(imin:imax,:).*conj(bfft(imin:imax,:)))/size(bfft(imin:imax,:),1));
+    B(m, :) = sqrt(sum(bfft(imin:imax,:).*conj(bfft(imin:imax,:)))/(imax - imin));
     %% make a cardiod if we don't want Z
     if (plotz == 0)
-      B(m, 4) = sqrt(sum(cfft(imin:imax).*conj(cfft(imin:imax)))/size(cfft(imin:imax),1)) / 2.0;
+      B(m, 4) = sqrt(sum(cfft(imin:imax).*conj(cfft(imin:imax)))/(imax - imin)) / 2.0;
     end
   end
   

calibration/src/sphear_bf_power_fons.m

@@ -67,10 +67,10 @@ function [B] = sphear_bf_power_fons(fmin, fmax, N, Fs, plotz)
     end
     %% return total power of signal in frequency range
     scaler = 8;
-    B(m, :) = sqrt(sum(bfft(imin:imax,:).*conj(bfft(imin:imax,:)))/size(bfft(imin:imax,:),1)) / scaler;
+    B(m, :) = sqrt(sum(bfft(imin:imax,:).*conj(bfft(imin:imax,:)))/(imax - imin)) / scaler;
     %% make a cardiod if we don't want Z
     if (plotz == 0)
-      B(m, 4) = sqrt(sum(cfft(imin:imax).*conj(cfft(imin:imax)))/size(cfft(imin:imax),1)) / (2 * scaler);
+      B(m, 4) = sqrt(sum(cfft(imin:imax).*conj(cfft(imin:imax)))/(imax - imin)) / (2 * scaler);
     end
   end
   

calibration/src/sphear_capsule_power.m

@@ -33,7 +33,9 @@
 function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, N, Npre)
   pkg load signal;
 
-  onames = cellstr(['00'; '90']);
+  capname = "EM200";
+  
+  onames = cellstr(['001'; '002']);
 
   %% all submitted measurements
   i = 1;
@@ -41,9 +43,9 @@ function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, N, Npre)
     %% orientation
     for o = [1, 2]
 
-      printf("%s\n", ["../data/2016.07.26/export/capsule-" sprintf("%02.0f", m) "-" char(onames(o)) "-1.wav"]);
+      printf("%s\n", ["../data/current/" capname "/export/capsule-" sprintf("%02.0f", m) "-" char(onames(o)) "-1.wav"]);
       
-      [capsule, Fs] = wavread(["../data/2016.07.26/export/capsule-" sprintf("%02.0f", m) "-"  char(onames(o)) "-1.wav"]);
+      [capsule, Fs] = wavread(["../data/current/" capname "/export/capsule-" sprintf("%02.0f", m) "-"  char(onames(o)) "-1.wav"]);
 
       %% indexes in array for frequency bins
       imin = round((N/2) * fmin / (Fs/2));
@@ -58,7 +60,7 @@ function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, N, Npre)
       fdom = fft(capsule(range), N);
       
       %% return total power in dB of signal in frequency range
-      R(i, o) = 20*log10(sqrt(sum(fdom(imin:imax).*conj(fdom(imin:imax)))/size(fdom(imin:imax),1)));
+      R(i, o) = 20*log10(sqrt(sum(fdom(imin:imax).*conj(fdom(imin:imax)))/(imax - imin)));
     end
     %% add 0 to 90 degree difference in dB (front to side ratio)
     R(i, 3) = R(i, 1) - R(i, 2);
@@ -74,7 +76,9 @@ function [RS Fs] = sphear_capsule_power(measurements, fmin, fmax, N, Npre)
   R(:,1) = R(:,1) - avg;
   R(:,2) = R(:,2) - avg;
   %% sort the results by directivity
-  [RS, I] = sort(R'(3,:));
+  %% [RS, I] = sort(R'(3,:));
+  %% sort the results by 0 degree gain
+  [RS, I] = sort(R'(1,:));
   RS = R'(:,I)';
   %% calculate the derivative of directivity
   i = 1;

calibration/src/sphear_fir_extract.m

@@ -1,3 +1,4 @@
+%% -*- mode: octave -*-
 %%
 %% *SpHEAR Project
 %%