Fast time integration troubles

% Importar datos de Excel
[data, ~] = xlsread('Pruebawavwav.xlsx', 'Pruebawav');

% Voltaje de referencia del ADC en V
Vref = 3.3;

% Ganancia del sensor
gain = 25;

% Convertir valores del ADC a voltajes
data_volts = (data / 1024) * Vref;

% Ajustar voltajes por la ganancia del sensor
data_volts_gain = data_volts / gain;

% Frecuencia de muestreo en Hz
fs = 44100;

% Coeficientes de la ponderación A según la norma ANSI S1.42
A_weighting_coeffs = [1.2589, -2.0082, 0.8315];

% Crear filtro de ponderación A con los coeficientes de la norma ANSI S1.42
A_weighting_filter = dsp.FIRFilter('Numerator', A_weighting_coeffs);

% Filtrar datos
data_filtered = A_weighting_filter(data_volts_gain);

% Coeficiente de suavizado del filtro de media móvil exponencial
alpha = 0.1;

% Coeficientes del filtro de media móvil exponencial
b = [alpha];

% Aplicar filtro de media móvil exponencial
data_fast_avg = filter(b, 1, data_filtered);

% Definir la presión sonora de referencia en Pascals
pref = 20e-6;

% Calcular el nivel de presión sonora en dB utilizando los datos filtrados en modo Fast
data_dB = 20 * log10(abs(data_fast_avg) / pref);

% Graficar señal filtrada en modo Fast
figure;
plot(data_dB);
title('Señal Filtrada Ponderación A (Modo Fast)');
xlabel('Ventanas de 125 ms');
ylabel('Nivel de Presión Sonora (dB)');

Im currently working in a project where itake data from a max4466 through a wemos d1 mini, an export it to an excel. After i upload the excel in matlab, the aim is that the code must work as a sound pressure level. For this i have to use an A weighting filte and a FAST time integration.

If someone knows about this topic, i would like if i can get some feedback about my code.

Thanks

[VBBV]

Can you share the data file in your code ? what is the error you get (if any) ?

[Federico]

I don't have any errors. My doubt is if what I am doing is theoretically correct. I would like to know if the temporal integration is well thought out.

[William Rose]

@Federico, I think your code for fast integration is not doing what you want. What you want is to smooth the signal with an exponential time constant of 125 msec. I also doubt that your FIR filter is acheiving the A weighting standard, but I am less certain about that. I also have doubts about the calibration to Pascals and then to dB(SPL). See details in my answer below. I agree with @VBBV that it would help if you posted Pruebawavwav.xlsx, or a fragment of it.

[Federico]

Hi @William Rose, in my last response to@VBBV i attached the file, please let me know if you can open it.

On the other hand, i looked at your last answer and will be anwsering it asap.

Answer 1

[William Rose]

@Federico,

The code you use to convert data to data_volts looks reasonable. Your code assumes that the ADC produces an output of 1024 at 3.3 volts and 0 at 0 volts. (A 10 bit ADC has a maximum numeric output=1023, but that differentce is negligible.) Then you divide data_volts by by gain to get data_volts_gain, which is the voltage before amplification. That looks reasonable to me.

Then you filter the data, by creating a FIR filter and applying it to data_volts_gain. I am surprised that an FIR filter with just 3 coefficients can effectively implement an A-weighting filter. If I were doing it, I would check that part, and I would compare the FIR filter results to results obtained with an IIR implementaiton of the A weighting standard. Here is a paper on getting the coefficients for an IIR version of the A weighting filter. Here is another paper. In your comments, you refer to ANSI standard S 1.42. The ANSI standards are very expensive to obtain, so I have not read it. Here is a site that summarizes the A weighting equations. It also refers to a IIR filter implemnation of the A standard in Matlab. Here is the older ANSI 1.4 standard, which does not cost anything, which describes the A weighting.

Then you create an IIR filter to do the Fast averaging. Fast averaging is an exponential decay filter with a 125 millisecond time constant. The Fast avergaging standard is defined in the IEC 61672-1 standard. This standard, like the ANSI standards, is very expensive, so I have not read it. It is not clear to me whether the exponential decay filter should be applied to the signal measured in Pascals (or its voltage equivalent), or to the sound pressure after it has been converted to dB, or to the intensity (Pascals squared). This site includes an image with Pascals squared on the vertical axis, when describing the Fast and Slow filters. If that is correct, the you should covert the signal to Pascals, then square it, then apply the Fast exponential decay filter. Your current code does not do this.

The exponential decay filter which you use to do the fast averaging does not look right to me. A discrete time version of an exponential decay filter, with time constant tau, is

fs=44100;               % sampling rate (Hz)
dt=1/fs;                % sampling interval (s)
tau=0.125;              % time constant for Fast filter, IEC 61672-1 (s)
b=dt/(tau+dt);          % numerator coefficients for IIR filter
a=[1,-tau/(tau+dt)];    % denominator coefficients for IIR filter
y=filter(b,a,x);        % apply filter to signal x

Your values for b and a are very different, and I think they will not give the desired result. I have used tau=0.125 s above, since that is the value for tau specified by IEC 61672-1, according to various online sources.

After doing the A-weight filtering and the fast filtering, you convert the signal to dB. I think you are trying to convert the signal to dB(SPL) rather than dBV or dBm. If I am correct, then the equation you use does not appear to make sense to me. You do

