Matlab實現(xiàn)WAV音頻文件計算聲品質(zhì)參數(shù)：dBA、響度、粗糙度、尖銳度、波動度

這篇具有很好參考價值的文章主要介紹了Matlab實現(xiàn)WAV音頻文件計算聲品質(zhì)參數(shù)：dBA、響度、粗糙度、尖銳度、波動度。希望對大家有所幫助。如果存在錯誤或未考慮完全的地方，請大家不吝賜教，您也可以點擊"舉報違法"按鈕提交疑問。

1.dBA

????????首先讀取WAV文件

[x,Fs] = audioread('pink.wav');   %讀取音頻文件

????????對時域信號進(jìn)行加窗劃分

function [dBA,dBZ,t,windowTime] = analyzeSignal(x,Fs)
responseType = 'fast';
C = 55;
t = 1/Fs:1/Fs:length(x)/Fs;
%% 確定傅里葉窗的大小
if strcmp(responseType,'slow')
   duration = 1.0;
else
   duration = 0.125;
end
N = ceil(duration*Fs);
N = 2^nextpow2(N);
%% 確定信號的dBA
windowStart = [1:N:(length(x)-N)];
dBA = zeros(length(windowStart),1);
dBZ = zeros(length(windowStart),1);
windowTime = t(windowStart+round((N-1)/2));
for i = [1:length(windowStart)]
   [dBA(i),dBZ(i)] = estimateLevel(x(windowStart(i)-1+[1:N]),Fs,C);
end

????????調(diào)用子函數(shù)計算dBA

function [dBA,dBZ] = estimateLevel(x,Fs,C)
X = abs(fft(x));
%% 避免產(chǎn)生對數(shù)取0的情況
X(find(X == 0)) = 1e-17;
%% 保證奈奎斯特采樣定律
f = (Fs/length(X))*[0:(length(X)-1)];
ind = find(f<Fs/2); f = f(ind); X = X(ind);
%% 通過A-weighting濾波器實現(xiàn)dBA計算
A = filterA(f);
XA = A'.*X;
%Z-weighting filter
Z = zeros(length(f),1);
Z(1:end) = 1;
XZ = Z.*X;
%% 用能量計算dBA
totalEnergy = sum(XA.^2)/length(XA);
meanEnergy = totalEnergy/((1/Fs)*length(x));
dBA = 10*log10(meanEnergy)+C;

totalEnergy1 = sum(XZ.^2)/length(XZ);
meanEnergy1 = totalEnergy1/((1/Fs)*length(x));
dBZ = 10*log10(meanEnergy1)+C;

????????A-weighted濾波器子函數(shù)

function A = filterA(f,plotFilter)
%% Define filter coefficients.
c1 = 3.5041384e16;
c2 = 20.598997^2;
c3 = 107.65265^2;
c4 = 737.86223^2;
c5 = 12194.217^2;

%% Evaluate A-weighting filter.
f(find(f == 0)) = 1e-17;
f = f.^2; num = c1*f.^4;
den = ((c2+f).^2) .* (c3+f) .* (c4+f) .* ((c5+f).^2);
A = num./den;

%% 畫dBA計權(quán)圖.
if exist('plotFilter') & plotFilter
    
   % Plot using dB scale.
   figure(2); clf;
   semilogx(sqrt(f),10*log10(A));
   title('A-weighting Filter');
   xlabel('Frequency (Hz)');
   ylabel('Magnitude (dB)');
   xlim([10 100e3]); grid on;
   ylim([-70 10]);
   
   % Plot using linear scale.
   figure(3); plot(sqrt(f),A);
   title('A-weighting Filter');
   xlabel('Frequency (Hz)');
   ylabel('Amplitude');
   xlim([0 44.1e3/2]); grid on;

end

2.響度

????????首先確定窗的大小

function [Loudness,Sharpness,Roughness,Fluctuation] = frameToCalculate(x,Fs);
%% 劃分時域幀的長度
responseType = 'fast';
%避免輸入信號不存在的情況
if ~exist('x')
   [x,Fs] = audioread();
   t = (1/Fs)*[0:(length(x)-1)]; t = t+81;
