Build A Weighting Table

README.md

Install

pkg install signal pkg load signal

Output

aweightingtableR.csv aweightingtableI.csv

build_Aweighting_table.matlab

% Build A-weighting filter gain table used to filter signal
% 
% copyright (c) Davide Gironi, 2013
% copyright (c) Pral2a, 2016 (Ported to GNU Octave)
%
% Released under GPLv3.
% Please refer to LICENSE file for licensing information.

%clear matlab workspace
clear;


function aspec(B,A,Fs); 
% ASPEC  Plots a filter characteristics vs. A-weighting specifications. 
%    ASPEC(B,A,Fs) plots the attenuation of the filter defined by B and A 
%    with sampling frequency Fs against the specifications for class 1 
%    and class 2 A-weighting filters. 
%
%    See also ADSGN, CDSGN, CSPEC. 

% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
%         couvreur@thor.fpms.ac.be
% Last modification: Aug. 26, 1997, 10:00am.

% References: 
%    [1] IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996. 


if (nargin ~= 3) 
  error('Invalide number of arguments.');
end

N = 512; 
W = 2*pi*logspace(log10(10),log10(Fs/2),N);
H = freqz(B,A,W/Fs);
fref = [10, 12.5 16 20, 25 31.5 40, 50 63 80, 100 125 160, 200 250 315, ...
    400 500 630, 800 1000 1250, 1600 2000 2500, 3150 4000 5000, ...
  6300 8000 10000, 12500 16000 20000 ];
Agoal = [-70.4, -63.4 -56.7 -50.5, -44.7 -39.4 -34.6, -30.2 -26.2 -22.5, ...
  -19.1 -16.1 -13.4, -10.9 -8.6 -6.6, -4.8 -3.2 -1.9, -.8 0 .6, ...
  1 1.2 1.3, 1.2, 1 .5, -.1 -1.1 -2.5, -4.3 -6.6 -9.3 ];
up1 = [ 3, .25 2 2, 2 1.5 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 0.7 1, ...
  1 1 1, 1 1 1.5, 1.5 1.5 2, 2 2.5 3 ];
low1 = [ inf, inf 4 2, 1.5 1.5 1, 1 1 1, 1 1 1, 1 1 1, 1 1 1, 1 0.7 1, ...
  1 1 1, 1 1 1.5, 2 2.5 3, 5 8 inf ];
up2 = [ 5, 5 5 3, 3 3 2, 2 2 2, 1.5 1.5 1.5, 1.5 1.5 1.5, 1.5 1.5 1.5, ...
  1.5 1 1.5, 2 2 2.5, 2.5 3 3.5, 4.5 5 5, 5 5 5 ];
low2 = [ inf, inf inf 3, 3 3 2, 2 2 2, 1.5 1.5 1.5, 1.5 1.5 1.5, ...
    1.5 1.5 1.5, 1.5 1 1.5, 2,2 2.5, 2.5 3 3.5, 4.5 5 inf, inf inf inf ];
clf;
figure(1)
semilogx(W/2/pi,20*log10(abs(H)),'-',fref,Agoal+up1,'--',fref,Agoal+up2,':'); 
hold on 
semilogx(W/2/pi,20*log10(abs(H)),'-',fref,Agoal-low1,'--',fref,Agoal-low2,':'); 
hold off
xlabel('Frequency [Hz]'); ylabel('Gain [dB]');
title(['A-weighting, Fs = ',int2str(Fs),' Hz']);
axis([fref(1)  max(fref(34),Fs/2) -75 5]);
legend('Filter','Class 1','Class 2',0);
endfunction

function [B,A] = adsgn(Fs); 
  
  % ADSGN  Design of a A-weighting filter.
  %    [B,A] = ADSGN(Fs) designs a digital A-weighting filter for 
  %    sampling frequency Fs. Usage: Y = FILTER(B,A,X). 
  %    Warning: Fs should normally be higher than 20 kHz. For example, 
  %    Fs = 48000 yields a class 1-compliant filter.
  %
  %    Requires the Signal Processing Toolbox. 
  %
  %    See also ASPEC, CDSGN, CSPEC. 
  
  % Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
  %         couvreur@thor.fpms.ac.be
  % Last modification: Aug. 20, 1997, 10:00am.
  
  % References: 
  %    [1] IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996. 
  
  
  % Definition of analog A-weighting filter according to IEC/CD 1672.
  f1 = 20.598997; 
  f2 = 107.65265;
  f3 = 737.86223;
  f4 = 12194.217;
  A1000 = 1.9997;
  pi = 3.14159265358979;
  NUMs = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];
  DENs = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]); 
  DENs = conv(conv(DENs,[1 2*pi*f3]),[1 2*pi*f2]);
  
  % Use the bilinear transformation to get the digital filter. 
  [B,A] = bilinear(NUMs,DENs,1/Fs);
