acousticFluctuation

Perceived fluctuation strength of acoustic signal

Syntax

fluctuation = acousticFluctuation(audioIn,fs)

fluctuation = acousticFluctuation(audioIn,fs,calibrationFactor)

fluctuation = acousticFluctuation(specificLoudnessIn)

fluctuation = acousticFluctuation(___,Name,Value)

[fluctuation,specificFluctuation] = acousticFluctuation(___)

[fluctuation,specificFluctuation,fMod] = acousticFluctuation(___)

acousticFluctuation(___)

Description

fluctuation = acousticFluctuation(audioIn,fs) returns fluctuation strength in vacil based on Zwicker et al. [1] and ISO 532-1 time-varying loudness [2].

fluctuation = acousticFluctuation(audioIn,fs,calibrationFactor) specifies a nondefault microphone calibration factor used to compute loudness.

fluctuation = acousticFluctuation(specificLoudnessIn) computes fluctuation using time-varying specific loudness.

fluctuation = acousticFluctuation(___,Name,Value) specifies options using one or more Name,Value pair arguments.

Example: fluctuation = acousticFluctuation(audioIn,fs,'SoundField','diffuse') returns fluctuation assuming a diffuse sound field.

[fluctuation,specificFluctuation] = acousticFluctuation(___) also returns specific fluctuation strength.

[fluctuation,specificFluctuation,fMod] = acousticFluctuation(___) also returns the dominant modulation frequency.

acousticFluctuation(___) with no output arguments plots fluctuation strength and specific fluctuation strength and displays the modulation frequency textually. If the input is stereo, the 3-D plot shows the sum of both channels.

Examples

Measure Acoustic Fluctuation

Measure acoustic fluctuation based on Zwicker et al [2] and ISO 532-1 [1]. Assume the recording level is calibrated such that a 1 kHz tone registers as 100 dB on an SPL meter.

[audioIn,fs] = audioread('WashingMachine-16-44p1-stereo-10secs.wav');

fluctuation = acousticFluctuation(audioIn,fs);

Fluctuation Measurements Using Calibrated Microphone

Set up an experiment as indicated by the diagram.

Create an audioDeviceReader object to read from the microphone and an audioDeviceWriter object to write to your speaker.

fs = 48e3;
deviceReader = audioDeviceReader(fs,"SamplesPerFrame",2048);
deviceWriter = audioDeviceWriter(fs);

Create an audioOscillator object to generate a 1 kHz sinusoid.

osc = audioOscillator("sine",1e3,"SampleRate",fs,"SamplesPerFrame",2048);

Create a dsp.AsyncBuffer object to buffer data acquired from the microphone.

dur = 5;
buff = dsp.AsyncBuffer(dur*fs);

For five seconds, play the sinusoid through your speaker and record using your microphone. While the audio streams, note the loudness as reported by your SPL meter. Once complete, read the contents of the buffer object.

numFrames = dur*(fs/osc.SamplesPerFrame);
for ii = 1:numFrames
    audioOut = osc();
    deviceWriter(audioOut);
    
    audioIn = deviceReader();
    write(buff,audioIn);
end

SPLreading = 60.4;

micRecording = read(buff);

To compute the calibration factor for the microphone, use the calibrateMicrophone function.

calibrationFactor = calibrateMicrophone(micRecording(fs+1:end,:),deviceReader.SampleRate,SPLreading);

You can now use the calibration factor you determined to measure the fluctuation of any sound that is acquired through the same microphone recording chain.

Perform the experiment again, this time, add 100% amplitude modulation at 4 Hz. To create the modulation signal, use audioOscillator and specify the amplitude as 0.5 and the DC offset as 0.5 to oscillate between 0 and 1.

mod = audioOscillator("sine",4,"SampleRate",fs, ...
    "Amplitude",0.5,"DCOffset",0.5,"SamplesPerFrame",2048);

dur = 5;
buff = dsp.AsyncBuffer(dur*fs);
numFrames = dur*(fs/osc.SamplesPerFrame);
for ii = 1:numFrames
    audioOut = osc().*mod();
    deviceWriter(audioOut);
    
    audioIn = deviceReader();
    write(buff,audioIn);
end

micRecording = read(buff);