else
   t = (1/Fs)*[0:(length(x)-1)];
end
%根據(jù)fast和slow進(jìn)行時間幀的劃分
if strcmp(responseType,'slow')
   duration = 1;
else
   duration = 0.125;
end
N = ceil(duration*Fs);
N = 2^nextpow2(N);
% Estimate signal level (within each windowed segment).
windowStart = [1:N:(length(x)-N)];
Loudness = zeros(length(windowStart),1);
Sharpness = zeros(length(windowStart),1);
Roughness = zeros(length(windowStart),1);
Fluctuation = zeros(length(windowStart),1);
windowTime = t(windowStart+round((N-1)/2));
%% 按照劃分的時間幀，分幀計算每幀的平均響度
for i = [1:length(windowStart)]
   Loudness(i) = estimateLoudness(x(windowStart(i)-1+[1:N]),Fs);
   Sharpness(i) = estimateSharpness(x(windowStart(i)-1+[1:N]),Fs);
   Roughness(i) = estimateRoughness(x(windowStart(i)-1+[1:N]),Fs);
   Fluctuation(i) = estimateFluctuation(x(windowStart(i)-1+[1:N]),Fs);
end
end

????????本文采用Zwicker模型，通過24Barks的特征響度進(jìn)行總響度計算：

matlab計算響度,聲品質(zhì),matlab,音頻

? matlab計算響度,聲品質(zhì),matlab,音頻

function loudness = estimateLoudness(x,Fs)
%% 24Barks劃分
fc = [20 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000;
    50 150 250 350 450 570 700 840 1000 1170 1370 1600 1850 2150 2500 2900 3400 4000 4800 5800 7000 8500 10500 13500;
    100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];
BandWidth = [80 100 100 100 110 120 140 150 160 190 210 240 280 320 380 450 550 700 900 1100 1300 1800 2500 3500];
%% 中心頻帶率，有何意義
z = zeros(1,length(fc(2,:)));
for i = 1:length(z)
    z(i) = 13.*atan(0.76*fc(2,i)/1000)+3.5.*(atan(fc(2,i)/1000)).^2; %?
end
%% 聽閾下限聲壓級，用getData在聽閾曲線上估計出來
Lspl = [48.72457794 28.15401982 16.78690581 10.97541703 7.18374514 3.645302897 1.116392826 -0.152013044 0.093691208 -0.165558761 -0.9278816 -1.94268155 -2.95710523 -4.475730591 -6.751034742 -7.765458422 -8.528157532 -8.030352439 -5.516869434 -1.490781387 4.047911702 6.813119278 8.569170952 11.83632259];
%% Calculate magnitude of FFT.
X = abs(fft(x));
X(find(X == 0)) = 1e-17;%避免產(chǎn)生對數(shù)取0的情況
% Retain frequencies below Nyquist rate.
f = (Fs/length(X))*[0:(length(X)-1)];
ind = find(f<Fs/2); f = f(ind); X = X(ind);
%% 按bark進(jìn)行聲壓級計算
SPL = zeros(1,24);
for i = 1:24
    location = find(f>=fc(1,i)&f<=fc(3,i));
    if ~isempty(location)
        totalEnergy = sum(X(location).^2)/length(location);
        SPL(i) = 10*log10(totalEnergy/((1/Fs)*length(x)))+55;
    else
        SPL(i) = 0;
    end
end
%% 按Bark進(jìn)行響度計算
N = zeros(1,length(Lspl));
N = 0.08.*((10.^(Lspl/10)).^0.23).*((0.5+0.5.*10.^((SPL-Lspl)/10)).^0.23-1);%得到24bark響度分布圖
loudness = sum(N);
end

3.尖銳度Sharpness

????????同樣選擇Zwicker模型

matlab計算響度,聲品質(zhì),matlab,音頻