endfunction

function buildFilterGainTable(B, A, gaintablesize, Fs);

  % Build a filter gain table used to filter signal
  % 
  % copyright (c) Davide Gironi, 2013
  %
  % Released under GPLv3.
  % Please refer to LICENSE file for licensing information.


  %build filter gaintable
  WindowSize = gaintablesize*2; %samples for gaintable creation / window length
  filter_gaintable = freqz(B, A, WindowSize, 'whole'); %built for the sized window
  filter_gaintableR = real(filter_gaintable);
  filter_gaintableI = imag(filter_gaintable);

  T = 1/Fs; %sample time
  L = WindowSize*20; %length of signal (20 * window)
  t = (0:L-1)*T; %time vector
  y = 0.7*sin(2*pi*50*t) + sin(2*pi*2000*t) + 2*randn(size(t)); 
  %sum of a 50Hz and a 2000 Hz sinusoid plus noise
  signal = y';

  %estimate signal level (within each windowed segment).
  windowIntervals = [1:WindowSize:(length(signal)-WindowSize)];
  windowTime = t(windowIntervals+round((WindowSize-1)/2));
  %will contain the full signal
  signalT = [];
  filteredsignalT = [];
  filteredsignalfreqrespT = [];
  %will contain the dBA computations
  dBA1 = zeros(length(windowIntervals),1);
  dBA2 = zeros(length(windowIntervals),1);
  dBA3 = zeros(length(windowIntervals),1);

  for i = [1:length(windowIntervals)]
    %take a part of the signal
      signalP = signal(windowIntervals(i)-1+[1:WindowSize]);
      
      %filter signal part by filter function
      filteredsignalP = filter(B, A, signalP); 

      %filter signal part by frequency response 
      NFFTP = length(signalP);
      fftsignalP = fft(signalP);
      fftfilteredsignalP = fftsignalP .* filter_gaintable;
      filteredsignalfreqrespP = ifft(fftfilteredsignalP, NFFTP);
      
      %built the full signal
      signalT = cat(1, signalT, signalP);
      filteredsignalT = cat(1, filteredsignalT, filteredsignalP);
      filteredsignalfreqrespT = cat(1, filteredsignalfreqrespT, filteredsignalfreqrespP);
      
      %calcuatate rms of the signal part
      rmssignalP = sqrt(mean(signalP.^2));
      rmsfilteredsignalP = sqrt(mean(filteredsignalP.^2));
      rmsfilteredsignalfreqrespP = sqrt(mean(filteredsignalfreqrespP.^2));
      
      % Estimate dBA value using Parseval's relation.
      dBA1(i) = 20*log10(rmssignalP);
      dBA2(i) = 20*log10(rmsfilteredsignalP);
      dBA3(i) = 20*log10(rmsfilteredsignalfreqrespP);
  end

  % %plot full signal, and filtered signal
   % figure(2)
   % plot([1:length(signalT)], signalT, '-', [1:length(filteredsignalT)], filteredsignalT, '-', [1:length(filteredsignalfreqrespT)], filteredsignalfreqrespT);
   % title(['compare signal (1 window), to filtered signal']);
   % legend('signal', 'filtered (using filter function)', 'filtered (using response table)');
   % grid on

  %plot out results
  figure(3)
  plot([1:length(dBA1)],dBA1,'-',[1:length(dBA2)],dBA2,'-',[1:length(dBA3)],dBA3);
  legend('dBA signal','dBA filtered (using filter function)', 'dBA filtered (using weight table)',3);
  grid on

  %save gaintable to xls
  CR(:,1) = filter_gaintableR(1:length(filter_gaintableR)/2); %retain frequencies below Nyquist rate
  CI(:,1) = filter_gaintableI(1:length(filter_gaintableI)/2); %retain frequencies below Nyquist rate
  csvwrite('aweightingtableR.csv', CR);
  csvwrite('aweightingtableI.csv', CI);
endfunction



%frequency used for filter creation
Fsf = 22050; 

%gaintable size
gaintablesize = 32;

%sampling frequency of signal
Fs = Fsf; %22050 

%build A-Weight Filter
[B,A] = adsgn(Fsf); 

%plot A-Weight and compare with Class I and Class II
aspec(B,A,Fsf);

%build the filter
buildFilterGainTable(B, A, gaintablesize, Fs);