Modulation Spectrum (MATLAB batch file)
Put all the files that require processing together in one folder, then run the following function. It will prompt you for the folder. Results (24 x 5 ) will be in the MATLAB workspace.
function [All_files,Total_rms,diff_rms,tA,CF_MFB] = Multifile_24x5_modulation_filterbank_SAF_FINAL(FSamp)
% DESCRIPTION :
% **** BEFORE YOU START PROGRAM ****
% Specify the sampling rate of the sound files you are going to process
% by typing at the command line: FSamp == 44100; or FSamp == 16000;
% if sampling rate is 48000 then down sample files by 3 first.
% **********************************************************
% Produces downsampled 5x24 spectro-temporal envelopes of audio (.wav) files. The 5 spectral
% channels span 100-7250 Hz, the 24 modulation rate bands span 0.9-40 Hz .
% INPUT VARIABLES :
% FSamp = Sampling frequency of wav files, in Hz. Only 44100 or 16000 accepted.
% OUTPUT VARIABLES :
% All_files is an array with each file's data separate (number of files) x 24 x 5
%
% CF = Arithmetic centre frequencies of the 5 spectral channels
% CF_MFB = Arithmetic centre frequencies of the 3 modulation rate channels (N.B. It may be more appropriate to report geometric CFs of 1.5 Hz, 5.5 Hz & 21.9 Hz)
% NSamp = Final downsampled frequency of the spectro-temporal envelopes (1050 or 1000 Hz)
% Total_rms is the array 1 x 24 x 5 of the mean of all the files.
%%% Analyse all of the signal, envelope is downsampled to 1050/1000 Hz.
% To use define the sampling rate i.e. FSamp = 44100
% then run filter by copying this to command window:
% All_files,Total_rms,diff_rms,tA,CF_MFB] = Multifile_24x5_modulation_filterbank(FSamp)
%%%%%%%%%%*******************************************************************************%%%%%
if FSamp == 44100;
ds1 = 3; %% if FSamp = 44.1 kHz, iSamp = 14700
ds2 = 14; %% NSamp = 1050 Hz
elseif FSamp == 16000;
ds1 = 1; %% iSamp = 16000
ds2 = 16; %% NSamp = 1000;
else
display('Incompatible sampling rate, please use 44100 or 16000 only')
end
iSamp = FSamp/ds1; %% Interim sampling frequency for spectral filtering
NSamp = iSamp/ds2;
N = 5;
edges = [10;100;300;700;1750;3900;7250]; %% Include first LPF dummy channel
[CF, bpfs] = MFB_coeffs(edges,iSamp,1);
CF = CF(2:end);
M = 24;
Fcor = [0.79; 0.93; 1.09; 1.27; 1.49; 1.74; 2.03; 2.38; 2.78; 3.25; 3.80; 4.45;...
5.20; 6.08; 7.11; 8.32; 9.72; 11.38; 13.30; 15.56; 18.20; 21.28; 24.89; 29.11;...
34.04; 39.81];
[CF_MFB,MFB_bpfs] = MFB_coeffs(Fcor,NSamp,2); %% Filter length ~7000
sil = zeros(length(MFB_bpfs)/2,1); %% Silence of half filter length to pad each side
CF_MFB = CF_MFB(2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
myDir = uigetdir; %gets directory
myFiles = dir(fullfile(myDir,'*.wav')); %gets all wav files in struct
L = length(myFiles);
All_files = zeros(L,24,5);
for k = 1:length(myFiles)
baseFileName = myFiles(k).name;
fullFileName = fullfile(myDir, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
[y,FSamp] = audioread(fullFileName);
y = y(1:ds1:end); %% Downsample sound once
%%%% 5 SPECTRAL BANDS
yfil = zeros(length(y),N);
for n = 2:N+1 %% Ignore first LPF channel
%%% Band-pass filters to form channel envelope signals
l_fil = bpfs(1,n); %% length of filter
fil_n = bpfs(2:1+l_fil,n); %% extract filter
nshift = floor(l_fil/2); %% sample shift to correct for FIR BPF delay
yfil = filter(fil_n, 1, y);
yfil = sample_advance(yfil, nshift, 1e-7); %% time-advance signal, preserve buffer length and pad with externally defined value
% ychan(:,n-1) = single(yfil(1:2:end)); %% Downsample by half and convert to single point precision
%%% Hilbert extraction
Ahil = abs(hilbert(yfil)); %% Hilbert envelope
Ahil = Ahil(1:ds2:end); %% Downsample envelope to NSamp (already downsmples by 2 in ychan)
clear yfil
%%%% 24 MODULATION RATE BANDS
Ahil = [sil; Ahil; sil];
Atmp = zeros(length(Ahil),1);
Afil = zeros(length(Ahil),M);
for m = 2:M+1
bpf_len = MFB_bpfs(1,m); % extract filter length from array
bpf_chan = MFB_bpfs(2:bpf_len+1,m); % extract filter
Atmp = filter(bpf_chan, 1, Ahil); % apply to hilbert envelope using MATLAB filter() function
dly_shift = floor(bpf_len/2); % compensate shift to time-align FIR filter outputs
len_A = length(Atmp) - dly_shift; % advance filter outputs
Afil(1:len_A,m-1) = Atmp(1+dly_shift:length(Atmp)); % time advance to preserve initial bit
Afil(1+len_A:length(Atmp),m-1) = 0; % kill the rest, since it contains nothing useful
clear Atmp
Atmp = zeros(length(Ahil),1);
end
A(:,:,n-1) = Afil(length(sil)+1:end-length(sil),:); % Crop silence
end
%Y = ychan;
tA = A;
clear A;
clear y;
clear Afil;
%%%%%%%%%%%%%%%%% April 16. S.F. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rmstA = rms(tA); % RMS of each modulation Channel (24 X 5))
overall_rms = rms(rmstA); % overall RMS of each spectral band
d = 5;
for S = 1:d
diff_rms(:,:,S) = (rmstA(:,:,S)/overall_rms(:,:,S)); % normalise data
end
clear S
All_files(k,:,:) = diff_rms(:,:,:);
end
intermediate = (mean(All_files,1)); % average the bands across files to give 1x24x5 output
Total_wavs_rms = 20*log10(intermediate);
big_intermediate = mean(intermediate,3); % Average across filesand bands to give 24x1 output
Total_rms = 20*log10(big_intermediate);
figure();
%semilogx(CF_MFB, Total_rms); grid on;
semilogx(CF_MFB,Total_wavs_rms(:,:,1),CF_MFB,Total_wavs_rms(:,:,2),CF_MFB,Total_wavs_rms(:,:,3),CF_MFB,Total_wavs_rms(:,:,4),CF_MFB,Total_wavs_rms(:,:,5),CF_MFB,Total_rms,'K-*')
grid on;
legend('Band 1:100-300 Hz','Band 2:300-700 Hz','Band 3:700 - 1750 Hz','Band 4:1750 - 3900 Hz ','Band 5:3900-7250 Hz','Average all bands');
xlabel('Modulation Rate (Hz)');
ylabel('Power(dB from mean) ');
%clear Afil
clear Atmp
%clear A
%clear ychan