Skip to content

Commit 436b453

Browse files
Andronicus1000Andronicus1000
authored
Files
1 parent a3cd31e commit 436b453

22 files changed

+3278
-0
lines changed

‎AddReceiver.m

+472
Large diffs are not rendered by default.

‎AddSource.m

+557
Large diffs are not rendered by default.

‎AirAbsorption.m

+56
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
function A = AirAbsorption(f,T,H,P)
2+
% AIRABSORPTION
3+
% Computes the frequency dependent air absorption coefficients in dB/m.
4+
% Uses ISO 9613-1:1993 standard.
5+
%
6+
% A = AirAbsorption(f)
7+
% A = AirAbsorption(f,T)
8+
% A = AirAbsorption(f,T,H)
9+
% A = AirAbsorption(f,T,H,P)
10+
%
11+
% where A - Frequency dependent air absorption (dB/m)
12+
% f - Vector of frequencies (Hz)
13+
% T - Temperature (Celcius) [Default: 20]
14+
% H - Relative humidity (0 to 1) [Default: 0.42]
15+
% P - Atmospheric pressure (kPa) [Default: 101.325]
16+
%
17+
% A. Wabnitz, N. Epain, 2011
18+
19+
% Set defaults
20+
if nargin < 4
21+
P = 101.325 ;
22+
if nargin < 3
23+
H = 0.42 ;
24+
if nargin < 2
25+
T = 20 ;
26+
end
27+
end
28+
end
29+
30+
% Constants
31+
P0 = 101.325 ; % Reference pressure (mean sea level, in kPa)
32+
T0 = 293.15 ; % Reference absolute temperature (in Kelvin)
33+
34+
% Absolute temperature
35+
Ta = T + 273.15 ;
36+
37+
% Relative pressure and temperature
38+
Tr = Ta / T0 ;
39+
Pr = P / P0 ;
40+
41+
% Molar concentration of water vapour (%)
42+
h = H * 100 * 10^(4.6151 - 6.8346 * (273.16/Ta)^1.261) / Pr ;
43+
44+
% Oxygen relaxation frequency
45+
frO = Pr * (24 + 4.04e4 * h * (0.02 + h)/(0.391 + h));
46+
47+
% Nitrogen relaxation frequency
48+
frN = Pr * Tr^(-1/2) * (9 + 280 * h * exp(-4.17*(Tr^(-1/3)-1)));
49+
50+
% Calculate the air absorption
51+
f2 = f.^2 ;
52+
A = 8.686 * f2 .* ( 1.84e-11 * (1/Pr) * sqrt(Tr) ...
53+
+ Tr^(-5/2) * ( 0.01275 * exp(-2239.1/Ta) ./ (frO+f2/frO) ...
54+
+ 0.1068 * exp(-3352/Ta) ./ (frN+f2/frN) ) ) ;
55+
56+

‎BandPassFilter.m

+103
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
function SigOut = BandPassFilter(SigInp,frq,smpFrq,bndOpt)
2+
%
3+
% SigOut = BandPassFilter(SigInp,frq,smpFrq,bndOpt) ;
4+
%
5+
% Band-pass filter a set of signals with octave-band or third-octave band
6+
% filters.
7+
%
8+
% Input: - SigInp is a set of signals ([nbSamp x nbChan] array)
9+
% - frq is the vector of the band-pass filter center frequencies.
10+
% - smpFrq is the sampling frequency (default value is 48000)
11+
%
12+
% Output: SigOut is the array of the band-pass filtered signals.
13+
% If nbChan = 1, SigOut is a [nbSamp x nbFreq] matrix.
14+
% Otherwise, SigOut is a [nbSamp x nbChan x nbFreq] 3D array.
15+
%
16+
% Options: bndOpt can be set to 'oct' (default) or '3rdoct' to filter with
17+
% octave-band or third-octave band filters.
18+
%
19+
% Note: Requires MATLAB's Signal Processing Toolbox
20+
%
21+
% N.Epain, 2011
22+
23+
% Default mode: 'oct' for octave band filtering
24+
if nargin < 4
25+
bndOpt = 'oct' ;
26+
end
27+
bndOpt = lower(bndOpt) ;
28+
29+
% Default sampling frequency: 48kHz
30+
if (nargin<3) || isempty(smpFrq)
31+
smpFrq = 48e3 ;
32+
end
33+
34+
% Min and max frequencies at which the order 6 and 8 Butterworth filters
35+
% can be used
36+
switch bndOpt
37+
case 'oct'
38+
minFrq = smpFrq/200 ;
39+
maxFrq = smpFrq/8 ;
40+
case '3rdoct'
41+
minFrq = smpFrq/80 ;
42+
maxFrq = smpFrq/8 ;
43+
end
44+
45+
% Number of frequencies
46+
nmbFrq = length(frq) ;
47+
48+
% Number of audio channels
49+
nmbChn = size(SigInp,2) ;
50+
51+
% Signal length
52+
sigLng = size(SigInp,1) ;
53+
54+
% Band-pass filter the signals
55+
SigOut = zeros(sigLng,nmbChn,nmbFrq) ;
56+
for I = 1 : nmbFrq
57+
58+
% Check if the current center frequency is in the minFrq -> maxFrq
59+
% interval. If not, the signal will be resample.
60+
if frq(I) < minFrq
61+
% The frequency is too low, the signal will be downsampled
62+
rat = [1 2^nextpow2(minFrq/frq(I))] ;
63+
elseif frq(I) > maxFrq
64+
% The frequency is too high, the signal will be upsampled
65+
rat = [2^nextpow2(frq(I)/maxFrq) 1] ;
66+
else
67+
% Keep the current smaple frequency
68+
rat = [1 1] ;
69+
end
70+
71+
% Band-pass filter for the current frequency
72+
switch bndOpt
73+
case 'oct'
74+
[num,den] = ...
75+
butter(3,[2^(-1/2) 2^(1/2)]*frq(I)/smpFrq*2*rat(2)/rat(1));
76+
case '3rdoct'
77+
[num,den] = ...
78+
butter(4,[2^(-1/6) 2^(1/6)]*frq(I)/smpFrq*2*rat(2)/rat(1));
79+
end
80+
81+
% Filter the signals
82+
for J = 1 : nmbChn
83+
84+
% Current signal
85+
sig = SigInp(:,J) ;
86+
87+
% Resample the signal
88+
sig = resample(sig,rat(1),rat(2)) ;
89+
90+
% Filter it
91+
sig = filter(num,den,sig) ;
92+
93+
% Resample back
94+
sig = resample(sig,rat(2),rat(1)) ;
95+
SigOut(:,J,I) = sig(1:sigLng) ;
96+
97+
end
98+
99+
end
100+
101+
% Squeeze the output array
102+
% (the output is a matrix if the input signal is one channel only)
103+
SigOut = squeeze(SigOut) ;

