rmeddis@38: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rmeddis@38: % This program is free software; you can redistribute it and/or modify rmeddis@38: % it under the terms of the GNU General Public License as published by rmeddis@38: % the Free Software Foundation; either version 2 of the License, or rmeddis@38: % (at your option) any later version. rmeddis@38: % rmeddis@38: % This program is distributed in the hope that it will be useful, rmeddis@38: % but WITHOUT ANY WARRANTY; without even the implied warranty of rmeddis@38: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the rmeddis@38: % GNU General Public License for more details. rmeddis@38: % rmeddis@38: % You can obtain a copy of the GNU General Public License from rmeddis@38: % http://www.gnu.org/copyleft/gpl.html or by writing to rmeddis@38: % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. rmeddis@38: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rmeddis@38: rmeddis@38: classdef cJob rmeddis@38: %CJOB Responsible for conversion of wav files into recogniser features rmeddis@38: % Please see the documentation located in a separate file for further rmeddis@38: % information. rmeddis@38: rmeddis@38: %% ********************************************************* rmeddis@38: % properties _ _ rmeddis@38: % | | (_) rmeddis@38: % _ __ _ __ ___ _ __ ___ _ __| |_ _ ___ ___ rmeddis@38: % | '_ \| '__/ _ \| '_ \ / _ \ '__| __| |/ _ \/ __| rmeddis@38: % | |_) | | | (_) | |_) | __/ | | |_| | __/\__ \ rmeddis@38: % | .__/|_| \___/| .__/ \___|_| \__|_|\___||___/ rmeddis@38: % | | | | rmeddis@38: % |_| |_| rmeddis@38: %************************************************************ rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % Public properties - can be set by user rmeddis@38: %************************************************************ rmeddis@38: properties(Access = public) rmeddis@38: wavFolder rmeddis@38: opFolder rmeddis@38: noiseFolder rmeddis@38: rmeddis@38: wavList rmeddis@38: todoStatus rmeddis@38: rmeddis@38: rmeddis@38: participant = 'Normal';%'DEMO2_multiSpont'; rmeddis@38: noiseName = '8TalkerBabble'; rmeddis@38: numWavs = 5; rmeddis@38: noiseLevToUse = -200; rmeddis@38: speechLevToUse = 50; rmeddis@38: rmeddis@38: speechLevStd = 0; rmeddis@38: noiseLevStd = 0; rmeddis@38: freezeNoise = 0; rmeddis@38: meanSNR = 20; rmeddis@38: speechDist = 'None'; rmeddis@38: noiseDist = 'None'; rmeddis@38: rmeddis@38: noisePreDur = 0.0; rmeddis@38: noisePostDur = 0.0; rmeddis@38: truncateDur = 0.0; %Dr. RF used 0.550 rmeddis@38: rmeddis@38: currentSpeechLevel rmeddis@38: currentNoiseLevel rmeddis@38: rmeddis@38: rmeddis@38: useSpectrogram = false; rmeddis@38: numCoeff = 9; rmeddis@38: rmeddis@38: %************************************************************ rmeddis@38: % plotting handles rmeddis@38: %************************************************************ rmeddis@38: probHaxes = []; rmeddis@38: probHaxesSM = []; rmeddis@38: sacfHaxes = []; rmeddis@38: sacfHaxesSM = []; rmeddis@38: featHaxes = []; rmeddis@38: reconHaxes = []; rmeddis@38: rmeddis@38: %************************************************************ rmeddis@38: % SACF params rmeddis@38: %************************************************************ rmeddis@38: useSACF = false; rmeddis@38: SACFacfTau = 2; % > 1 = Wiegrebe mode rmeddis@38: SACFnBins = 128; rmeddis@38: SACFminLag = 1 / 4000; rmeddis@38: SACFmaxLag = 1 / 50; rmeddis@38: SACFlambda = 10e-3; rmeddis@38: rmeddis@38: %************************************************************ rmeddis@38: % MAP params rmeddis@38: %************************************************************ rmeddis@38: MAProot = fullfile('..'); rmeddis@38: MAPplotGraphs = 0; rmeddis@38: MAPDRNLSave = 0; rmeddis@38: rmeddis@38: MAPopLSR = 0; rmeddis@38: MAPopMSR = 0; rmeddis@38: MAPopHSR = 1; rmeddis@38: rmeddis@38: MAPparamChanges = {}; rmeddis@38: rmeddis@38: %************************************************************ rmeddis@38: % HTK stuff - writing to HTK recogniser format rmeddis@38: %************************************************************ rmeddis@38: frameshift = 10; % shift between frames (ms) rmeddis@38: sampPeriodFromMsFactor = 10000; % appropriate for 10ms frame rate rmeddis@38: paramKind = 9; % HTK USER format for user-defined features (p73 of HTK book) rmeddis@38: rmeddis@38: removeEnergyStatic = false; rmeddis@38: doCMN = false; % Cepstral mean normalisation rmeddis@38: rmeddis@38: %************************************************************ rmeddis@38: % Portable EssexAid params rmeddis@38: %************************************************************ rmeddis@38: bwOct = 1/1; rmeddis@38: filterOrder = 2; rmeddis@38: rmeddis@38: mainGain = [ 1; 1; 1; 1; 1]; % gain in linear units rmeddis@38: TCdBO = [40; 40; 40; 40; 40]; %Compression thresholds (in dB OUTPUT from 2nd filt) rmeddis@38: TMdBO = [10; 10; 10; 10; 10]; %MOC thresholds (in dB OUTPUT from 2nd filt) rmeddis@38: DRNLc = [ 0.2; 0.2; 0.2; 0.2; 0.2;] rmeddis@38: rmeddis@38: ARtau = 0.03; % decay time constant rmeddis@38: ARthresholddB = 85; % dB SPL (input signal level) =>200 to disable rmeddis@38: MOCtau = 0.3; rmeddis@38: MOCfactor = 0.5; %dB per dB OUTPUT rmeddis@38: rmeddis@38: numSamples = 1024; %MAX=6912 rmeddis@38: useAid = 0; rmeddis@38: end rmeddis@38: rmeddis@38: %% ********************************************************* rmeddis@38: % Protected properties - inter/ra class communication rmeddis@38: %************************************************************ rmeddis@38: properties(Access = protected) rmeddis@38: jobLockFid % File identifier for mutex rmeddis@38: rmeddis@38: %Nick C comment on this: rmeddis@38: %OK. The big-endian thing used to work because in the config file rmeddis@38: %'config_tr_zcpa12' there was a flag called NATURALREADORDER that was set to rmeddis@38: %FALSE and thus appeared to override x86 standard: rmeddis@38: %little-endian. Endianess has **!nothing!** to do with win vs *nix rmeddis@38: byteOrder = 'le'; % byte order is big endian rmeddis@38: end rmeddis@38: properties(Dependent = true) rmeddis@38: jobLockTxtFile rmeddis@38: end rmeddis@38: rmeddis@38: %% ********************************************************* rmeddis@38: % methods _ _ _ rmeddis@38: % | | | | | | rmeddis@38: % _ __ ___ ___| |_| |__ ___ __| |___ rmeddis@38: %| '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __| rmeddis@38: %| | | | | | __/ |_| | | | (_) | (_| \__ \ rmeddis@38: %|_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/ rmeddis@38: %************************************************************ rmeddis@38: rmeddis@38: methods rmeddis@38: %% ********************************************************** rmeddis@38: % Constructor rmeddis@38: %************************************************************ rmeddis@38: function obj = cJob(LearnOrRecWavs, jobFolder) rmeddis@38: if nargin < 1 rmeddis@38: warning('job:ambiguouscorpus',... rmeddis@38: 'We need to know whether the (L)earning or (R)ecognition corpus is being used, assuming ''L'''); rmeddis@38: LearnOrRecWavs = 'L'; rmeddis@38: end rmeddis@38: if nargin > 1 rmeddis@38: obj.opFolder = jobFolder; rmeddis@38: else rmeddis@38: if isunix rmeddis@38: if ismac rmeddis@38: obj.opFolder = '~/ASR/exps/_foo'; rmeddis@38: else rmeddis@38: obj.opFolder = '/scratch/nrclark/exps/_foo'; rmeddis@38: end rmeddis@38: else rmeddis@38: obj.opFolder = 'D:\exps\_foo'; rmeddis@38: end rmeddis@38: end rmeddis@38: rmeddis@38: obj = obj.assignWavPaths(LearnOrRecWavs); rmeddis@38: obj = obj.initMAP; rmeddis@38: rmeddis@38: rmeddis@38: end % ------ OF CONSTRUCTOR rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % Set Wav Paths rmeddis@38: %************************************************************ rmeddis@38: function obj = assignWavPaths(obj, LearnOrRecWavs) rmeddis@38: if isunix rmeddis@38: if ismac rmeddis@38: lWAVpath = fullfile('demo_wavs', 'TrainingData-Clean'); rmeddis@38: rWAVpath = fullfile('demo_wavs', 'TripletTestData'); rmeddis@38: obj.noiseFolder = fullfile('demo_wavs', 'noises'); rmeddis@38: else rmeddis@38: lWAVpath = fullfile('demo_wavs', 'TrainingData-Clean'); rmeddis@38: rWAVpath = fullfile('demo_wavs', 'TripletTestData'); rmeddis@38: obj.noiseFolder = fullfile('demo_wavs', 'noises'); rmeddis@38: end rmeddis@38: else rmeddis@38: lWAVpath = 'D:\AURORA digits (wav)\TrainingData-Clean'; rmeddis@38: rWAVpath = 'D:\AURORA digits (wav)\TripletTestData'; rmeddis@38: obj.noiseFolder = 'D:\AURORA digits (wav)\noises'; rmeddis@38: end rmeddis@38: rmeddis@38: if strcmpi(LearnOrRecWavs, 'l') rmeddis@38: obj.wavFolder = lWAVpath; rmeddis@38: elseif strcmpi(LearnOrRecWavs, 'r') rmeddis@38: obj.wavFolder = rWAVpath; rmeddis@38: else rmeddis@38: error('First argument to constructor must be ''L'' or ''R'''); rmeddis@38: end rmeddis@38: end rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % mutex related _ rmeddis@38: % | | rmeddis@38: % _ __ ___ _ _| |_ _____ __ rmeddis@38: % | '_ ` _ \| | | | __/ _ \ \/ / rmeddis@38: % | | | | | | |_| | || __/> < rmeddis@38: % |_| |_| |_|\__,_|\__\___/_/\_\ rmeddis@38: %************************************************************ rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % lockJobList - File mutex protecting from race conditions rmeddis@38: %************************************************************ rmeddis@38: function obj = lockJobList(obj) rmeddis@38: lockON = false; rmeddis@38: while (~lockON) rmeddis@38: if numel(dir(obj.jobLockTxtFile)) %Check to see if lock already in place rmeddis@38: wTime = randi(750)+250; %3,000-10,000 ms rmeddis@38: disp(['File mutex in place. Retrying in ' num2str(wTime) ' ms']) rmeddis@38: pause(wTime/1000); rmeddis@38: else rmeddis@38: obj.jobLockFid = fopen(obj.jobLockTxtFile,'w'); rmeddis@38: disp('Job locked'); rmeddis@38: pause(50/1000); rmeddis@38: lockON = true; rmeddis@38: end rmeddis@38: end rmeddis@38: fclose(obj.jobLockFid); rmeddis@38: end% ------ OF LOCKJOBLIST rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % unlockJobList - unlocks for other workers rmeddis@38: %************************************************************ rmeddis@38: function obj = unlockJobList(obj) rmeddis@38: lockON = logical(numel(dir(obj.jobLockTxtFile))); rmeddis@38: while(lockON) rmeddis@38: try rmeddis@38: delete(obj.jobLockTxtFile); rmeddis@38: disp('Job unlocked'); rmeddis@38: pause(50/1000); rmeddis@38: lockON = false; rmeddis@38: catch %#ok rmeddis@38: disp('Unjamming lock - retrying immediately') rmeddis@38: pause(200/1000) rmeddis@38: end rmeddis@38: end rmeddis@38: end% ------ OF UNLOCKJOBLIST rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % storeSelf - save a copy of this object in opFolder rmeddis@38: %************************************************************ rmeddis@38: function storeSelf(obj) rmeddis@38: doneSaving = false; rmeddis@38: while(~doneSaving) rmeddis@38: try rmeddis@38: save(fullfile(obj.opFolder, 'jobObject'), 'obj'); rmeddis@38: doneSaving = true; rmeddis@38: catch %#ok rmeddis@38: wTime = randi(3800)+200; %200-4000 ms rmeddis@38: disp(['Save collision (THIS IS VERY VERY BAD - have you not used the mutex?). Retrying in ' num2str(wTime) ' ms']); rmeddis@38: pause(wTime/1000); rmeddis@38: end rmeddis@38: end rmeddis@38: end % ------ OF STORESELF rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % loadSelf - load a copy of this object from opFolder rmeddis@38: %************************************************************ rmeddis@38: function obj = loadSelf(obj) rmeddis@38: doneLoading = false; rmeddis@38: while(~doneLoading) rmeddis@38: try rmeddis@38: load(fullfile(obj.opFolder,'jobObject')); rmeddis@38: doneLoading = true; rmeddis@38: catch %#ok rmeddis@38: wTime = randi(3800)+200; %200-4000 ms rmeddis@38: disp(['Load collision (THIS IS VERY VERY BAD - have you not used the mutex?). Retrying in ' num2str(wTime) ' ms']) rmeddis@38: pause(wTime/1000); rmeddis@38: end rmeddis@38: end rmeddis@38: end% ------ OF LOADSELF rmeddis@38: rmeddis@38: rmeddis@38: %% ********************************************************* rmeddis@38: % _ _ _ rmeddis@38: % | | | | (_) rmeddis@38: % | |__ ___ _ _ ___ ___| | _____ ___ _ __ _ _ __ __ _ rmeddis@38: % | '_ \ / _ \| | | / __|/ _ \ |/ / _ \/ _ \ '_ \| | '_ \ / _` | rmeddis@38: % | | | | (_) | |_| \__ \ __/ < __/ __/ |_) | | | | | (_| | rmeddis@38: % |_| |_|\___/ \__,_|___/\___|_|\_\___|\___| .__/|_|_| |_|\__, | rmeddis@38: % | | __/ | rmeddis@38: % |_| |___/ rmeddis@38: %************************************************************ rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % get method for dependtent var jobLockTxtFile rmeddis@38: %************************************************************ rmeddis@38: function value = get.jobLockTxtFile(obj) rmeddis@38: %Name the MUTEX file here rmeddis@38: value = [fullfile(obj.opFolder, 'jobLock') '.txt']; rmeddis@38: end rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % checkStatus - see how much of the current job is complete rmeddis@38: %************************************************************ rmeddis@38: function checkStatus(obj) rmeddis@38: NOWopen = sum(obj.todoStatus==0); %Starting from R2010b, Matlab supports enumerations. For now, we resort to integers for compatibility. rmeddis@38: NOWpend = sum(obj.todoStatus==1); rmeddis@38: NOWdone = sum(obj.todoStatus==2); rmeddis@38: rmeddis@38: zz = clock; rmeddis@38: disp([num2str(zz(3)) '-' num2str(zz(2)) '-' num2str(zz(1)) ' ' num2str(zz(4)) ':' num2str(zz(5))]) rmeddis@38: disp(['CURRENT JOB:' obj.opFolder]); rmeddis@38: disp(' ') rmeddis@38: disp(['open - ' num2str(NOWopen) ' || pending - ' num2str(NOWpend) ' || complete - ' num2str(NOWdone)]) rmeddis@38: rmeddis@38: % ---- PROGRESS BAR ---- rmeddis@38: pcDone = 100*NOWdone/obj.numWavs; rmeddis@38: progBarLength = 40; rmeddis@38: charBars = repmat('=',1,floor(pcDone/100 * progBarLength)); rmeddis@38: charWhiteSpace = repmat(' ',1, progBarLength - numel(charBars)); rmeddis@38: disp(' ') rmeddis@38: disp([' -[' charBars charWhiteSpace ']- ' num2str(pcDone, '%0.1f') '%']) rmeddis@38: disp(' ') rmeddis@38: disp(' ') rmeddis@38: % -- END PROGRESS BAR --- rmeddis@38: end% ------ OF CHECKSTATUS rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % initMAP - add MAP stuff to path rmeddis@38: %************************************************************ rmeddis@38: function obj = initMAP(obj) rmeddis@38: addpath(...fullfile(obj.MAProot, 'modules'),... rmeddis@38: fullfile(obj.MAProot, 'utilities'),... rmeddis@38: fullfile(obj.MAProot, 'MAP'),... rmeddis@38: fullfile(obj.MAProot, 'parameterStore')); rmeddis@38: end % ------ OF INIT MAP rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % assign files to testing and training sets rmeddis@38: %************************************************************ rmeddis@38: function obj = assignFiles(obj) rmeddis@38: speechWavs = dir(fullfile(obj.wavFolder, '*.wav')); rmeddis@38: assert(obj.numWavs <= size(speechWavs, 1) ,... rmeddis@38: 'not enough files available in the corpus'); % make sure we have enough wavs rmeddis@38: rmeddis@38: randomWavs = rand(1, size(speechWavs, 1)); rmeddis@38: [~,b] = sort(randomWavs); rmeddis@38: trainFileIdx = b(1:obj.numWavs); rmeddis@38: rmeddis@38: obj.wavList = speechWavs(trainFileIdx); %This is a record of all of the wavs that should be done rmeddis@38: rmeddis@38: %Starting from R2010b, Matlab should support enumerated types. For now we rmeddis@38: %use integers for compatibility. rmeddis@38: %0=open, 1=processing, 2=done rmeddis@38: obj.todoStatus = zeros(numel(obj.wavList), 1); rmeddis@38: rmeddis@38: end % ------ OF ASSIGN FILES rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % generate feature rmeddis@38: %************************************************************ rmeddis@38: function obj = genFeat(obj, currentWav) rmeddis@38: fprintf(1,'Processing: %s \n', currentWav); rmeddis@38: if strcmpi(obj.speechDist,'Gaussian') rmeddis@38: tempSpeechLev = obj.speechLevToUse + obj.speechLevStd*randn; rmeddis@38: elseif strcmpi(obj.speechDist,'Uniform') rmeddis@38: % for a uniform distribution, the standard deviation is rmeddis@38: % range/sqrt(12) rmeddis@38: % http://mathforum.org/library/drmath/view/52066.html rmeddis@38: tempSpeechLev = obj.speechLevToUse - obj.speechLevStd*sqrt(12)/2 + obj.speechLevStd*sqrt(12)*rand; rmeddis@38: elseif strcmpi(obj.speechDist,'None') rmeddis@38: tempSpeechLev = obj.speechLevToUse; rmeddis@38: end rmeddis@38: rmeddis@38: if strcmpi(obj.noiseDist,'Gaussian') rmeddis@38: tempNoiseLev = speechLev - obj.meanSNR + obj.noiseLevStd*randn; rmeddis@38: elseif strcmpi(obj.noiseDist,'Uniform') rmeddis@38: tempNoiseLev = tempSpeechLev - obj.meanSNR - obj.noiseLevStd*sqrt(12)/2 + obj.noiseLevStd*sqrt(12)*rand; rmeddis@38: elseif strcmpi(obj.noiseDist,'None') rmeddis@38: tempNoiseLev = obj.noiseLevToUse; rmeddis@38: end rmeddis@38: rmeddis@38: disp(['Current speech level = ' num2str(tempSpeechLev)]); rmeddis@38: disp(['Current noise level = ' num2str(tempNoiseLev)]); rmeddis@38: rmeddis@38: obj.currentSpeechLevel = tempSpeechLev; rmeddis@38: obj.currentNoiseLevel = tempNoiseLev; rmeddis@38: [finalFeatures, ~] = processWavs(obj, currentWav); %discard the output from ANprobabilityResponse and method using ~ rmeddis@38: opForHTK(obj, currentWav, finalFeatures); rmeddis@38: end % ------ OF GENFEAT rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % write o/p in HTK friendly format rmeddis@38: %************************************************************ rmeddis@38: function obj = opForHTK(obj, currentWav, featureData) rmeddis@38: rmeddis@38: featureName = strrep(currentWav, '.wav','.map'); rmeddis@38: targetFilename = fullfile(obj.opFolder, featureName); rmeddis@38: rmeddis@38: % write in a format HTK compliant for the recogniser to use rmeddis@38: obj.writeHTK(... rmeddis@38: targetFilename,... rmeddis@38: featureData,... rmeddis@38: size(featureData,2),... rmeddis@38: obj.frameshift*obj.sampPeriodFromMsFactor,... rmeddis@38: size(featureData,1)*4,... rmeddis@38: obj.paramKind,... rmeddis@38: obj.byteOrder); rmeddis@38: end % ------ OF opForHTK rmeddis@38: rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % _ _ _ rmeddis@38: % (_) | | (_) rmeddis@38: % ___ _ __ _ _ __ __ _| | _ __ _ __ ___ ___ ___ ___ ___ _ _ __ __ _ rmeddis@38: % / __| |/ _` | '_ \ / _` | | | '_ \| '__/ _ \ / __/ _ \/ __/ __| | '_ \ / _` | rmeddis@38: % \__ \ | (_| | | | | (_| | | | |_) | | | (_) | (_| __/\__ \__ \ | | | | (_| | rmeddis@38: % |___/_|\__, |_| |_|\__,_|_| | .__/|_| \___/ \___\___||___/___/_|_| |_|\__, | rmeddis@38: % __/ | | | __/ | rmeddis@38: % |___/ |_| |___/ rmeddis@38: %************************************************************ rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % getStimulus - what it says on the tin rmeddis@38: %************************************************************ rmeddis@38: function [stimulus, sampleRate] = getStimulus(obj, currentWav) rmeddis@38: rmeddis@38: % getStimulus.m - NC Aug 2010 rmeddis@38: % Modified version of Rob's script to include: rmeddis@38: % rmeddis@38: % 1)Signal and noise samples with different sample rate. The component with rmeddis@38: % lower sample rate is upsampled to the rate of that with the rmeddis@38: % higher rate. rmeddis@38: % 2) Clearer level setting rmeddis@38: % 3) Parameter to change noise intro duration rmeddis@38: % 4) Noise padding at end of stimulus rmeddis@38: % rmeddis@38: % ORIGINAL HEADER: rmeddis@38: % getStimulus.m rmeddis@38: % rmeddis@38: % Robert T. Ferry rmeddis@38: % 13th May 2009 rmeddis@38: rmeddis@38: % Set levels rmeddis@38: [speech speechSampleRate] = wavread(fullfile(obj.wavFolder, currentWav )); rmeddis@38: speech = speech./sqrt(mean(speech.^2)); %Normalize RMS to 1 rmeddis@38: speech = speech * 20e-6 * 10^(obj.currentSpeechLevel/20); %Convert RMS to pascals at desired level rmeddis@38: %disp(20*log10(sqrt(mean(speech.^2))/20e-6)) rmeddis@38: rmeddis@38: [noise noiseSampleRate] = wavread(fullfile(obj.noiseFolder, obj.noiseName )); rmeddis@38: noise = noise./sqrt(mean(noise.^2)); %Normalize RMS to 1 rmeddis@38: noise = noise * 20e-6*10^(obj.currentNoiseLevel/20); %Convert RMS to pascals at desired level rmeddis@38: %disp(20*log10(sqrt(mean(noise.^2))/20e-6)) rmeddis@38: rmeddis@38: % Do sample rate conversion if needed rmeddis@38: % Will always convert stimulus component with lower rate up to that with rmeddis@38: % higher rate. rmeddis@38: if speechSampleRate > noiseSampleRate rmeddis@38: % disp('S>N') rmeddis@38: [p,q] = rat(speechSampleRate/noiseSampleRate,0.0001); rmeddis@38: noise = resample(noise, p, q); rmeddis@38: noiseSampleRate = speechSampleRate; rmeddis@38: elseif noiseSampleRate > speechSampleRate rmeddis@38: % disp('N>S') rmeddis@38: [p,q] = rat(noiseSampleRate/speechSampleRate,0.0001); rmeddis@38: speech = resample(speech, p, q); rmeddis@38: speechSampleRate = noiseSampleRate; %#ok rmeddis@38: else rmeddis@38: %DO NOTHING BUT ASSERT rmeddis@38: assert(speechSampleRate == noiseSampleRate); rmeddis@38: end rmeddis@38: sampleRate = noiseSampleRate; rmeddis@38: dt = 1/sampleRate; rmeddis@38: rmeddis@38: % mix stimuli rmeddis@38: % Everything from here down (mostly) is RTF's original rmeddis@38: silenceStart = floor(obj.noisePreDur*sampleRate); rmeddis@38: silenceEnd = floor(obj.noisePostDur*sampleRate); rmeddis@38: rmeddis@38: silencePointsStart = zeros(silenceStart,1); rmeddis@38: silencePointsEnd = zeros(silenceEnd,1); rmeddis@38: rmeddis@38: speech = [silencePointsStart; speech; silencePointsEnd]; rmeddis@38: rmeddis@38: stimLength = length(speech); rmeddis@38: noiseLength = length(noise); rmeddis@38: if obj.freezeNoise rmeddis@38: idx = 1; rmeddis@38: else rmeddis@38: idx = ceil(rand*(noiseLength-stimLength)); rmeddis@38: end rmeddis@38: noise = noise(idx:idx+stimLength-1); rmeddis@38: rmeddis@38: stimulus = speech+noise; rmeddis@38: rmeddis@38: % add ramps to noise rmeddis@38: stimInNoiseTime = dt:dt:dt*length(stimulus); rmeddis@38: rampDuration = 0.100; rmeddis@38: rampTime = dt : dt : rampDuration; rmeddis@38: ramp = [0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(stimInNoiseTime)-length(rampTime))]; rmeddis@38: stimulus = stimulus'.*ramp; rmeddis@38: ramp = fliplr(ramp); rmeddis@38: stimulus = stimulus.*ramp; rmeddis@38: stimulus = stimulus'; rmeddis@38: %disp(20*log10(sqrt(mean(stimulus.^2))/20e-6)) rmeddis@38: rmeddis@38: % add additional silent period to start of stimulus for model to 'settle down' rmeddis@38: additionalSilenceLength = round(0.050*sampleRate); rmeddis@38: additionalSilencePoints = zeros(additionalSilenceLength,1); rmeddis@38: stimulus = [additionalSilencePoints; stimulus]'; %now rotated. rmeddis@38: %disp(20*log10(sqrt(mean(stimulus.^2))/20e-6)) rmeddis@38: end% ------ OF GETSTIMULUS rmeddis@38: rmeddis@38: rmeddis@38: %% ********************************************************** rmeddis@38: % processWavs - do all of the MAP signal processing rmeddis@38: %************************************************************ rmeddis@38: function [finalFeatures, ANprobabilityResponse] = processWavs(obj, currentWav) rmeddis@38: rmeddis@38: %********************************************************** rmeddis@38: % FIRST GET THE STIMULUS rmeddis@38: %********************************************************** rmeddis@38: [stimulus, sampleRate] = obj.getStimulus(currentWav); rmeddis@38: rmeddis@38: %********************************************************** rmeddis@38: % NOW TO LOAD IN THE HEARING AID rmeddis@38: %********************************************************** rmeddis@38: if obj.useAid rmeddis@38: stimulus = [stimulus; stimulus]'; %EsxAid requires stereo stim rmeddis@38: stimulus = EssexAid_guiEmulatorWrapper(stimulus, sampleRate, obj); rmeddis@38: stimulus = stimulus(1,:); %convert back to mono rmeddis@38: end rmeddis@38: rmeddis@38: AN_spikesOrProbability = 'probability'; rmeddis@38: rmeddis@38: if obj.useSpectrogram rmeddis@38: lowestBF=100; highestBF= 4500; numChannels=30; rmeddis@38: F=round(logspace(log10(lowestBF),log10(highestBF),numChannels)); rmeddis@38: rmeddis@38: nfft = 1024; rmeddis@38: hopSamples = 64; rmeddis@38: noverlap = nfft - hopSamples; rmeddis@38: dt = hopSamples/sampleRate; rmeddis@38: [~,~,~,P] = spectrogram(stimulus,nfft,noverlap,F,sampleRate); rmeddis@38: rmeddis@38: ANprobabilityResponse = 10*log10( abs(P) / ((20e-6)^2) ); %now correct [(a^2)/(b^2) = (a/b)^2] rmeddis@38: rmeddis@38: else rmeddis@38: [ANprobabilityResponse, dt, myBFlist] = MAPwrap(stimulus, sampleRate, -1, obj.participant, AN_spikesOrProbability, obj.MAPparamChanges); rmeddis@38: end rmeddis@38: nChannels = numel(myBFlist); rmeddis@38: rmeddis@38: rmeddis@38: time_ANresponse = dt:dt:dt*size(ANprobabilityResponse,2); rmeddis@38: idx = time_ANresponse > obj.truncateDur; %RTF had this @ 0.550 rmeddis@38: ANprobabilityResponse = ANprobabilityResponse(:, idx); rmeddis@38: rmeddis@38: rmeddis@38: if ~obj.useSpectrogram rmeddis@38: if obj.MAPopLSR && ~obj.MAPopHSR rmeddis@38: ANprobabilityResponse = ANprobabilityResponse(1:nChannels, :); %use LSR rmeddis@38: end rmeddis@38: if ~obj.MAPopLSR && obj.MAPopHSR rmeddis@38: ANprobabilityResponse = ANprobabilityResponse(end-nChannels+1:end, :); %or use HSR rmeddis@38: end rmeddis@38: if obj.MAPopMSR rmeddis@38: assert(0,'Not working with MSR at the mo') rmeddis@38: end rmeddis@38: end rmeddis@38: rmeddis@38: % OPTIONAL PLOTTING rmeddis@38: YTickIdx = 1:floor(numel(myBFlist)/6):numel(myBFlist); rmeddis@38: YTickIdxRev = numel(myBFlist)+1-YTickIdx; rmeddis@38: if ~isempty(obj.probHaxes) rmeddis@38: axes(obj.probHaxes); %#ok rmeddis@38: imagesc(ANprobabilityResponse); colorbar('peer', obj.probHaxes) rmeddis@38: set(obj.probHaxes, 'YTick', YTickIdx); rmeddis@38: set(obj.probHaxes, 'YTickLabel', num2str( myBFlist(YTickIdxRev)' )); rmeddis@38: ylabel('cf in Hz') rmeddis@38: end rmeddis@38: rmeddis@38: % OPTIONAL PLOTTING SMOOTHED rmeddis@38: if ~isempty(obj.probHaxesSM) rmeddis@38: axes(obj.probHaxesSM); %#ok rmeddis@38: anSM=flipud(obj.makeANsmooth(ANprobabilityResponse, 1/dt)); rmeddis@38: imagesc((1:size(anSM,2))./100,1:size(ANprobabilityResponse,1),anSM); rmeddis@38: set(obj.probHaxesSM, 'YTick', YTickIdx); rmeddis@38: set(obj.probHaxesSM, 'YTickLabel', num2str( myBFlist(YTickIdxRev)' )); rmeddis@38: shading(obj.probHaxesSM, 'flat'); colorbar('peer', obj.probHaxesSM) rmeddis@38: ylabel('cf (Hz)') rmeddis@38: xlabel('Time (s)') rmeddis@38: end rmeddis@38: rmeddis@38: rmeddis@38: %********************************************************** rmeddis@38: % optional SACF stage rmeddis@38: %********************************************************** rmeddis@38: if obj.useSACF rmeddis@38: rmeddis@38: % A slightly ugly copying is needed rmeddis@38: SACFparams.lambda = obj.SACFlambda; rmeddis@38: SACFparams.acfTau = obj.SACFacfTau; rmeddis@38: SACFparams.lags = logspace(log10(obj.SACFminLag),log10(obj.SACFmaxLag),obj.SACFnBins); rmeddis@38: SACFparams.lags = linspace(obj.SACFminLag, obj.SACFmaxLag,obj.SACFnBins ); rmeddis@38: rmeddis@38: SACFmethod.dt = dt; rmeddis@38: SACFmethod.nonlinCF = myBFlist; rmeddis@38: rmeddis@38: %This is slightly misleading as the ANprob is now a SACF rmeddis@38: [ANprobabilityResponse, ~, ~, ~] = filteredSACF(ANprobabilityResponse, SACFmethod, SACFparams); rmeddis@38: rmeddis@38: % OPTIONAL PLOTTING rmeddis@38: YTickIdx = 1:floor(obj.SACFnBins/6):obj.SACFnBins; rmeddis@38: %YTickIdxRev = obj.SACFnBins+1-YTickIdx; rmeddis@38: if ~isempty(obj.sacfHaxes) rmeddis@38: axes(obj.sacfHaxes); %#ok rmeddis@38: imagesc(flipud(ANprobabilityResponse)); shading(obj.sacfHaxes, 'flat'); colorbar('peer', obj.sacfHaxes) rmeddis@38: set(obj.sacfHaxes, 'YTick', YTickIdx); rmeddis@38: set(obj.sacfHaxes, 'YTickLabel', num2str( 1./SACFparams.lags(YTickIdx)' ,'%0.1f' )); rmeddis@38: ylabel('Pitch in Hz') rmeddis@38: end rmeddis@38: rmeddis@38: % OPTIONAL PLOTTING SMOOTHED rmeddis@38: if ~isempty(obj.sacfHaxesSM) rmeddis@38: axes(obj.sacfHaxesSM); %#ok rmeddis@38: imagesc(flipud(obj.makeANsmooth(ANprobabilityResponse, 1/dt))); shading(obj.sacfHaxesSM, 'flat'); colorbar('peer', obj.sacfHaxesSM) rmeddis@38: set(obj.sacfHaxesSM, 'YTick', YTickIdx); rmeddis@38: set(obj.sacfHaxesSM, 'YTickLabel', num2str( 1./SACFparams.lags(YTickIdx)' ,'%0.1f' )); rmeddis@38: ylabel('Pitch in Hz') rmeddis@38: end rmeddis@38: rmeddis@38: end rmeddis@38: rmeddis@38: rmeddis@38: finalFeatures = obj.makeANfeatures( ... rmeddis@38: obj.makeANsmooth(ANprobabilityResponse, 1/dt), obj.numCoeff ); rmeddis@38: rmeddis@38: if obj.removeEnergyStatic rmeddis@38: finalFeatures = finalFeatures(2:end,:); rmeddis@38: % disp(size(finalFeatures)) rmeddis@38: end rmeddis@38: rmeddis@38: if obj.doCMN rmeddis@38: m = mean(finalFeatures,2); rmeddis@38: finalFeatures = finalFeatures - repmat(m,1,size(finalFeatures,2)); rmeddis@38: end rmeddis@38: rmeddis@38: % OPTIONAL PLOTTING rmeddis@38: if ~isempty(obj.featHaxes) rmeddis@38: pcolor(obj.featHaxes, finalFeatures); shading(obj.featHaxes, 'flat'); colorbar('peer', obj.featHaxes) rmeddis@38: end rmeddis@38: if ~isempty(obj.reconHaxes) rmeddis@38: reconsData = idct(finalFeatures,obj.SACFnBins); rmeddis@38: axes(obj.reconHaxes); %#ok rmeddis@38: imagesc(flipud( reconsData )); rmeddis@38: end rmeddis@38: rmeddis@38: opForHTK(obj, currentWav, finalFeatures); rmeddis@38: end % ------ OF PROCESSWAVS rmeddis@38: rmeddis@38: end % ------ OF METHODS rmeddis@38: rmeddis@38: %% ********************************************************* rmeddis@38: % _ _ _ _ _ _ rmeddis@38: % | | | | (_) | | | | | | rmeddis@38: % ___| |_ __ _| |_ _ ___ _ __ ___ ___| |_| |__ ___ __| |___ rmeddis@38: % / __| __/ _` | __| |/ __| | '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __| rmeddis@38: % \__ \ || (_| | |_| | (__ | | | | | | __/ |_| | | | (_) | (_| \__ \ rmeddis@38: % |___/\__\__,_|\__|_|\___| |_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/ rmeddis@38: %************************************************************ rmeddis@38: rmeddis@38: methods(Static) rmeddis@38: %% ******************************************************** rmeddis@38: % makeANsmooth - smooth the AN response into hanning windowed chunks rmeddis@38: %********************************************************** rmeddis@38: function ANsmooth = makeANsmooth(ANresponse, sampleRate, winSize, hopSize) rmeddis@38: if nargin < 3 rmeddis@38: winSize = 25; %default 25 ms window rmeddis@38: end rmeddis@38: if nargin < 4 rmeddis@38: hopSize = 10; %default 10 ms jump between windows rmeddis@38: end rmeddis@38: rmeddis@38: winSizeSamples = round(winSize*sampleRate/1000); rmeddis@38: hopSizeSamples = round(hopSize*sampleRate/1000); rmeddis@38: rmeddis@38: % smooth rmeddis@38: hann = cJob.NRC_hanning(winSizeSamples); rmeddis@38: rmeddis@38: ANsmooth = [];%Cannot pre-allocate a size as it is unknown until the enframing rmeddis@38: for chan = 1:size(ANresponse,1) rmeddis@38: f = cJob.enframe(ANresponse(chan,:), hann, hopSizeSamples); rmeddis@38: ANsmooth(chan,:) = mean(f,2)'; %#ok see above comment rmeddis@38: end rmeddis@38: end% ------ OF makeANsmooth rmeddis@38: rmeddis@38: %% ******************************************************** rmeddis@38: % makeANfeatures - dct wizardry rmeddis@38: %********************************************************** rmeddis@38: function ANfeatures = makeANfeatures(ANrate, numCoeff) rmeddis@38: % make feature vectors rmeddis@38: features = cJob.GJB_dct(ANrate); rmeddis@38: ANfeatures = features(1:numCoeff,:); rmeddis@38: end % ------ OF makeANfeatures rmeddis@38: rmeddis@38: %% ************************************************************************ rmeddis@38: % enframe - AVOID SIGNAL PROCESSING TOOLBOX buffer function rmeddis@38: %************************************************************************** rmeddis@38: function f=enframe(x,win,inc) rmeddis@38: %ENFRAME split signal up into (overlapping) frames: one per row. F=(X,WIN,INC) rmeddis@38: % rmeddis@38: % F = ENFRAME(X,LEN) splits the vector X(:) up into rmeddis@38: % frames. Each frame is of length LEN and occupies rmeddis@38: % one row of the output matrix. The last few frames of X rmeddis@38: % will be ignored if its length is not divisible by LEN. rmeddis@38: % It is an error if X is shorter than LEN. rmeddis@38: % rmeddis@38: % F = ENFRAME(X,LEN,INC) has frames beginning at increments of INC rmeddis@38: % The centre of frame I is X((I-1)*INC+(LEN+1)/2) for I=1,2,... rmeddis@38: % The number of frames is fix((length(X)-LEN+INC)/INC) rmeddis@38: % rmeddis@38: % F = ENFRAME(X,WINDOW) or ENFRAME(X,WINDOW,INC) multiplies rmeddis@38: % each frame by WINDOW(:) rmeddis@38: rmeddis@38: % Copyright (C) Mike Brookes 1997 rmeddis@38: % Version: $Id: enframe.m,v 1.4 2006/06/22 19:07:50 dmb Exp $ rmeddis@38: % rmeddis@38: % VOICEBOX is a MATLAB toolbox for speech processing. rmeddis@38: % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html rmeddis@38: % rmeddis@38: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rmeddis@38: % This program is free software; you can redistribute it and/or modify rmeddis@38: % it under the terms of the GNU General Public License as published by rmeddis@38: % the Free Software Foundation; either version 2 of the License, or rmeddis@38: % (at your option) any later version. rmeddis@38: % rmeddis@38: % This program is distributed in the hope that it will be useful, rmeddis@38: % but WITHOUT ANY WARRANTY; without even the implied warranty of rmeddis@38: % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the rmeddis@38: % GNU General Public License for more details. rmeddis@38: % rmeddis@38: % You can obtain a copy of the GNU General Public License from rmeddis@38: % http://www.gnu.org/copyleft/gpl.html or by writing to rmeddis@38: % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. rmeddis@38: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rmeddis@38: rmeddis@38: nx=length(x(:)); rmeddis@38: nwin=length(win); rmeddis@38: if (nwin == 1) rmeddis@38: len = win; rmeddis@38: else rmeddis@38: len = nwin; rmeddis@38: end rmeddis@38: if (nargin < 3) rmeddis@38: inc = len; rmeddis@38: end rmeddis@38: nf = fix((nx-len+inc)/inc); rmeddis@38: f=zeros(nf,len); rmeddis@38: indf= inc*(0:(nf-1)).'; rmeddis@38: inds = (1:len); rmeddis@38: f(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:)); rmeddis@38: if (nwin > 1) rmeddis@38: w = win(:)'; rmeddis@38: f = f .* w(ones(nf,1),:); rmeddis@38: end rmeddis@38: end % ------ OF ENFRAME rmeddis@38: rmeddis@38: %% ************************************************************************ rmeddis@38: % GJB_dct - AVOID SIGNAL PROCESSING TOOLBOX rmeddis@38: %************************************************************************** rmeddis@38: function b=GJB_dct(a,n) rmeddis@38: rmeddis@38: if nargin == 0, rmeddis@38: error('Not enough input arguments.'); rmeddis@38: end rmeddis@38: rmeddis@38: if isempty(a) rmeddis@38: b = []; rmeddis@38: return rmeddis@38: end rmeddis@38: rmeddis@38: % If input is a vector, make it a column: rmeddis@38: do_trans = (size(a,1) == 1); rmeddis@38: if do_trans, a = a(:); end rmeddis@38: rmeddis@38: if nargin==1, rmeddis@38: n = size(a,1); rmeddis@38: end rmeddis@38: m = size(a,2); rmeddis@38: rmeddis@38: % Pad or truncate input if necessary rmeddis@38: if size(a,1) by anonymous Sue (2001) rmeddis@38: %************************************************************************** rmeddis@38: function retcode = writeHTK(filename, htkdata, nFrames, sampPeriod, SampSize, ParamKind, byte_order) rmeddis@38: % Write an HTK format file. rmeddis@38: % rmeddis@38: % Input parameters: rmeddis@38: % filename HTK data file rmeddis@38: % htkdata HTK data read: an m x n matrix with rmeddis@38: % m = no. of channels rmeddis@38: % n = no. of frames rmeddis@38: % The following are from the HTK header (see HTK manual): rmeddis@38: % nFrames no. of frames (samples) rmeddis@38: % sampPeriod sample period (in 100 ns units?) rmeddis@38: % SampSize sample size rmeddis@38: % ParamKind parameter kind code rmeddis@38: % rmeddis@38: % byteorder 'be' for big-endian (typical for Unix) (default) rmeddis@38: % 'le' for little-endian (typical for MSWindows) rmeddis@38: % rmeddis@38: % Output parameters: rmeddis@38: % retcode 0 if successful rmeddis@38: rmeddis@38: % Written by Sue 17/12/01 rmeddis@38: rmeddis@38: retcode=-1; % initialise in case of error rmeddis@38: if nargin < 6 rmeddis@38: fprintf('Usage: %s(filename, htkdata, nFrames, sampPeriod, SampSize, ParamKind [, byte_order])', mfilename); rmeddis@38: end; rmeddis@38: rmeddis@38: % Default to big-endian (HTK format) rmeddis@38: if nargin < 7 rmeddis@38: byte_order = 'be'; rmeddis@38: end; rmeddis@38: rmeddis@38: fid = fopen (filename, 'w', sprintf('ieee-%s', byte_order)); rmeddis@38: if fid < 1 rmeddis@38: fprintf('%s: can''t open output file %s\n', mfilename, filename); rmeddis@38: return rmeddis@38: end rmeddis@38: rmeddis@38: % Write header rmeddis@38: fwrite (fid, nFrames, 'int32'); %nSamples in HTK rmeddis@38: fwrite (fid, sampPeriod, 'int32'); rmeddis@38: fwrite (fid, SampSize, 'int16'); rmeddis@38: fwrite (fid, ParamKind, 'int16'); rmeddis@38: rmeddis@38: % Write actual data rmeddis@38: fwrite(fid, htkdata, 'float32'); rmeddis@38: rmeddis@38: fclose(fid); rmeddis@38: rmeddis@38: retcode=0; rmeddis@38: end% ------ OF WRITEHTK rmeddis@38: rmeddis@38: %% ************************************************************************ rmeddis@38: % readHTK - just incase you ever want to go backwards rmeddis@38: %************************************************************************** rmeddis@38: function [htkdata,nframes,sampPeriod,sampSize,paramKind] = readHTK(filename,byte_order) rmeddis@38: rmeddis@38: if nargin<2 rmeddis@38: byte_order = 'be'; rmeddis@38: end rmeddis@38: rmeddis@38: fid = fopen(filename,'r',sprintf('ieee-%s',byte_order)); rmeddis@38: rmeddis@38: nframes = fread(fid,1,'int32'); rmeddis@38: sampPeriod = fread(fid,1,'int32'); rmeddis@38: sampSize = fread(fid,1,'int16'); rmeddis@38: paramKind = fread(fid,1,'int16'); rmeddis@38: rmeddis@38: % read the data rmeddis@38: rmeddis@38: htkdata = fread(fid,nframes*(sampSize/4),'float32'); rmeddis@38: htkdata = reshape(htkdata,sampSize/4,nframes); rmeddis@38: fclose(fid); rmeddis@38: end % ------ OF READHTK rmeddis@38: rmeddis@38: end % ------ OF STATIC METHODS rmeddis@38: rmeddis@38: end % ------ OF CLASS