pref = 20e-6; % Definir la presión sonora de referencia en Pascals
data_dB = 20 * log10(abs(data_fast_avg) / pref);

The code above would make sense if the signal in data_fast_avg already has units of Pascals, because 20e-6 Pascals corresponds to 0 dB SPL. But I am not confident that data_fast_avg does have those units. I recommend that you confirm the calibration factors.

If the signal has units of Pa, and then you square it, before doing fast averaging, which would be consistent with the plot here, then the output from fast averaging will have units of Pa^2. Suppose the output from fast averaging is dataFast , and suppose it has units of Pa^2. Then, to convert to dB SPL, you would do

pref=20e-6;  % pressure in Pa that corresponds to 0 dB SPL
data_dBSPL = 10*log10(dataFast/pref^2);  % convert signal from Pa^2 to dB SPL.

Note that I squared pref in the denominator, and I multipled by 10, not 20.

Good luck.

[William Rose]

@Federico,

I have looked at the xlsx file you uploaded. I see that there are about 23k points, i.e. about 0.5 seconds of data. The values are integers between 273 and 709, with a mean of 542. When we apply the conversions in your code (which are to multiply by 3.3/1024, then divide by 25), the mean value becomes 0.070. If the filtering does not affect the mean value in a major way, then this level would correspond to about 71 dB(SPL). This is a believable value. Therefore I am no longer as concerned as I was before about the calibration of the signal. I still think the A weighting filter needs attention - see the two articles I cited. I recommend you use the approach I described for the Fast filter.

I always check my filtering code by creating simulated signals with known frequency content. Then I see if my code does what I expect to the signals.

[William Rose]

@Federico, Here is the frequency repsonse of the A_weighting_filter in your script.

I plotted the frequency response with the command

fvtool(A_weighting_filter,'Fs',fs)

The frequency response above is not the A_weighting filter response. Therefore I recommend implementing the A weighting filter by following one of the two references which I provided in my answer.

I said earlier that one source shows that the pressure signal is squared before applying the Fast filter. I now realize another reason why one should square the signal before applying the Fast filter. The reason is that, after applying the A filter, the signal will have rapid fluctuations and zero mean. If you smooth a zero mean signal, you get very little. So you square it first, to make the negative parts positive. Then you smooth it with the Fast filter. This is analogous to full or half wave rectification of a signal before lowpass filtering - like a simple AM radio receiver. It is also analogous to extracting the envelope of the electromyogram (EMG) signal.

[William Rose]

@Federico, the A-weighting filter is built in to the Matlab Audio toobox. Therefore, if you have the Audio Toolbox, do

AWeighting = weightingFilter('A-weighting',fs);
audioOut = Aweighting(audioIn)

I don't have the Audio toolbox, so the weightingFilter() command does not work for me.

[William Rose]

If you don't have the Audio Toolbox, search for "A-weighting" on the Matlab File Exchange. Multiple users have uploaded scripts and libraries that implement the A-weighting filter.

[William Rose]

@Federico, try attached script, and see comments in the script. I think it does what you want.

[Federico]

Hi @William Rose im still trying to understand the code and also testing it. I don t know if it is your area of expertise, but i would like to know if it is a way to convert this code into an arduino compatible one. In a future i would like to be able to measure in real time. The tricky part is the .ino code needs to be able to run in a Wemos D1 Mini.

PS: Thanks for your time and help!

[Federico]

Regarding the A-weighting filter, i do have the Audio Toolbox. My doubt is if im using the Excel data, would it be compatible to the weightingFilter() command?

[Federico]

Another observation is the plotted values of data_dB in the third subplot:

The values shown are quite low compared to a SPL phone app. I still dont get why the data_dB equation its like this:

data_dB = 10 * log10(data_fast_avg / pref^2);

instead of:

data_dB = 20 * log10(data_fast_avg / pref^2);

With the second one i get the following:

Here, the dB make more sense comared with the measurements taken with the phone app.

I know this has a lot to do with the fact that I now operate in Pa but i dont understad all the way.

[William Rose]

@Federico,

You wrote:

I would like to know if it is a way to convert this code into an arduino compatible one. In a future i would like to be able to measure in real time. The tricky part is the .ino code needs to be able to run in a Wemos D1 Mini.

I am not able to assist with converting the code to run on an Arduino.

You wrote:

I still dont get why the data_dB equation its like this:

data_dB = 10 * log10(data_fast_avg / pref^2);

instead of:

data_dB = 20 * log10(data_fast_avg / pref^2);

Note that data_fast_avg, in my code, is filtered version of the squared pressure.

The equation for SPL in dB is

where p is the root mean square pressure in Pa, and .

Due to the mathematics of logarithms, this is equivalent to

and data_fast_avg is the smothed, or averaged, versino of the squared pressure, which means data_fast_avg is analogous to mean square pressure.

[William Rose]

@Federico,

You wrote:

I do have the Audio Toolbox. My doubt is if im using the Excel data, would it be compatible to the weightingFilter() command?

I am glad to hear you have the Audio Toolbox. You can use data imported from Excel with the weightingFilter() funciton. Once you have imported the data, it doesn't matter whether it came from a .txt file or an Excel file or a .csv file or some other file type.

I'm sorry if I misunderstood your question.