function sharpness = estimateSharpness(x,Fs)
%% 24Barks劃分
fc = [20 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000;
    50 150 250 350 450 570 700 840 1000 1170 1370 1600 1850 2150 2500 2900 3400 4000 4800 5800 7000 8500 10500 13500;
    100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];
BandWidth = [80 100 100 100 110 120 140 150 160 190 210 240 280 320 380 450 550 700 900 1100 1300 1800 2500 3500];
%% 中心頻帶率，有何意義？
z = zeros(1,length(fc(2,:)));
for i = 1:length(z)
    z(i) = 13.*atan(0.76*fc(2,i)/1000)+3.5.*(atan(fc(2,i)/1000)).^2; %?
end
%% 根據(jù)中心頻帶率z加上響度計權(quán)函數(shù)
gz = zeros(1,length(z));
for i = 1:length(z)
    if z(i)<=16
        gz(i) = z(i);
    else
        gz(i) = 0.06.*exp(0.17*z(i));
    end
end
%% 聽閾下限聲壓級，用getData在聽閾曲線上估計出來
Lspl = [48.72457794 28.15401982 16.78690581 10.97541703 7.18374514 3.645302897 1.116392826 -0.152013044 0.093691208 -0.165558761 -0.9278816 -1.94268155 -2.95710523 -4.475730591 -6.751034742 -7.765458422 -8.528157532 -8.030352439 -5.516869434 -1.490781387 4.047911702 6.813119278 8.569170952 11.83632259];
%% Calculate magnitude of FFT.
X = abs(fft(x));
X(find(X == 0)) = 1e-17;%避免產(chǎn)生對數(shù)取0的情況
% Retain frequencies below Nyquist rate.
f = (Fs/length(X))*[0:(length(X)-1)];
ind = find(f<Fs/2); f = f(ind); X = X(ind);
%% 按bark進(jìn)行聲壓級計算
SPL = zeros(1,24);
for i = 1:24
    location = find(f>=fc(1,i)&f<=fc(3,i));
    if ~isempty(location)
        totalEnergy = sum(X(location).^2)/length(location);
        SPL(i) = 10*log10(totalEnergy/((1/Fs)*length(x)))+55;
    else
        SPL(i) = 0;
    end
end
%% 按Bark進(jìn)行響度計算
N = zeros(1,length(Lspl));
N = 0.08.*((10.^(Lspl/10)).^0.23).*((0.5+0.5.*10.^((SPL-Lspl)/10)).^0.23-1);%得到24bark響度分布圖
%% 計算sharpness
sharpness = 0.104.*sum(N.*gz.*z)./sum(N);
end

4.粗糙度

????????粗糙度計算公式

matlab計算響度,聲品質(zhì),matlab,音頻

?由于掩蔽深度無法用數(shù)學(xué)公式進(jìn)行描述，在這里用一個系數(shù)乘以響度簡單代替（希望各位能教教我更準(zhǔn)確的數(shù)值計算方法）

function roughness = estimateRoughness(x,Fs)
%% 24Barks劃分
fc = [20 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000;
    50 150 250 350 450 570 700 840 1000 1170 1370 1600 1850 2150 2500 2900 3400 4000 4800 5800 7000 8500 10500 13500;
    100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];
BandWidth = [80 100 100 100 110 120 140 150 160 190 210 240 280 320 380 450 550 700 900 1100 1300 1800 2500 3500];
%% 中心頻帶率
z = zeros(1,length(fc(2,:)));
for i = 1:length(z)
    z(i) = 13.*atan(0.76*fc(2,i)/1000)+3.5.*(atan(fc(2,i)/1000)).^2; %