Call acousticFluctuation with the microphone recording, sample rate, and calibration factor. The fluctuation reported from acousticFluctuation uses the true acoustic loudness measurement as specified by 532-1. Display the average fluctuation strength over the 5 seconds.

fluctuation = acousticFluctuation(micRecording,deviceReader.SampleRate,calibrationFactor);
fprintf('Average fluctuation = %d (vacil)',mean(fluctuation(501:end,:)))

Average fluctuation = 1.413824e+00 (vacil)

acousticFluctuation(micRecording,deviceReader.SampleRate,calibrationFactor)

Measure Fluctuation from Specific Loudness

Read in an audio file.

[audioIn,fs] = audioread("Engine-16-44p1-stereo-20sec.wav");

Call acousticLoudness to calculate the specific loudness.

[~,specificLoudness] = acousticLoudness(audioIn,fs,'TimeVarying',true);

Call acousticSharpness without any outputs to plot the acoustic sharpness.

acousticSharpness(specificLoudness,'TimeVarying',true)

Call acousticFluctuation without any outputs to plot the acoustic fluctuation.

acousticFluctuation(specificLoudness)

Effect of Frequency Modulation on Acoustic Fluctuation

Generate a pure tone with a 1500 Hz center frequency and approximately 700 Hz frequency deviation at a modulation frequency of 0.25 Hz.

fs = 48e3;

fMod = 0.25;
dur = 20;

numSamples = dur*fs;
t = (0:numSamples-1)/fs;

tone = sin(2*pi*t*fMod)';

fc = 1500;
excursionRatio = 0.47;

excursion = 2*pi*(fc*excursionRatio/fs);
audioIn = modulate(tone,fc,fs,'fm',excursion);

Listen to the first 5 seconds of the audio and plot the spectrogram.

sound(audioIn(1:5*fs),fs)
spectrogram(audioIn(1:5*fs),hann(512,'periodic'),256,1024,fs,'yaxis')

Call acousticFluctuation with no output arguments to plot the acoustic fluctuation strength.

acousticFluctuation(audioIn,fs);

Specify Known Modulation Frequency

The acousticFluctuation function enables you to specify a known fluctuation frequency. If you do not specify a known fluctuation frequency, the function auto-detects the fluctuation.

Create a dsp.AudioFileReader object to read in an audio signal frame-by-frame. Create an audioOscillator object to create a modulation wave. Apply the modulation wave to the audio file.

fileReader = dsp.AudioFileReader('Engine-16-44p1-stereo-20sec.wav');

fmod = 10.8;
amplitude = 0.15;

osc = audioOscillator('sine',fmod, ...
    "DCOffset",0.5, ...
    "Amplitude",amplitude, ...
    "SampleRate",fileReader.SampleRate, ...
    "SamplesPerFrame",fileReader.SamplesPerFrame);

testSignal = [];
while ~isDone(fileReader)
    x = fileReader();
    testSignal = [testSignal;osc().*fileReader()];
end

Listen to two seconds of the test signal and plot its waveform.

samplesToView = 1:2*fileReader.SampleRate;
sound(testSignal(samplesToView,:),fileReader.SampleRate);

plot(samplesToView/fileReader.SampleRate,testSignal(samplesToView,:))
xlabel('Time (s)')

Plot the acoustic fluctuation. The detected frequency of the modulation is displayed textually.

acousticFluctuation(testSignal,fileReader.SampleRate);

Specify the known modulation frequency and then plot the acoustic fluctuation again.

acousticFluctuation(testSignal,fileReader.SampleRate,'ModulationFrequency',fmod)

Input Arguments

audioIn — Audio input

column vector | two-column matrix

Audio input, specified as a column vector (mono) or matrix with two columns (stereo).

Tip

To measure fluctuation strength given any modulation frequency, the recommended minimum signal duration is 10 seconds.

Data Types: single | double

fs — Sample rate (Hz)

positive scalar

Sample rate in Hz, specified as a positive scalar. The recommended sample rate for new recordings is 48 kHz.

Note

The minimum acceptable sample rate is 8 kHz.

Data Types: single | double

calibrationFactor — Microphone calibration factor

sqrt(8) (default) | positive scalar

