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

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

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

FS = 48e3;
t = (1:2*FS)/FS;
s = rng('default');
xtone = 0.046*sin(2*pi*t*1000).' + 1e-2*pinknoise(2*FS);
rng(s)

splMeterReading = 90.1; 

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

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

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

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

Загрузите вертолет, записывающий, который содержит пространственную информацию. Только необходимо сохранить первый канал, который соответствует всенаправленному сигналу (похожий на моно запись с нормальным ненаправленным микрофоном). Создайте соответствующий временной вектор.

[x,FS] = audioread("Heli_16ch_ACN_SN3D.wav");
x = x(:,1);
t = (1:size(x,1))/FS;

Создайте splMeter возразите и выберите C-weighting, быстрое время, взвесив и 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 Helicopter Noise')
xlabel('Time (seconds)')
ylabel('SPL (dB)')
ylim([68 100])
grid on

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

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

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

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

acousticLoudness(x,FS,calib)

Громкость составляет 42,5 сона с пиком в 6 (Шкала коры). Преобразуйте пиковое значение громкости в Гц с помощью bark2hz

fprintf("Peak loudness frequency: %d Hz\n",round(bark2hz(6),-1));

Peak loudness frequency: 630 Hz

По умолчанию уровни находятся в сонах. Может быть легче изучить это измерение, если это сравнилось с SPL (дБ) чтение. Единица фонов громкости является соответствующей к ссылочной частоте 1 кГц. Например, любой сигнал с громкостью 60 фонов, как воспринимают, является столь же громким как тон на 1 кГц, который указывает на уровне 60 дБ на метре SPL (1 кГц незатронут взвешиванием частоты). Следовательно, преобразовывая 42,5 сона в фоны с помощью sone2phon сообщает, что вертолет воспринят столь же громкий как тон на 94 дБ 1 кГц.

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

Equivalent 1 kHz SPL: 94 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

Можно также построить изменяющуюся во времени громкость и определенную громкость, чтобы анализировать звук движущегося вертолета. В этом случае шум является стационарным, но можно наблюдать distictive "импульсивную" природу вертолета.

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

Постройте определенную громкость с частотой в Гц (логарифмическая шкала).

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

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

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

sharp = acousticSharpness(spec)

sharp = 1.4507

Розовый шум имеет резкость 2 acums, таким образом, вы теперь знаете, что вертолетный шум склоняется к низким частотам.

Постройте изменяющуюся во времени резкость.

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')

Заметьте, что частота модуляции достигает максимума на уровне 18,8 Гц. Ниже 30 Гц модуляция воспринята как колебание (в противоположность шероховатости).

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

acousticFluctuation(x,FS,calib)

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

[vacil,spec,fMod] = acousticFluctuation(tvspec);
clf; % do not reuse previous subplot
flucHz = bark2hz(0.5:0.5:23.5);
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

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

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

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

с:

В этом примере резкость была меньше 1.75, таким образом, это не влияющий фактор. Шум в основном модулируется частотой ниже, чем 30 Гц, поэтому считайте шероховатость незначительной для этого примера. Поэтому можно установить $\mathit{R}$ и ${\mathit{w}}_{\mathit{s}}$ обнулять.

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

N5 = perc(2);

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

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

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

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

pa = 46.7171

Эффект защиты уха

Измерьте удар защиты органов слуха на измеренном SPL и воспринятом шуме. Характеристики затухания 24 дБ, оцененных, слыша средство защиты:

Частота (Гц) 125 250 500 1000 2000 3000 4000 6000 8000

Затухание (дБ) 14.4 21.6 27.0 31.9 36.1 39.5 41.9 39.8 37.0

Интерполируйте данные о защите органов слуха.

att = [  125   250   500  1000  2000  3000  4000  6000  8000 ; ...
        14.4  21.6  27.0  31.9  36.1  39.5  41.9  39.8  37.0 ].';

att = [ 20,0; att ]; % Assume zero attenuation at 20 Hz

cf = getCenterFrequencies(splMeter("Bandwidth","1/3 octave"));

p = interp1(att(:,1),att(:,2),cf,"pchip");

Отобразите интерполированное затухание.

semilogx(cf,-p,att(:,1),-att(:,2),'rx')
title('Attenuation from Hearing Protection')
ylabel('Relative SPL (dB)')
xlabel('Frequency (Hz)')
xlim(cf([1,end]))
grid on

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

Спроектируйте graphicEQ объект симулировать эффект защиты органов слуха.

geq = graphicEQ("SampleRate",FS,"Bandwidth","1/3 octave","EQOrder",14,"Gains",-p);

Отобразите частотную характеристику 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 Hearing Protection Simulation Filter')
ylabel('Relative SPL (dB)')
xlabel('Frequency (Hz)')
xlim(cf([1,end]))
grid on

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

xhp = geq(x);

Сравните SPL с и без защиты органов слуха. Снова используйте те же настройки метра SPL, но используйте отфильтрованную запись. Установите пределы Y как прежде, чтобы сделать визуальное сравнение легче с более ранним графиком.

reset(spl)
[LCFnew,~,LCpeaknew] = spl(xhp);
plot(t,LCpeak,t,LCF,t,LCpeaknew,t,LCFnew)
legend('LCpeak (original)', 'LCF (original)', ...
       'LCpeak (with ear protection)', ...
       'LCF (with ear protection)', ...
       'Location','southeast')
title('SPL Measurement of Helicopter Noise')
xlabel('Time (seconds)')
ylabel('SPL (dB)')
ylim([68 100])
grid on

Затухание может быть меньше ожидаемого путем рассмотрения спецификации, но как вы видели с измерением резкости, heliocopter шум сконцентрирован в более низких частотах, где любая защита органов слуха менее эффективна.

Посмотрите на воспринятую громкость, чтобы получить лучшее представление о том, как громкий вертолет теперь звучит.

acousticLoudness(xhp,FS,calib)

Громкость пошла от 42,5 до 6,2 сонов. Однако может быть легче посмотреть на громкость в фонах вместо сонов, потому что фоны эквивалентны SPL тона на 1 кГц. Преобразуйте единицы сона к фонам с помощью sone2phon.

fprintf("Loudness without hearing protection:\t%.1f phons\n",sone2phon(42.5));

Loudness without hearing protection:	94.1 phons

fprintf("Loudness with hearing protection:\t%.1f phons\n",sone2phon(6.2));

Loudness with hearing protection:	66.3 phons

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

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

Эта защита органов слуха оценивается как 24 дБ, но вертолетный шум воспринят как на приблизительно 28 дБ более тихий при износе их.

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

[~,spechp,perchp] = acousticLoudness(xhp,FS,calib,"TimeVarying",true);
shp = acousticSharpness(spechp,'TimeVarying',true);
new_sharp = mean(shp(501:end))

new_sharp = 0.7258

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

vacilhp = acousticFluctuation(spechp);
fhp = mean(vacilhp(501:end,1));

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

N5hp = perchp(2); % N5 with hearing protection
pahp = N5hp * (1 + abs(2.18/(N5hp^.4)*(.4*fhp)))

pahp = 7.0230