Калибровка уровня записи

Психоакустические измерения дают наиболее точные результаты с калиброванным входным уровнем микрофона. Как использовать calibrateMicrophone чтобы соответствовать уровню записи показанию SPL-счетчика, можно использовать источник тонального сигнала 1 кГц (например, генератор тонального сигнала или приложение сотового телефона) и калиброванный SPL-счетчик. SPL калибровочного тонального сигнала 1 кГц должен быть достаточно громким, чтобы доминировать в любом фоновом шуме. Пример калибровки, в котором в качестве источника тонального сигнала 1 кГц используется MATLAB, см. calibrateMicrophone.

Симулируйте запись тонального сигнала и включите некоторый фоновый шум. Предположим, что значение счетчика SPL составляет 83,1 дБ (взвешенный по C).

FS = 48e3;
t = (1:2*FS)/FS;
s = rng('default');
testTone = 0.46*sin(2*pi*t*1000).' + .1*pinknoise(2*FS);
rng(s)

splMeterReading = 83.1;

Чтобы вычислить уровень калибровки цепи записи, вызовите calibrateMicrophone и укажите тестовый тональный сигнал, частоту дискретизации, показания SPL и частотное взвешивание SPL-счетчика. Чтобы компенсировать возможный фоновый шум и получить точный уровень калибровки, совпадайте с настройкой взвешивания частоты SPL-счетчика.

calib = calibrateMicrophone(testTone,FS,splMeterReading,"FrequencyWeighting","C-weighting");

Уровни звукового давления (SPL)

Если у вас есть коэффициент калибровки для вашей цепи записи, вы можете сделать акустические измерения. При использовании физического счетчика вы ограничены настройками, выбранными во время измерения. С splMeter объект можно изменить после записи. Это облегчает эксперименты с различными опциями взвешивания по времени и частоте.

Загрузите запись двигателя и создайте соответствующий временной вектор.

[x,FS] = audioread('Engine-16-44p1-stereo-20sec.wav');
x = x(1:8*FS,1); % use only channel 1 and keep only 8 seconds.
t = (1:size(x,1))/FS;

Создайте splMeter объект и выберите C-взвешивание, быстрое взвешивание во времени и интервал 0,2 секунды для пикового измерения SPL.

spl = splMeter("CalibrationFactor",calib, ...
               "FrequencyWeighting","C-weighting", ...
               "TimeWeighting","Fast", ...
               "TimeInterval",0.2, ...
               "SampleRate",FS);

Постройте графики SPL и пиковых SPL.

[LCF,~,LCpeak] = spl(x);
plot(t,LCpeak,t,LCF)
legend('LCpeak','LCF','Location',"southeast")
title('SPL Measurement of Engine Noise')
xlabel('Time (seconds)')
ylabel('SPL (dB)')
ylim([70 95])
grid on

Психоакустические метрики

Уровень громкости

Мониторинг SPL важен для податливости безопасности труда. Однако измерения SPL не отражают громкость, воспринимаемую фактическим прослушивателем. acousticLoudness измеряет уровни громкости, воспринимаемые человеческим прослушивателем с нормальным слухом (без нарушений слуха). The acousticLoudness функция также показывает, какие полосы частот больше всего способствуют восприятию громкости.

Используя тот же уровень калибровки, что и прежде, и принимая запись в свободное поле (по умолчанию), постройте стационарную громкость.

acousticLoudness(x,FS,calib)

Громкость составляет 23,8 сона, а большая часть шума - ниже 3,3 (шкала Барка). Перевести Барк 3,3 в Гц используя bark2hz

fprintf("Loudness frequency of 3.3 Bark: %d Hz\n",round(bark2hz(3.3),-1));

Loudness frequency of 3.3 Bark: 330 Hz

The acousticLoudness функция возвращает воспринимаемую громкость в сонах. Чтобы понять измерение сигнала, сравните его с показанием SPL (дБ). Сигнал с громкостью 60 фонов воспринимается таким же громким, как тон 1 кГц, измеренный в 60 дБ SPL. Перевод 23.8 sones в телефоны с использованием sone2phon демонстрирует громкость восприятия шума двигателя так же громко, как тон 1 кГц, измеренный в 86 дБ SPL.

fprintf("Equivalent 1 kHz SPL: %d phons\n", round(sone2phon(23.8)));

Equivalent 1 kHz SPL: 86 phons

Сделайте свой собственный график с модулями в фонах и частотой в Гц на журнал шкале.

[sone,spec] = acousticLoudness(x,FS,calib);
barks = 0.1:0.1:24; % Bark scale for ISO 532-1 loudness
hz = bark2hz(barks);
specPhon = sone2phon(spec);
semilogx(hz,specPhon)
title('Specific Loudness')
subtitle(sprintf('Loudness = %.1f phons',sone2phon(sone)))
xlabel('Frequency (Hz)')
ylabel('Loudness (phons/Bark)')
xlim(hz([1,end]))
grid on