Microphone calibration factor, specified as a positive scalar. The default calibration factor corresponds to a full-scale 1 kHz sine wave with a sound pressure level of 100 dB (SPL). To compute the calibration factor specific to your system, use the calibrateMicrophone function.

Data Types: single | double

specificLoudnessIn — Specific loudness (sones/Bark)

T-by-240-by-C

Specific loudness in sones/Bark, specified as a T-by-240-by-C array, where:

• T is one per 2 ms of the time-varying signal.

• 240 is the number of Bark bins in the domain for specific loudness. The Bark bins are 0.1:0.1:24.

• C is the number of channels.

You can use the acousticLoudness function to calculate time-varying specific loudness using this syntax:

[~,specificLoudness] = acousticLoudness(audioIn,fs,'TimeVarying',true);

Data Types: single | double

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: acousticFluctuation(audioIn,fs,'ModulationFrequency',50)

ModulationFrequency — Known modulation frequency (Hz)

'auto-detect' (default) | scalar or two-element vector with values in the range [0.1,100]

Known modulation frequency in Hz, specified either 'auto-detect' or as a scalar or two-element vector with values in the range [0.1,100]. If ModulationFrequency is set to 'auto-detect', then the function limits the search range to between 0.2 and 64 Hz. If the input is mono, then the modulation frequency must be specified as a scalar. If the input is stereo, then the modulation frequency can be specified as either a scalar or two-element vector.

Data Types: single | double | char | string

SoundField — Sound field

'free' (default) | 'diffuse'

Sound field of audio recording, specified as 'free' or 'diffuse'.

Data Types: char | string

PressureReference — Reference pressure (Pa)

20e-6 (default) | positive scalar

Reference pressure for dB calculation in pascals, specified as a positive scalar. The default value, 20 micropascals, is the common value of air.

Data Types: single | double

Output Arguments

fluctuation — Fluctuation strength (vacil)

K-by-1 | K-by-2

Fluctuation strength in vacil, returned as a K-by-1 column vector or K-by-2 matrix of independent channels. K corresponds to the time dimension.

Data Types: single | double

specificFluctuation — Specific fluctuation strength (vacil/Bark)

K-by-47 matrix | K-by-47-by-2 array

Specific fluctuation strength in vacil/Bark, returned as a K-by-47 matrix or a K-by-47-by-2 array. The first dimension of specificFluctation, K, corresponds to the time dimension and matches the first dimension of fluctuation. The second dimension of specificFluctation, 47, corresponds to bands on the Bark scale, with centers from 0.5 to 23.5, inclusive, in 0.5 increments. The third dimension of specificFluctation corresponds to the number of channels and matches the second dimension of fluctuation.

Data Types: single | double

fMod — Dominant modulation frequency (Hz)

scalar (mono input) | 1-by-2 vector (stereo input)

Dominant modulation frequency in Hz, returned as a scalar for mono input or a 1-by-2 vector for stereo input.

Data Types: single | double

Algorithms

Acoustic fluctuation strength is a perceptual measurement of slow modulations in amplitude or frequency. The acoustic loudness algorithm is described in [1] and implemented in the acousticLoudness function. The acoustic fluctuation calculation is described in [2]. The algorithm for acoustic fluctuation is outlined as follows.

$$fluctuation=\frac{0.008{\displaystyle \underset{z=0}{\overset{24}{\int }}\Delta Ldz}}{\left(\frac{f_{\mathrm{mod}}}{4}\right)+\left(\frac{4}{f_{\mathrm{mod}}}\right)}$$

Where f_mod is the detected or known modulation frequency and ΔL is the perceived modulation depth. If the modulation frequency is not specified when calling acousticFluctuation, it is auto-detected by peak-picking a frequency-domain representation of the acoustic loudness. The perceived modulation depth, ΔL, is calculated by passing rectified specific loudness bands through ½ octave filters centered around f_mod, followed by a lowpass filter to determine the envelope.

References

[1] ISO 532-1:2017(E). "Acoustics – Methods for calculating loudness – Part 1: Zwicker method." International Organization for Standardization.

[2] Zwicker, Eberhard, and H. Fastl. Psychoacoustics: Facts and Models. 2nd updated ed, Springer, 1999.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Version History

Introduced in R2020b