‎Clap.m

+44
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
function snd = Clap(ImpRsp,smpFrq)
2+
%
3+
% snd = Clap(ImpRsp,smpFrq) ;
4+
%
5+
% Play a clap convolved with an impulse response.
6+
%
7+
% Input: - ImpRsp is an impulse reponse ([Nx1] vector) or array of impulse
8+
% responses ([NxM] matrix).
9+
% - smpFrq is the sampling frequency.
10+
% This parameter can be omitted, the default value is 48kHz.
11+
%
12+
% Output: - if ImpRsp is a vector, snd is a vector whose elements are a
13+
% clap sound convolved with the input impulse response. If ImpRsp
14+
% is a matrix then snd is the convolution of the clap with the 2
15+
% first columns of ImpRsp.
16+
%
17+
%
18+
% N.Epain, 2011
19+
20+
% Default sampling frequency
21+
if nargin < 2
22+
smpFrq = 48000 ;
23+
end
24+
25+
% Load the clap
26+
load MCRoomSim_clap.mat
27+
28+
% Resample the clap if the sampling frequency is not 48kHz
29+
if smpFrq ~= 48000
30+
clp = resample(clp,smpFrq,48000) ;
31+
end
32+
33+
% Convolve the clap with the impulse response(s)
34+
if size(ImpRsp,2) == 1
35+
snd = fftfilt(ImpRsp,[clp;zeros(size(ImpRsp,1),1)]) ;
36+
else
37+
snd = fftfilt(ImpRsp(:,1:2),[clp;zeros(size(ImpRsp,1),1)]) ;
38+
end
39+
40+
% Normalise the sound peak level
41+
snd = .99*snd/max(max(abs(snd))) ;
42+
43+
% Play the sound
44+
sound(snd,smpFrq) ;

‎EarlyDecayTime.m