Можно также построить график изменяющейся во времени громкости и определенной громкости, чтобы анализировать звук двигателя, если он изменяется со временем. Это может отображаться с другими соответствующими изменяющимися во времени данными, такими как обороты двигателя в минуту (RPM). В этом случае шум является стационарным, но можно наблюдать импульсивный характер шума.

acousticLoudness(x,FS,calib,'TimeVarying',true,'TimeResolution','high')

Постройте график конкретной громкости с частотой в Гц (журнал).

[tvsoneHD,tvspecHD,perc] = acousticLoudness(x,FS,calib,'TimeVarying',true,'TimeResolution','high');
tvspec = tvspecHD(1:4:end,:,:); % for standard resolution measurements
spectimeHD = 0:5e-4:5e-4*(size(tvspecHD,1)-1); % time axis for loudness output
clf; % do not reuse the previous subplot
surf(spectimeHD,hz,sone2phon(tvspecHD).','EdgeColor','interp');
set(gca,'View',[0 90],'YScale','log','YLim',hz([1,end]));
title('Specific Loudness (HD)')
zlabel('Specific Loudness (phons/Bark)')
ylabel('Frequency (Hz)')
xlabel('Time (seconds)')
colorbar

Уровень резкости

Воспринимаемая резкость звука может существенно способствовать его общему уровню раздражения. Оцените воспринимаемый уровень резкости акустического сигнала, используя acousticSharpness функция.

sharp = acousticSharpness(spec)

sharp = 1.1512

Розовый шум имеет резкость 2 акума. Это означает, что шум двигателя смещен к низким частотам.

Постройте график изменяющейся во времени резкости.

acousticSharpness(x,FS,calib,'TimeVarying',true);

Сила колебаний

В случае шума двигателя низкочастотные модуляции способствуют воспринимаемому уровню раздражения.

Во-первых, посмотрите, какие частоты модуляции присутствуют в сигнале.

N = 2^nextpow2(size(x,1));
xa = abs(x); % Use the rectified signal
pspectrum(xa-mean(xa),FS,'FrequencyLimits',[0 80],'FrequencyResolution',1)
title('Modulation Frequencies')

Частота модуляции достигает пика 24,9 Гц. Ниже 30 Гц модуляция воспринимается доминирующим образом как колебание. Есть второй пик на 49,7 Гц, который находится в области значений шероховатостей.

Использование acousticFluctuation для вычисления воспринимаемой силы колебаний. Шум двигателя относительно постоянен в этой записи, поэтому у нас есть алгоритм, автоматически обнаруживающий наиболее слышимую частоту колебаний (fMod).

acousticFluctuation(x,FS,calib)

Интерпретируйте результаты в Hertz вместо Bark. Чтобы уменьшить расчеты, повторно используйте ранее вычисленную изменяющуюся во времени специфическую громкость. Кроме того, можно также задать частоту модуляции, которая вас интересует.

[vacil,spec,fMod] = acousticFluctuation(tvspec,'ModulationFrequency',24.9);
clf; % do not reuse previous subplot
flucHz = bark2hz(0.5:0.5:23.5);
spectime = 0:2e-3:2e-3*(size(spec,1)-1);
surf(spectime,flucHz,spec.','EdgeColor','interp');
set(gca,'View',[0 90],'YScale','log','YLim',flucHz([1,end]));
title('Specific Flucutation Strength')
zlabel('Specific Flucutation Strength (vacils/Bark)')
ylabel('Frequency (Hz)')
xlabel('Time (seconds)')
colorbar

Грубость

Используйте acousticRoughness функция для вычисления воспринимаемой шероховатости сигнала. Пусть алгоритм автоматически обнаруживает самую слышимую частоту модуляции (fMod).

acousticRoughness(x,FS,calib)

Интерпретируйте результаты в Hertz вместо Bark. Чтобы уменьшить расчеты, повторно используйте ранее вычисленную изменяющуюся во времени специфическую громкость. Задайте частоту модуляции.

[asper,specR,fModR] = acousticRoughness(tvspecHD,'ModulationFrequency',49.7);
clf; % do not reuse previous subplot
rougHz = bark2hz(0.5:0.5:23.5);
surf(spectimeHD,rougHz,specR.','EdgeColor','interp');
set(gca,'View',[0 90],'YScale','log','YLim',rougHz([1,end]));
title('Specific Roughness')
zlabel('Specific Roughness (aspers/Bark)')
ylabel('Frequency (Hz)')
xlabel('Time (seconds)')
colorbar

Качество звука

Для общей оценки качества звука объедините предыдущие метрики, чтобы получить метрику психоакустического раздражения (заданную Zwicker и Fastl). Отношение следующее:

Количественное описание было разработано с использованием результатов психоакустических экспериментов:

с:

В этом примере резкость была меньше 1,75, поэтому она не является способствующим фактором. Поэтому можно задать ${\mathit{w}}_{\mathit{s}}$ в нуль.

Процентильная громкость, ${\mathit{N}}_5$, - второе значение, возвращаемое третьим выходом acousticLoudness когда "TimeVarying" установлено в true.

N5 = perc(2);

Вычислите среднюю силу флуктуации, игнорируя первую секунду сигнала.

f = mean(vacil(501:end,1));

Вычислите среднюю шероховатость, игнорирующую первую секунду сигнала.

r = mean(asper(2001:end,1));

Вычислите психоакустическую метрику раздражения.

pa = N5 * (1 + abs(2.18/(N5^.4)*(.4*f+.6*r)))

pa = 25.8848

Эффект улучшенной звукоизоляции

Измерьте влияние улучшенной звукоизоляции на измеренный SPL и воспринимаемый шум.

Симуляция с использованием графического банка EQ-фильтров

Проектирование graphicEQ объект для моделирования ослабления предлагаемой звукоизоляции. Низкие частоты сложнее ослабить, поэтому мы создаем модель, которая лучше всего превышает 200 Гц.

geq = graphicEQ("Bandwidth","1 octave","SampleRate",FS,"Gains",[-0.5 -1.25 -3.4 -7 -8.25 -8.4 -8 -7 -6.4 -5.6]);
cf = getCenterFrequencies(splMeter("Bandwidth","1 octave"));

Отобразите частотную характеристику graphicEQ объект.

[B,A] = coeffs(geq);
sos = [B;[ones(1,size(A,2));A]].';
[H,w] = freqz(sos,2^16,FS);
semilogx(w,db(abs(H)))
title('Frequency Response of Soundproofing Simulation Filter')
ylabel('Relative SPL (dB)')
xlabel('Frequency (Hz)')
xlim(cf([1,end]))
grid on

Отфильтруйте запись двигателя с помощью графического EQ, чтобы симулировать звукоизоляцию.

x2 = geq(x);

Сравните SPL со звукоизоляцией и без нее. Повторно используйте те же настройки SPL-счетчика, но используйте отфильтрованную запись.

reset(spl)
[LCFnew,~,LCpeaknew] = spl(x2);
plot(t,LCpeak,t,LCF,t,LCpeaknew,t,LCFnew)
legend('LCpeak (original)', 'LCF (original)', ...
       'LCpeak (with soundproofing)', ...
       'LCF (with soundproofing)', ...
       'Location','southeast')
title('SPL Measurement of Engine Noise')
xlabel('Time (seconds)')
ylabel('SPL (dB)')
ylim([70 95])
grid on

Сравните воспринимаемые измерения громкости до и после звукоизоляции.

acousticLoudness(x2,FS,calib)

Громкость снизилась с 23,8 до 16,3 соны. Однако может быть легче интерпретировать громкость в фонах. Преобразуйте модули sone в телефоны с помощью sone2phon.

fprintf("Loudness without soundproofing:   \t%.1f phons\n",sone2phon(23.8));

Loudness without soundproofing:   	85.7 phons

fprintf("Loudness with added soundproofing:\t%.1f phons\n",sone2phon(16.3));

Loudness with added soundproofing:	80.3 phons

fprintf("Perceived noise reduction:\t\t%.1f phons (dB SPL at 1 kHz)\n",sone2phon(23.8)-sone2phon(16.3));

Perceived noise reduction:		5.5 phons (dB SPL at 1 kHz)

После звукоизоляции, acousticLoudness показывает, что восприятие шума двигателя примерно на 5,5 дБ тише. Человеческое восприятие звука ограничено на очень низких частотах, где большая часть шума двигателя. Звукоизоляция более эффективна на более высоких частотах.

Вычислите уменьшение фактора психоакустического раздражения. Начните с вычисления среднего значения акустической резкости.

[~,spec2hd,perc2] = acousticLoudness(x2,FS,calib,"TimeVarying",true,"TimeResolution","high");
spec2 = spec2hd(1:4:end,:,:);
shp = acousticSharpness(spec2,'TimeVarying',true);
new_sharp = mean(shp(501:end))

new_sharp = 1.0796

Резкость уменьшилась, потому что звукоизоляция более эффективна при высокой частоте ослаблении. Он ниже порога 1,75, поэтому его игнорируют за коэффициент раздражения.

Теперь вычислите среднее значение прочности и шероховатости колебаний.

vacil2 = acousticFluctuation(spec2);
f2 = mean(vacil2(501:end,1));
asper2 = acousticRoughness(spec2hd);
r2 = mean(asper2(2001:end,1));

Вычислите новый коэффициент психоакустического раздражения. Он сократился, с 25,9 до 17,7.

N5hp = perc2(2); % N5 with soundproofing
pahp = N5hp * (1 + abs(2.18/(N5hp^.4)*(.4*f2+.6*r2)))

pahp = 17.7364

Ссылки

[1] Zwicker, Eberhard, and Hugo Fastl. Психоакустика: факты и модели. Том 22. Springer Science & Business Media, 2013.