end
%% 聽閾下限聲壓級，用getData在聽閾曲線上估計出來
Lspl = [48.72457794 28.15401982 16.78690581 10.97541703 7.18374514 3.645302897 1.116392826 -0.152013044 0.093691208 -0.165558761 -0.9278816 -1.94268155 -2.95710523 -4.475730591 -6.751034742 -7.765458422 -8.528157532 -8.030352439 -5.516869434 -1.490781387 4.047911702 6.813119278 8.569170952 11.83632259];
%% Calculate magnitude of FFT.
X = abs(fft(x));
X(find(X == 0)) = 1e-17;%避免產(chǎn)生對數(shù)取0的情況
% Retain frequencies below Nyquist rate.
f = (Fs/length(X))*[0:(length(X)-1)];
ind = find(f<Fs/2); f = f(ind); X = X(ind);
%% 按bark進(jìn)行聲壓級計算
SPL = zeros(1,24);
for i = 1:24
    location = find(f>=fc(1,i)&f<=fc(3,i));
    if ~isempty(location)
        totalEnergy = sum(X(location).^2)/length(location);
        SPL(i) = 10*log10(totalEnergy/((1/Fs)*length(x)))+55;
    else
        SPL(i) = 0;
    end
end
%% 按Bark進(jìn)行粗糙度計算
N = zeros(1,length(Lspl));
N = 0.08.*((10.^(Lspl/10)).^0.23).*((0.5+0.5.*10.^((SPL-Lspl)/10)).^0.23-1);%得到24bark響度分布圖
fmod = 300/1000;  %調(diào)制聲頻率
k = 0.01;   %掩蔽深度轉(zhuǎn)換為響度的系數(shù)
L = k.*N;   %掩蔽深度
roughness = 0.3.*fmod.*sum(L);
end

5.波動度

根據(jù)計算模型

matlab計算響度,聲品質(zhì),matlab,音頻

function fluctuation = estimateFluctuation(x,Fs)
%% 24Barks劃分
fc = [20 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000;
    50 150 250 350 450 570 700 840 1000 1170 1370 1600 1850 2150 2500 2900 3400 4000 4800 5800 7000 8500 10500 13500;
    100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];
BandWidth = [80 100 100 100 110 120 140 150 160 190 210 240 280 320 380 450 550 700 900 1100 1300 1800 2500 3500];
%% 中心頻帶率
z = zeros(1,length(fc(2,:)));
for i = 1:length(z)
    z(i) = 13.*atan(0.76*fc(2,i)/1000)+3.5.*(atan(fc(2,i)/1000)).^2; %?
end
%% 聽閾下限聲壓級，用getData在聽閾曲線上估計出來
Lspl = [48.72457794 28.15401982 16.78690581 10.97541703 7.18374514 3.645302897 1.116392826 -0.152013044 0.093691208 -0.165558761 -0.9278816 -1.94268155 -2.95710523 -4.475730591 -6.751034742 -7.765458422 -8.528157532 -8.030352439 -5.516869434 -1.490781387 4.047911702 6.813119278 8.569170952 11.83632259];
%% Calculate magnitude of FFT.
dBZ(dBZ == 0) = 1e-17;%避免產(chǎn)生對數(shù)取0的情況
%% 按bark進(jìn)行聲壓級計算
SPL = zeros(1,24);
for i = 1:24
    location = find(f>=fc(1,i)&f<=fc(3,i));
    if ~isempty(location)
        SPL(i) = sum(X(location))/length(location);
    else
        SPL(i) = 0;
    end
end
%% 按Bark進(jìn)行響度計算
N = zeros(1,length(Lspl));
N = 0.08.*((10.^(Lspl/10)).^0.23).*((0.5+0.5.*10.^((SPL-Lspl)/10)).^0.23-1);%得到24bark響度分布圖
k = 0.01;   %掩蔽深度轉(zhuǎn)換為響度的系數(shù)
L = k.*N;   %掩蔽深度
fmod = 300/1000;  %調(diào)制頻率
fluctuation = 0.008.*sum(L)/((fmod/4)+(4/fmod));
end

（對于掩蔽深度的計算暫未找到合適的數(shù)學(xué)模型）文章來源地址http://www.zghlxwxcb.cn/news/detail-810487.html

到了這里，關(guān)于Matlab實現(xiàn)WAV音頻文件計算聲品質(zhì)參數(shù)：dBA、響度、粗糙度、尖銳度、波動度的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！