A-weighting Sound Filter

I am trying to obtain an A-weighted sound power level (SPL) in dB by using the fdesign.audioweighting function available in the DSP System toolbox. My beginning data are pressure values from multiple microphones which recorded the data at 25.6 kHz. I know all of the equations to calculate unweighted SPL and my calculated values are correct for that case. However, I have tried applying the A-weighting filter at multiple points in the conversion of data from pressure to SPL with no success. I have also tried calculating an FFT and applying the filter to the data in the frequency domain. All of the help descriptions and tutorials end at using fvtool to visualize the filter, but do not describe how to implement the filter on real data. Any insights would be extremely helpful.

Thanks in advance.

[Daniel Shub]

Can you post some example code showing what you are trying, what you expect, and what you get.

[Kevin]

Below is a simplified version of my code. Rather than analyzing the raw pressure data, I model it as a sine wave (close to reality in my application) with the same time between samples and roughly the same magnitude.

As you can see, I tried applying the filter to the data in multiple places using the sosMatrix produced by the audioweighting function. If I'm looking at the weighting correctly, at a frequency of 50 Hz, A-weighting should decrease my sound pressure by ~35dB. Since the unweighted sound pressure should be ~68dB for a pressure of 0.05Pa, the A-weighted sound pressure should be ~33dB. I do not get this result for any of the cases I have tried. Also, I do not see energy conservation between the time and frequency domains (MY vs. My in the code) due to subtracting out the mean value in the FFT (I see energy conservation if I eliminate the offset in the sine wave.) If the mean is not subtracted, the energy is not zero at zero frequency, which it should be.

The main questions I have is how to use the filter constructed using fdesign.audioweighting? I know the contents of the filter, but cannot find any relevant information on how to use the contents during implementation. Also, how do I add the mean pressure into the FFT values so this energy is not lost?

Thank you for the help. Please let me know if you need any more information.

-Kevin

clear; clc;

%%Constants
Fs  = 25600;                    % Sampling frequency (Hz)
dt  = 1/Fs;                     % Sample time (sec)
Ns  = 512;                      % Length of signal (one wavelength)
T   = dt*Ns;                    % Total time
t   = (1:Ns)*dt;                % Time vector
d   = 1.5;                      % Distance from source to microphone (m)
S0 = 1;                        % Reference surface area (m^2)
S1 = 4 * pi() * d^2;           % Measurement surface area (m^2)
Pref = 2e-5                    % Reference pressure

%%Unweighted Sound Power
y = .05 + 0.01*sin(2*pi*50*t);              % Idealized pressure wave (Pa)
y_dB = 10 .* log10(y.^2 ./ Pref^2);         % wave in units of dB (sound pressure)

Lp = 10 * log10((sum(y .^2)/Ns)/Pref^2);    % Sound pressure level (dB)
Lw = Lp + 10 * log10(S1/S0);                % Sound power level (dB)

%%FFT
NFFT = 2^nextpow2(Ns*10);               % Zero-padding for FFT
f = Fs/2*linspace(1/NFFT,1,NFFT/2+1);   % FFT frequency vector
Y = abs((fft(y-mean(y),NFFT))/Ns);      % Frequency distribution of signal

%%Audio Weighting
HawfA = fdesign.audioweighting('WT,Class','A',1,25.6e3);
Af = design(HawfA);
y_filt    = sosfilt(Af.sosMatrix,y);          % Filter pressure data (time domain)
ydB_filt  = sosfilt(Af.sosMatrix,y_dB);       % Filter sound pressure data
Y_filt    = sosfilt(Af.sosMatrix,Y);          % Filter data (freq domain)

% Sound pressure of filtered data
Lp_y   = 10 * log10((sum(y_filt   .^2)/Ns)/Pref^2);  
Lp_ydB = 10 * log10((sum(ydB_filt .^2)/Ns)/Pref^2);
Lp_Y   = 10 * log10((sum(Y_filt   .^2)/Ns)/Pref^2);

% Sound power of filtered data
Lw_y   = Lp_y   + 10 * log10(S1/S0);              
Lw_ydB = Lp_ydB + 10 * log10(S1/S0);
Lw_Y   = Lp_Y   + 10 * log10(S1/S0);

% Check energy conservation between tim and frequency domain
My = sum(y.*conj(y))*dt;
MY = sum(Y.*conj(Y))*Fs/NFFT;

%%Plots
figure; plot(t,y,'r',t,y_filt,'b');                 % Plot ideal wave (Pa)
title('Raw Signal'); xlabel('time (seconds)')
figure; plot(f,Y(1:NFFT/2+1),f,Y_filt(1:NFFT/2+1)); % Plot FFT
title('FFT'); xlabel('frequency (Hz)')

Answer 1

[Daniel Shub]

Your example code is for a single period of a sine wave. The filter you design, however, rings for a period of time. You want to look at the long term average. Changing your definition of Ns to

Ns  = 512*1e3;

I think gives you what you would expect ...

[Kevin]

Thanks for the prompt reply! I tried that and it seems to change the data that was A-weighted in the frequency domain (Lp_Y) to the expected value of ~33dB. However, this seems to be a coincidence (or something special about 1000 wavelengths) since the weighted sound pressures do not converge to this value. Is there a specific reason 1000 wavelengths of data (rather than 100 or 500 or 2000) works properly? Is weighting in the frequency domain (rather than time) the proper way to use this filter?

Also, while I agree that looking at the wave for a longer duration would improve accuracy due to the numerical methods used, I wouldn't have guessed that the calculated values be orders of magnitude different as they are in this case. Do I have to use a large number of wavelengths to use the filter? Because it is not possible for me to get that length for my actual data.

Thanks again for your help. -Kevin

[Davide]

i do not know i can do this, if not please forget this post :) anyway, i've a issue on a-weightening, you can find it here http://www.mathworks.it/matlabcentral/answers/85788-a-weighting-filter-to-lookup-table can you look at this? thank you, davide