function [ STI_val, ALcons, STI_val_approx, ALcons_approx ] = STI( ImpulseResponse, fs, OctaveFilters ) % Calculation of the Speech Transmission Index (STI) % % Syntax: [STI_val, ALcons, STI_val_approx, ALcons_approx] = STI(ImpulseResponse, fs, OctaveFilters) % This function calculates the Speech Transmission Index from a given % impulse response. Octave band filters (125Hz-8kHz) can be passed to % the function as an array for faster computation. % % Inputs: % ImpulseResponse - Impulse response of the channel to test % fs - Sampling frequency of the impulse response % OctaveFilters - Array of MATLAB octave band filters % % Outputs: % STI_val - The Speech Transmission Index result % ALcons - Articulation loss of consonants result % STI_val_approx - STI approximation from computable bands if fs too low % ALcons_approx - ALcons approximation from computable bands if fs too low % % Example: % y=sinc(-7999:8000); % fs=44100; % [STIval,ALcons,STIval_,ALcons_]=STI(y,fs) % fs=16000; % [STIval,ALcons,STIval_,ALcons_]=STI(y,fs) % fs=8000; % [STIval,ALcons,STIval_,ALcons_]=STI(y,fs) % % See also: STI_BandFilters.m, extractIR.m, synthSweep.m % Author: Jacob Donley % University of Wollongong % Email: jrd089@uowmail.edu.au % Version: 1.0 (12 April 2017) % Version: 0.1 (30 September 2015) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 1. Take the full bandwidth impulse response, run it through octave band filters BandsPerOctave = 1; N = 6; % Filter Order F0 = 1000; % Center Frequency (Hz) f = fdesign.octave(BandsPerOctave,'Class 1','N,F0',N,F0,fs); F0 = validfrequencies(f); F0(F0<125 | F0>min([8000,fs/2]))=[]; % Only keep bands in the range required for the STI calculation Nfc = length(F0); % Create MTF function handle getMTF = @(ir) ... abs( ... Take the square root of the sum of the real part squared and the imaginary part squared fft(ir .^ 2) ... Compute the Fourier transform of the squared impulse response / sum(ir .^ 2) );% Integrate the squared impulse response to get the total energy and normalize the envelope spectrum for i=1:Nfc if nargin < 3 f.F0 = F0(i); H = design(f,'butter'); ir_filtered = filter(H, ImpulseResponse); else ir_filtered = filter(OctaveFilters(i), ImpulseResponse); end MTF_octband(:,i) = getMTF(ir_filtered); end % 2. Take the amplitude at 14 modulation frequencies %modulation_freqs = [ 0.63, 0.80, 1.00, 1.25, 1.60, 2.00, 2.50, 3.15, 4.00, 5.00, 6.30, 8.00, 10.00, 12.50]; modulation_freqs = 2 .^ ([-2:11]./3); freqs = linspace(0, fs/2, size(MTF_octband,1)/2+1); freqs(end)=[];%No nyquist frequency for i=1:Nfc m(i,:) = interp1(freqs,MTF_octband(1:end/2,i),modulation_freqs); end good_freqs = ~any(isnan(m),2); m(~good_freqs,:)=[]; % 3. Convert each of the 98 m values into an “apparent signal-to-noise ratio” (S/N) in dB SNR_apparent = pow2db( m ./ (1-m) ); % 4. Limit the Range SNR_apparent( SNR_apparent > 15 ) = 15; SNR_apparent( SNR_apparent < -15 ) = -15; % 5. Compute the mean (S/N) for each octave band SNR_avg = mean(SNR_apparent,2); % 6. Weight the octave mean (S/N) values W = [0.13, 0.14, 0.11, 0.12, 0.19, 0.17, 0.14]; % Values from the STI standard SNR_Wavg = sum(SNR_avg' .* W(good_freqs)); SNR_Wavg_approx = sum(SNR_avg' .* W(good_freqs)/sum(W(good_freqs))); % 7. Convert the overall mean (S/N) to an STI value STI_val = (SNR_Wavg + 15) / 30; STI_val_approx = (SNR_Wavg_approx + 15) / 30; % (Optional) Convert STI to ALcons sti2alcons = @(x) 170.5405*exp(-5.419*x); ALcons = sti2alcons(STI_val); ALcons_approx = sti2alcons(STI_val_approx); end