+153
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
function [EDT,frq] = EarlyDecayTime(ImpRsp,smpFrq,bndOpt,pltEDT,pltSchCur)
2+
%
3+
% [EDT,frq] = EarlyDecayTime(ImpRsp,smpFrq,bndOpt,pltRvbTme,pltSchCur) ;
4+
%
5+
% Calculate the Early Decay Times from a set of impulse responses.
6+
%
7+
% Input: - ImpRsp is an impulse reponse ([Nx1] vector) or array of impulse
8+
% responses ([NxM] matrix).
9+
% - smpFrq is the sampling frequency.
10+
% This parameter can be omitted, the default value is 48kHz.
11+
%
12+
% Output: - if ImpRsp is a vector, EDT is the vector of the EDTs for
13+
% different frequency bands. If ImpRsp is a matrix, EDT is the
14+
% matrix of the EDTs for each impulse response.
15+
% - frq is the vector of the frequency band center frequencies.
16+
%
17+
% Options: - bndOpt can be set to 'oct' (default) to obtain the EDTs for 9
18+
% octave bands with center frequency from 62.5 to 16kHz, or
19+
% '3rdoct' in which case the EDTs are calculated for 25
20+
% third-octave bands.
21+
% - if pltEDT is set to true (default), the EDTs are plotted.
22+
% - if pltSchCur is set to true, the Schroeder curves are plotted.
23+
%
24+
% Note: Requires MATLAB's Signal Processing Toolbox
25+
%
26+
% N.Epain, 2011
27+
28+
% Don't plot the Schroeder curve by default
29+
if nargin < 5
30+
pltSchCur = false ;
31+
end
32+
33+
% Plot the EDTs by default
34+
if nargin < 4
35+
pltEDT = true ;
36+
end
37+
38+
% Default "frequency band option": 'oct'
39+
% (octave bands with center frequencies from 62.5 to 16 kHz)
40+
if nargin < 3
41+
bndOpt = 'oct' ;
42+
end
43+
bndOpt = lower(bndOpt) ;
44+
45+
% Default sampling frequency: 48kHz
46+
if (nargin<2) || isempty(smpFrq)
47+
smpFrq = 48e3 ;
48+
end
49+
50+
% Length of the impulse responses
51+
nmbSmp = size(ImpRsp,1) ;
52+
53+
% Number of impulse responses
54+
nmbImp = size(ImpRsp,2) ;
55+
56+
% Vector of the center frequencies
57+
switch bndOpt
58+
case 'oct'
59+
frq = 1000 * 2.^(-4:4)' ;
60+
case '3rdoct'
61+
frq = 1000 * 2.^(-4:1/3:4)' ;
62+
end
63+
nmbFrq = length(frq) ;
64+
65+
% Colors used for the reverb time and Schroeder curve plots
66+
clr = lines(nmbImp) ;
67+
68+
% Initialise the output
69+
EDT = zeros(nmbFrq,nmbImp) ;
70+
71+
% Create a new figure for the Schroeder curves
72+
if pltSchCur == true
73+
figure
74+
end
75+
76+
% Calculate the Schroeder curves
77+
SchCur = SchroederCurve(ImpRsp,smpFrq,bndOpt,false) ;
78+
if nmbImp == 1
79+
SchCur = permute(SchCur,[1 3 2]) ;
80+
end
81+
82+
% Loop on the impulse responses
83+
for J = 1 : nmbImp
84+
85+
% Loop on the center frequencies
86+
for I = 1 : nmbFrq
87+
88+
% Estimate the slope of the Schroeder curve between -5 and -35dB
89+
fst = find(SchCur(:,J,I)>=-1e-2,1,'last') ;
90+
lst = find(SchCur(:,J,I)<=-10,1,'first') ;
91+
coe = LinearRegression((fst:lst)'/smpFrq,SchCur(fst:lst,J,I)) ;
92+
93+
% Plot the Shroeder curve
94+
if pltSchCur == true
95+
switch bndOpt
96+
case 'oct'
97+
subplot(3,3,I), hold on
98+
case '3rdoct'
99+
subplot(5,5,I), hold on
100+
end
101+
plot((1:nmbSmp)/(smpFrq),SchCur(:,J,I), ...
102+
'color',clr(J,:),'linewidth',2)
103+
plot((1:nmbSmp)/(smpFrq), ...
104+
coe(2)+coe(1)*(1:nmbSmp)/(smpFrq), ...
105+
'color',clr(J,:))
106+
title(['f = ' num2str(round(frq(I))) ' Hz']) ;
107+
xlabel('Time [s]')
108+
ylabel('Energy [dB]')
109+
axis([0 1.05*nmbSmp/smpFrq -95 5])
110+
grid on, box on
111+
end
112+
113+
% EDT
114+
EDT(I,J) = -60/coe(1) ;
115+
116+
end
117+
118+
end
119+
120+
% Plot the EDTs
121+
if pltEDT == true
122+
figure('color','white')
123+
semilogx(frq,EDT,'-o','linewidth',2)
124+
xlim([frq(1)*2^(-1/3) frq(end)*2^(1/3)])
125+
title('Early Decay Time (EDT)','fontsize',14) ;
126+
xlabel('Frequency [Hz]','fontsize',14)
127+
ylabel('EDT [s]','fontsize',14)
128+
set(gca,'YGrid','on')
129+
set(gca,'Xtick',[],'fontsize',14)
130+
switch lower(bndOpt)
131+
case 'oct'
132+
set(gca,'Xtick',[frq(1)/2;frq;frq(end)*2])
133+
case '3rdoct'
134+
set(gca,'Xtick',[frq(1)/2;frq(1:3:end);frq(end)*2])
135+
end
136+
box on
137+
end
138+
139+
140+
% Linear regression subroutine
141+
function coe = LinearRegression(tme,eng)
142+
143+
% Number of points to fit
144+
nmb = length(tme) ;
145+
146+
% Matrix we need to invert
147+
Mat = [tme(:) ones(nmb,1)] ;
148+
149+
% Linear coefficients
150+
coe = pinv(Mat) * eng ;
151+
152+
153+

0 commit comments

Comments
 (0)