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