2018-01-08
235% Часть 1. РАСЧЕТ РАЗБОРЧИВОСТИ РЕЧИ ПО ИЗВЕСТНЫМ СПЕКТРАЛЬНЫМ УРОВНЯМ РЕЧЕВОГО СИГНАЛА В ОКТАВНЫХ ПОЛОСАХ
i=7; %i-количество октавных полос
fc=[125 250 500 1000 2000 4000 8000] % Среднегеометрическая частота полосы,Гц
fN=90; % нижняя частота диапазона речи, Гц
fV=11200; % верхняя частота диапазона речи, Гц
% используем в расчете справочные параметры акустического сигнала АС
% LS — интегральный уровень АС, измеряемый на расстоянии 1 м
LS=64; % спокойная речь, дБ
% LS7Т — интегральные уровни АС в спектральных полосах при LS=64
LS7T=[47 60 60 55 50 47 43]
LS7=LS7T % здесь - используем теоретические:
% далее могут использоваться уровни LS7,полученные экспериментальным или расчетным путем
% ЗАДАЙТЕ интегральные уровни АС LS7f в спектральных полосах /с фильтров/
%LS7=LS7f
% или
% РАСЧИТАЙТЕ изменения LS при возрастании расстояния r от источника
% r=3
% LS7r=LS7T-20*log10(r)
% LS7=LS7r
% или учитываем изменения LS при прохождении через конструкции
% Z=30 – коэффициент звукоизоляции двери ДСП, дБ
% если расстояния r от источника также задано
|
|
|
% r=3
% LS7rp=LS7T-20*log10(r)-Z
% LS7=LS7rp
% Расчет W - словесной разборчивости речи
% используем справочные данные:
% k — весовые коэффициенты, характеризующие вероятность наличия формант в октавной полосе
k=[0.01 0.03 0.12 0.20 0.30 0.26 0.07];
% A - средний спектральный модальный уровень формант в спектральной полосе, дБ;
A=[25 18 14 9 6 5 4]% из таблицы
% Рассчитаем DА — формантный параметр, характеризующий энергетическую избыточность дискретной составляющей АС
DA=LS7-A
% ОПРЕДЕЛИМ p – матрицу коэффициентов восприятия формант слуховым аппаратом
% человека, которые представляют собой вероятное относительное количество
% формантных составляющих речи, которые будут иметь уровни интенсивности выше
% порогового значения восприятия.
% УЧТЕМ влияние шумов (непреднамеренных либо маскирующих)
% LN(i) - уровень шума (помехи) в i-й частотной полосе, дБ.
% ВВЕСТИ уровни шумов
LN7 = zeros(1,7) %(ПОКА равны 0)
% LN7=[///////] спектральные уровни шума (например, с фильтров LN7f)
q=LS7-LN7 % отношение «уровень речевого сигнала/уровень шума» в i-й частотной полосе, дБ
Q=LS7-DA-LN7
% поэлементно расчет матрицы р:
p = zeros(size(Q))
for i =1:length(Q)
if Q(i)>0
p(i) = 1-(0.78+5.46*exp(-0.0043*(27.3-Q(i)^2)))/(1+(10^0.1*abs(Q(i))));
else p(i) = (0.78+5.46*exp(-0.0043*(27.3-Q(i)^2)))/(1+(10^0.1*abs(Q(i))));
end
end
% как вариант (без учета условия Q<>0) возможен расчет не поэлементно, а матрицы в целом
% pp=1-(0.78+5.46*exp(-0.0043*(27.3-(Q).^2)))./(1+(10^0.1*abs(Q)))
p % вывод результата расчета матрицы р
% расчет R - интегрального индекса артикуляции речи (формантная разборчивость)
% Ri - спектральный индекс понимаемости речи
R7=k.*p
R=sum(R7)
% расчет словесной разборчивости речи W
|
|
|
if R<0.15
W=1.54*R^0.25*[1-exp(-11*R)]
else W=1-exp(-11*R/(1+0.7*R))
end
Результаты расчета:
>> oktava Расчет для r=1 м без учета преднамеренных маскирующих шумов
fc =
125 250 500 1000 2000 4000 8000
LS7T =
47 60 60 55 50 47 43
LS7 =
47 60 60 55 50 47 43
A =
25 18 14 9 6 5 4
DA =
22 42 46 46 44 42 39
LN7 =
0 0 0 0 0 0 0
q =
47 60 60 55 50 47 43
Q =
25 18 14 9 6 5 4
p =
0 0 0 0 0 0 0
p =
-1.2211 0.1405 0.3526 0.3789 0.2461 0.1519 0.0091
R7 =
-0.0122 0.0042 0.0423 0.0758 0.0738 0.0395 0.0006
R =
0.2241
W =
0.8812
% Часть 2. РАСЧЕТ СПЕКТРАЛЬНЫХ УРОВНЕЙ РЕЧЕВОГО СИГНАЛА В ОКТАВНЫХ И ДОЛЕОКТАВНЫХ ПОЛОСАХ
clc; clear; close all;
filename='f:\matlab\02.wav';
[y,Fs] = audioread(filename); % Считывание из файла "*.wav" и запись информации в матрицу y
player = audioplayer(y(1:(1*Fs)), Fs); % Создание объекта типа audioplayer
play(player); % Проигрывание
%создание октавных и долеоктавных фильтров
BandsPerOctave = 3; % 1/3 октавный
N = 6; % Filter Order
F0 = 1000; % Center Frequency (Hz)
% Fs = 32000; % Sampling Frequency (Hz)
f = fdesign.octave(BandsPerOctave,'Class 1','N,F0',N,F0,Fs);
F0 = validfrequencies(f);% можно посмотреть все возможные центральные частоты
F0 = F0(F0 > 90 & F0 < Fs/2 & F0 < 11201); % получилось - 21 фильтр
fm = F0';
Nfc = length(F0);
for i = 1:Nfc,
f.F0 = F0(i);
Hd(i) = design(f,'butter'); %#ok<SAGROW>
end
%
fvtool(Hd(11),'FrequencyScale','log','color','white'); % Анализ фильтра
%11го поддиапазона (1 000 кГц)
hfvt = fvtool(Hd,'FrequencyScale','log','color','white'); % Анализ всей
%группы фильтров
% Nx = 100000;
% SA2 = dsp.SpectrumAnalyzer('SpectralAverages',50,'SampleRate',Fs,...
% 'PlotAsTwoSidedSpectrum',false,'FrequencyScale','Log',...
% 'RBWSource','Property','RBW',2000);
% yw = zeros(Nx,Nfc);
% tic,
% while toc < 15
% % Run for 15 seconds
% xw = randn(Nx,1);
% for i=1:Nfc,
% yw(:,i) = filter(Hd(i),xw);
% end
% step(SA2,yw);
% end
% figure
% plot(xw)
SA = dsp.SpectrumAnalyzer('SpectralAverages',50,'SampleRate',Fs,...
'PlotAsTwoSidedSpectrum',false,'FrequencyScale','Log',...
'YLimits', [-250 -0]);
% YLimits - Границы по оси Y
step1 = Fs/2;
% step1 = length(y);
n_step = floor(length(y)/step1);
Gi = zeros(n_step,Nfc);
Gn = zeros(n_step,Nfc);
G = zeros(n_step,1);
for i = 1:n_step
z = y((step1*(i-1)+1):step1*i);
G(i) = sum(z.^2/length(z)); % Средняя плотность энергии во всей полосе частот
for j = 1:Nfc
reset(Hd(j));
zf = filter(Hd(j),z); % Фильтрация сигнала
Gi(i,j) = sum(zf.^2/length(zf)); % Средняя плотность энергии во одной полосе
% Спектр сигнала и спектр одной полосы
step(SA,[zf; z]');
% axis([0.0125 20 -150 -10]) % не работают для спектроанализатора
% pause(0.1) % можно вставить паузу заданной длины
% pause() % можно вставить ожидание нажатия любой клавиши
reset(SA); % Сброс спектроанализатора
end
i
pause(0.1)
end
for i = 1:n_step
Gn(i,:) = Gi(i,:)/G(i);
end
%--------------------------------------------------------------------------
% Построение распределения по частотам для всех участков звука
figure
h = axes();
surf((1:n_step)*step1/Fs,fm,Gn');
set(h,'Xlim',[1 n_step]*step1/Fs);
set(h,'Ylim',[0 fm(end)]);
set(h,'YScale','Log');
%--------------------------------------------------------------------------
% Построение распределения по частотам для 1 участка звука
figure
h1 = axes();
stem(fm,10*log10(Gn(1,:))) % Вместо 1 можно любой другой № кусочка звука
set(h1,'XScale','Log');
xlabel('Частота, Гц')
ylabel('Уровень в полосе, Дб')
