Mercurial > hg > emotion-detection-top-level
changeset 0:ea0c737c6323
first commit
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Collation/get_PRAAT.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,171 @@ +function [] = get_PRAAT( dirName, statsFileID ) + +% this function collates the results of all values calculated using the PRAAT software + + % identify the speaker in the stats file + fprintf( statsFileID, '%s ', dirName ); + +% -------------- get the jitter metrics ------------------- +% JITTER: ddp \t local \t ppq5 \t rap \t + FileName = [ dirName '_jitter_ddp.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING JITTER DDP METRICS FILE'); + fprintf( statsFileID, '\t ddp missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + FileName = [ dirName '_jitter_local.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING JITTER LOCAL METRICS FILE'); + fprintf( statsFileID, '\t local missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + FileName = [ dirName '_jitter_ppq5.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING JITTER PPQ5 METRICS FILE'); + fprintf( statsFileID, '\t ppq5 missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + FileName = [ dirName '_jitter_rap.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING JITTER RAP METRICS FILE'); + fprintf( statsFileID, '\t rap missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + + +%-------------- get the shimmer metrics ---------------------- +% SHIMMER: local \t dda \t apq3 \t apq5 \t apq11 + FileName = [ dirName '_shimmer_local.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING SHIMMER METRICS FILE'); + fprintf( statsFileID, '\t local missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + FileName = [ dirName '_shimmer_dda.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING SHIMMER METRICS FILE'); + fprintf( statsFileID, '\t dda missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + FileName = [ dirName '_shimmer_apq3.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING SHIMMER METRICS FILE'); + fprintf( statsFileID, '\t apq3 missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + FileName = [ dirName '_shimmer_apq5.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING SHIMMER METRICS FILE'); + fprintf( statsFileID, '\t apq5 missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + + FileName = [ dirName '_shimmer_apq11.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING SHIMMER METRICS FILE'); + fprintf( statsFileID, '\t apq11 missing'); + else + shimmer = fscanf( FileID, '%f', inf ); + fprintf( statsFileID, '\t %f ', shimmer ); + end + fclose(FileID); + +%-------------- get the formant metrics ---------------------- + + % need to discard all formant information for unvoiced frames. + + FileName = [ dirName '_Formant.txt']; + FileID = fopen( FileName ); + if( FileID <= 0 ) %does the file exist? + % if not + disp('WARNING: MISSING FORMANT METRICS FILE'); + fprintf( statsFileID, '\t formant missing'); + else + % file format is not straight forward + noOfValues = 0; + formants = []; + while( ~(feof(FileID)) ) + + % search for numberOfFormants + thestr = fgetl(FileID);%, '%s', 1); + + if( strfind( thestr , 'numberOfFormants' ) > 0 ) + noOfValues = noOfValues + 1; + %numberOfFormants found + pos = find( thestr == '=' ); + numberOfFormants = str2num(thestr(pos+2:end)); + formants( noOfValues, 1 ) = numberOfFormants; + % discard the 'formant []' line + thestr = fgetl(FileID); + % now read the formant positions + for (i=0:numberOfFormants-1) + thestr = fgetl(FileID); + pos = find( thestr == '=' ); + formants( noOfValues, i+2 ) = str2num(thestr( pos+2 : end )); + end + + % discard the 'bandwidth []' line + thestr = fgetl(FileID); + % now read the formant bandwidths + for (i=0:numberOfFormants-1) + thestr = fgetl(FileID); + pos = find( thestr == '=' ); + formants( noOfValues, i+2+numberOfFormants ) = str2num(thestr( pos+2 : end )); + end + + end + end + fclose(FileID); + end + + + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/Matlab/Common/ShortTimeEnergy.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,12 @@ +function E = ShortTimeEnergy(signal, windowLength, step); +signal = signal / max(max(signal)); +curPos = 1; +L = length(signal); +numOfFrames = floor((L-windowLength)/step) + 1; +%H = hamming(windowLength); +E = zeros(numOfFrames,1); +for (i=1:numOfFrames) + window = (signal(curPos:curPos+windowLength-1)); + E(i) = (1/(windowLength)) * sum(abs(window.^2)); + curPos = curPos + step; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/Matlab/Common/calculate_Silence.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,109 @@ +function [Limits] = calculate_Silence( x, fs, frameLength ) + +% Convert mono to stereo +if (size(x, 2)==2) + x = mean(x')'; +end +[m n]=size(x); +if(m>n) + x=x'; +end +% Window length and step (in seconds): +win = frameLength/fs; +step = win; + +%%%%%%%%%%%%%%%%%%%%%%%%%%% +% THRESHOLD ESTIMATION +%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Weight = 10; % used in the threshold estimation method + +% Compute short-time energy and spectral centroid of the signal: +Eor = ShortTimeEnergy(x, win*fs, step*fs); + + +% Apply median filtering in the feature sequences (twice), using 5 windows: +% (i.e., 250 mseconds) +E = medfilt1(Eor, 5); E = medfilt1(E, 5); + +% normalise + +E = E*(1/max(E)); + +% Get the average values of the smoothed feature sequences: +E_mean = mean(E); +% +%Find energy threshold: +% [HistE, X_E] = hist(E, round(length(E) / 10)); % histogram computation +% [MaximaE, countMaximaE] = findMaxima(HistE, 3); % find the local maxima of the histogram + + +T_E = 0.0005; % determined empirically + +% Thresholding: +Flags1 = (E>=T_E); + +flags = Flags1; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SPEECH SEGMENTS DETECTION +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +count = 1; +WIN = 0; +Limits = []; +while (count < length(flags)) % while there are windows to be processed: + % initilize: + curX = []; + countTemp = 1; + % while flags=1: + while ((flags(count)==1) && (count < length(flags))) + if (countTemp==1) % if this is the first of the current speech segment: + Limit1 = round((count-WIN)*step*fs)+1; % set start limit: + if (Limit1<1) Limit1 = 1; end + end + count = count + 1; % increase overall counter + countTemp = countTemp + 1; % increase counter of the CURRENT speech segment + end + + if (countTemp>1) % if at least one segment has been found in the current loop: + Limit2 = round((count+WIN)*step*fs); % set end counter + if (Limit2>length(x)) + Limit2 = length(x); + end + + Limits(end+1, 1) = Limit1; + Limits(end, 2) = Limit2; + end + count = count + 1; % increase overall counter +end + +%%%%%%%%%%%%%%%%%%%%%%% +% POST - PROCESS % +%%%%%%%%%%%%%%%%%%%%%%% + +% find the frame index of the start and end frames of non-silent segments +Limits = floor( Limits / frameLength ); +pos = find(Limits == 0); +if(pos) + Limits(pos) = 1; +end + +% A. MERGE OVERLAPPING SEGMENTS: +RUN = 1; +while (RUN==1) + RUN = 0; + for (i=1:size(Limits,1)-1) % for each segment + if (Limits(i,2)>=Limits(i+1,1)) + RUN = 1; + Limits(i,2) = Limits(i+1,2); + Limits(i+1,:) = []; + break; + end + end +end + +% B. Get final segments: +segments = {}; +for (i=1:size(Limits,1)) + segments{end+1} = x(Limits(i,1)*frameLength:Limits(i,2)*frameLength); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/Matlab/Common/calculate_VoicedUnvoicedDecision.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,103 @@ +function [ vuv ] = calculate_VoicedUnvoicedDecision( sampleName, x, fs, frameLength, noOfFrames, silentFrames ) + +DEBUG = 0; % do we want to plot the output? + +y = buffer( x, frameLength ); + +% % open original annotation file for PhD speech +% fileName = [sampleName '.txt']; +% % read metrics from file +% fileID = fopen( fileName ); +% data = fscanf( fileID, '%d', inf ); +% annotatedPitch = data(3:4:end); +% fclose( fileID ); + +% open HNR file +fileName = [sampleName '_HNR.txt']; +% read metrics from file +fileID = fopen( fileName ); +data = fscanf( fileID, '%f', inf ); +harmonic2noise = data(2:2:end); +fclose( fileID ); + +largestH2NValue = max( [max(harmonic2noise) abs(min(harmonic2noise))] ); +%normalise +harmonic2noise = harmonic2noise/largestH2NValue; +%round to 1dp +harmonic2noise= (round(harmonic2noise*10))/10; + +% open audio power file +fileName = [sampleName '_AP.txt']; +% read metrics from file +fileID = fopen( fileName ); +data = fscanf( fileID, '%f', inf ); +audioPower = data(2:2:end); +fclose( fileID ); +maxPower = max(audioPower); + + +minFreq = ceil(fs/getVariables('getMinFreq')); %default minimum fundatmental frequency +maxFreq = ceil(fs/getVariables('getMaxFreq')); %default maximum fundamental frequency + +% open pitch file +fileName = [sampleName '_YIN_pitch.txt']; +% read metrics from file +fileID = fopen( fileName ); +data = fscanf( fileID, '%f', inf ); +pitch = data(2:2:end); +fclose( fileID ); + +% open band power file +% fileName = [sampleName '_ASBP_norm.txt']; +% % read metrics from file +% fileID = fopen( fileName ); +% data = fscanf( fileID, '%f', inf ); +% bandPowerRatio = data(2:2:end); +% fclose( fileID ); +% bandPowerRatioRMS = sqrt( mean(( bandPowerRatio/max(bandPowerRatio) ).^2 )); +% + + +HNRThresh = 0; +if(DEBUG) hold off; end +minLength = min([length(harmonic2noise) length(audioPower) length(pitch)]); +for( i=1:minLength ) + + + if( silentFrames(i) == 0 ) % harmonic2noise(i) == 0 ) + %silent + if(DEBUG) plot( ((i-1)*frameLength)+1:(i*frameLength), y(:,i), 'y' );hold on; end; + vuv(i) = 0; + elseif( ((harmonic2noise(i) <= HNRThresh))... %|| (bandPowerRatio(i)/max(bandPowerRatio) < 0.001)) ... + || (pitch(i) > minFreq) ... + || (pitch(i) < maxFreq)) %... +% || (audioPower(i)/maxPower < 0.0005)) + % unvoiced + if(DEBUG) plot( ((i-1)*frameLength)+1:(i*frameLength), y(:,i), 'r' );hold on; end + vuv(i) = 2; + else + %voiced + if(DEBUG) plot( ((i-1)*frameLength)+1:(i*frameLength), y(:,i), 'g' );hold on; end + vuv(i) = 1; + end +end + +if(DEBUG) + + plot( frameLength/2:frameLength:(frameLength*length(harmonic2noise)), harmonic2noise, 'b' );hold on; +% plot( frameLength/2:frameLength:(frameLength*length(audioPower)), audioPower/maxPower, 'm' );hold on; + plot( frameLength/2:frameLength:(frameLength*length(pitch)), pitch/max(pitch), 'w' );hold on; +% plot( frameLength/2:frameLength:(frameLength*length(bandPowerRatio)), bandPowerRatio/max(bandPowerRatio), 'c' );hold on; + + L=line([0 length(x)], [maxFreq/max(pitch) maxFreq/max(pitch)]); + set(L,'color',[1 0 0]); + L=line([0 length(x)], [minFreq/max(pitch) minFreq/max(pitch)]); + set(L,'color',[1 0 0]); +% L=line([0 length(x)], [0.0001 0.0005]); +% set(L, 'color', [0 1 1] ); + %HNR threshold + L=line([0 length(x)], [HNRThresh HNRThresh]); + set(L, 'color', [0 0 1] ); +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/Matlab/Common/detect_Silence.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,37 @@ +function [ silenceSegments ] = detect_Silence( currentSampleName, OVERWRITE ); + +%returns an array containing the start and end frames for non-silent segments + +% open original silence calculation +fileName = [ currentSampleName '_silence.txt']; +% read pitch metrics from file +fileID = fopen( fileName ); + +if( (fileID <= 0) || (OVERWRITE) ) %does the file exist? + % no + disp('WARNING: MISSING SILENCE FILE'); + %calculate it + [x, fs, frameLength, noOfFrames] = openFile( [ sampleFileName '.wav' ] ); + limits = calculate_Silence( x, fs, frameLength ); + % create voicing metrics file + fileID = fopen( fileName, 'w'); + fprintf( fileID, 'non-silent start frame \t non-silent end frame \n '); + rowNum = size(limits); + rowNum = rowNum(1); + if( rowNum == 0 ) + % there are no silent frames + fprintf( fileID, '%d \t %d \n', 0, 0 ); + end + for i=1:rowNum + fprintf( fileID, '%d \t %d \n', limits(i,1), limits(i,2) ); + end + fclose( fileID ); + fileID = fopen( fileName ); + end +% +% strip off header +junk = fscanf( fileID, '%s', 6 ); +silence = fscanf( fileID, '%d', inf ); +silenceSegments = [silence(1:2:end) silence(2:2:end)]; +fclose( fileID ); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/Matlab/Common/detect_VoicedUnvoiced.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,36 @@ +function [vuv] = Detect_VoicedUnvoiced( sampleWavFileName, x, fs, frameLength, noOfFrames ) + +sampleFileName = sampleWavFileName( 1 : length( sampleWavFileName ) - 4 ); +fileName = [ sampleFileName '_VUV.txt']; +fileID = fopen( fileName ); +% +if( fileID < 0 ) + %---------------- GET THE SILENT FRAME VALUES ------------------- + + % only wish to consider pitch values from voiced frames. + % silent and unvoiced frames will produce pitch values that + % are random and therefore will bias our results + segmentFrames = detect_Silence( sampleFileName, 0 ); + + % remove the silent frames + [ silentFrames ] = removeSilentData( segmentFrames, noOfFrames ); + + % [vuv] = voicingByClustering( nonSilentAudio, fs, noOfFrames, frameLength ); + [vuv] = calculate_VoicedUnvoicedDecision( sampleFileName, x, fs, frameLength, noOfFrames, silentFrames ); + + noOfValidFrames = length(vuv); + +% plotAudioFrameByType( nonSilentAudio, noOfValidFrames, vuv, frameLength ); +% fclose( fileID ); + fileID = fopen( fileName, 'w'); + for i = 1 : length(vuv) + fprintf(fileID,'%d %d \n', i, vuv(i)); + end + fclose( fileID ); + fileID = fopen( fileName ); +end + +vuv = fscanf( fileID, '%d', inf ); +vuv = vuv( 2:2:end ); + +fclose( fileID ); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/Matlab/Common/detect_pitch.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,32 @@ +function [ pitch ] = pitch_Detection( sampleWavFileName, OVERWRITE ) + +% open original pitch calculation +sampleFileName = sampleWavFileName( 1 : length( sampleWavFileName ) - 4 ); +fileName = [ sampleFileName '_YIN_pitch.txt']; +% read pitch metrics from file +fileID = fopen( fileName ); + +if( (fileID <= 0) || (OVERWRITE) ) %does the file exist? + % no + disp('WARNING: MISSING PITCH FILE'); + %calculate it + [x, fs, frameLength, noOfFrames] = openFile( [ sampleFileName '.wav' ] ); + % perform pitch detection + framePrd = calculate_pitchYIN( sampleFileName, x, fs, frameLength ); + + fileID = fopen( fileName, 'w'); + for( i=1:length( framePrd )) + fprintf(fileID, '%d %f \n', framePrd(i,1), framePrd(i,2) ); + end + fclose( fileID ); + fileID = fopen( fileName ); +end + +pitch = fscanf( fileID, '%f', inf ); +pitch = pitch(2:2:end); +fclose( fileID ); + + + +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/Matlab/Common/removeSilentData.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,40 @@ +function [ silentFrames ] = removeSilentData( silentSegments, noOfFrames ) + +%---------------- GET THE SILENT FRAME VALUES ------------------- +% only wish to consider values from voiced frames. +% silent frames will produce data values that +% are random and therefore will bias our results +% returns an array of values where '1' indicates a non-silent frame + + silentFrames = zeros(1,noOfFrames); + [m n] = size(silentSegments); + if( silentSegments(1,1) == 0 && silentSegments(1,2) == 0 ) + % we have only one segment - the entire file + silentSegments(1,1) = 1; + silentSegments(1,2) = floor(length(x)/frameLength); + [m n] = size(silentSegments); + end + + % check for start at sample 0 + if( silentSegments(1,1) == 0) + silentSegments(1,1) = 1; + end + + for i=1:m + start = silentSegments(i,1); + stop = silentSegments(i,2); + if( stop > noOfFrames ) + stop = noOfFrames; + silentSegments(i,2) = noOfFrames; + end + + if( start == 0 ) + start = 1; + end + + + silentFrames( start : stop ) = 1; + end + +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/calculate_pitchYIN.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,51 @@ +function [framePrd] = calculate_pitchYIN( sampleFileName, x, fs, frameLength, noOfFrames ) + +% sampleFileName = sampleWavFileName( 1 : length( sampleWavFileName ) - 4 ); +% pitchFileName = [ sampleFileName '_YIN_pitch.txt']; +% pitchFileID = fopen( pitchFileName, 'w'); + +[periodInSamples, fileinfo] = yin ([ '../' sampleFileName '.wav']); + +% YIN produces one pitch estimate for every 32 samples. For our purposes +% this must be scaled to one sample every frame (approx 10ms but also a +% power of two). + + frameLength = 2^(ceil(log2(fileinfo.sr/100))); + noOfFrames = floor(fileinfo.nsamples / frameLength); %skip the final samples if necessary + +yinPrdSamplesPerFrame = length(periodInSamples) / noOfFrames; + + +start = 1; +stop = length(periodInSamples); +framePrd = []; + +for i=1:noOfFrames + % find the first and last periodInSamples frames to be averaged for + % every master frame of length frameLength + start = floor((i-1) * (length(periodInSamples) / noOfFrames)) + 1; + stop = ceil(i * (length(periodInSamples) / noOfFrames)); + if( stop <= length(periodInSamples) ) + % exclude NaN frames + prdChunk = periodInSamples( start : stop ); + goodPrdIdx = find( isnan( prdChunk )==0 ); + if( length( goodPrdIdx ) > 0 ) + framePrd(i,2) = mean( prdChunk(goodPrdIdx)); + else + framePrd(i,2) = NaN; + end + framePrd(i,1) = i; % record the frame number + % might need to know this in case of frame mismatch + prdFrames(i) = stop - start; + % write to file +% fprintf(pitchFileID, '%d %f \n', framePrd(i,1), framePrd(i,2) ); + else + break; + end +end + +% plot(periodInSamples,'x'); hold on; +% xaxisdiv=((length(periodInSamples)/length(framePrd))/2):(length(periodInSamples)/length(framePrd)):length(periodInSamples) + ((length(periodInSamples)/length(framePrd))/2); +% plot(xaxisdiv(1:length(framePrd)),framePrd(:,2),'rx'); + +% fclose( pitchFileID ); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/cumnorm_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,9 @@ +% cumnorm_inplace(x) - cumulative mean-normalization +% +% x: matrix to normalize +% +% Each column of x is normalized by dividing each sample by the mean of +% samples with indices smaller or equal to it. +% +% Mex function. +% Warning: the input argument is assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/dftoperiod.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,13 @@ +% prd=dftoperiod(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) is the period estimate. +% +% Mex function.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/dftoperiod2.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,17 @@ +% prd=dftoperiod2(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) +% is the period estimate. +% +% This version differs from dftoperiod in that the threshold is +% incremented by the global minimum of the difference function. +% +% Mex function.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/minparabolic.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,9 @@ +% [y,idx] = minparabolic(x) - locate minima using parabolic interpolation +% +% x: input vector +% y: same as x, with local minima replaced by interpolated minima +% idx: position of interpolated minima (fractionary) +% +% At each local minimum of x, that sample and its neighbors are used to +% determine a parabola. The minimum of the parabola is used as the new +% minimum value, the abscissa of the minimum is used as its position.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/rdiff_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,19 @@ +% rdiff_inplace(x,y,d,lags,n) - in place running cross-difference function +% +% x: column vector +% y: column vector +% r: result matrix (time X lag) +% lags: 2-column matrix of lags +% n: (samples) frame-rate & window size (default=1) +% +% Vectors x and y are each delayed by amounts specified by lags and subtracted +% sample to sample. The difference is squared and added over a time window +% of n samples. This processing is repeated every n samples, as many times +% as there are columns in r. +% +% A positive lag applied to x causes it to be delayed with respect to y. +% A positive lag applied to y causes it to be delayed with respect to x. +% The number of rows of r must match that of lags. +% +% Mex function. +% Beware: input arguments are assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/rsmooth.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,18 @@ +%y=rsmooth(x,smooth,npasses,trim) - smooth by running convolution +% +% X: input matrix +% SMOOTH: samples - size of square smoothing window +% NPASSES: number of smoothing passes (default=1) +% TRIM: if true, clip Y to same size as X +% +% Y: output matrix +% +% RSMOOTH smooths each column of matrix X by convolution with a square window +% followed by division by the window size. +% Multiple passes allow smoothing with a triangular window (npasses=2), or +% window shapes that approach a gaussian (npasses large). Convolution is +% implemented as a running sum for speed. +% +% Y has NPASSES*(SMOOTH-1) more rows than X unless TRIM is set. +% +% mex function
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/rsum_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +% rsum_inplace(x,N) - in place running sum +% +% x: matrix to sum +% N: samples - size of running sum window +% +% Matrix x is summed row-wise. The last (N-1) samples are not valid. +% N may be fractionary (handled by linear interpolation). +% +% Mex function. +% Beware: input argument is assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/sf_cleanup.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,29 @@ +function i = sf_cleanup(i) +% i=sf_cleanup(i) - cleanup after use (close i.fd if open) + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin ; error('sf_info: no input arguement'); end +if ~isa(i, 'struct') error('sf_info: expected struct'); end + +% close file if open +if isfield(i, 'fd') & fopen(i.fd) + fclose(i.fd); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/sf_format.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,213 @@ +function i = sf_format(i) +% i=sf_format(i) - try to guess file format from +% magic words or file suffix + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin; error('no argument'); end +if ~isa(i, 'struct') + i.fname=i; + i=sf_format(i); + return +end +if ~isfield(i, 'fname'); error('no fname field'); end + +if isa(i.fname, 'double') + i.format = 'matrix'; return +end + +if isa(i.fname, 'char') + disp(['guessing format of ', i.fname, '...']); +else + disp('guessing format...'); +end + +% Well-known file extensions +a = findstr('.', i.fname); +if a + suff = i.fname( a(size(a,2))+1 : size(i.fname,2)); + suff = lower(suff); + switch suff + case {'short', 'long', 'float', 'double', 'ascii', 'mat'} + i.format = suff; + return +% the following are commented out to exercise header magic +% case {'aiff'} +% i.format = 'AIFF'; +% return +% case {'aifc'} +% i.format = 'AIFC'; +% return +% case 'wav' +% i.format = 'WAV'; +% return +% case 'au' +% i.format = 'AU'; +% case 'sdif' +% i.format = 'sdif'; + % others ? + end +end + +% close file if open +if exist('i.fd') & fopen(i.fd) + fclose(i.fd); +end + + +% open little-endian, check for .wav magic +[i.fd, msg] = fopen(i.fname, 'r', 'l'); +if i.fd == -1 + error(['could not open: >', i.fname, '<']); + if ~isempty(msg) + error(msg); + end +end + +% file size in bytes +if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; +i.nbytes = ftell(i.fd); +if i.nbytes == 0 & strcmp(computer, 'MAC2') + disp('file is empty: perhaps its a macintosh SND resource?'); + i.format = 'MACSND'; % macintosh sound resource (maybe...) + return +end +if i.nbytes < 4; error('file is less than 4 bytes long'); end; +fseek(i.fd, 0, -1); +magic = char(fread(i.fd, 4, 'uchar'))'; +if strcmp(magic, 'RIFF') & i.nbytes >=12 + dummy = fread(i.fd, 1, 'ulong'); % chunk size + magic = char(fread(i.fd, 4, 'uchar'))'; + if strcmp(magic, 'WAVE'); + i.format='WAV'; + return + end +end + +% reopen for standard binary, look for magic +fclose(i.fd); +i.fd = fopen(i.fname, 'r'); +if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; +i.nbytes = ftell(i.fd); +fseek(i.fd, 0, -1); +% magic strings? +magic = char(fread(i.fd, 4, 'uchar'))'; +switch magic +case '.snd' + i.format = 'AU'; + return; +case 'FORM' + fseek(i.fd, 8, -1); + i.format = char(fread(i.fd, 4, 'uchar'))'; + if ~strcmp(i.format, 'AIFF') & ~strcmp(char(i.format),'AIFC') + error (['expected AIFF or AIFC, found ' i.format]); + end; + return; +case 'NIST' + fseek(i.fd, 0, -1); + line = fscanf(i.fd, '%s' , 1); + if ~strcmp(line, 'NIST_1A'); + error(['expected NIST_1A, found ', magic]); + end; + disp(i.format); + return; +case 'SDIF' + i.format = 'SDIF'; + return; +end +% wff magic number? +wff_magic = 195894762; +fseek(i.fd, 0, -1); +magic = fread(i.fd, 1, 'long'); +if magic == wff_magic + i.format = 'wff'; % Nottingham IHR wff format + return +end +% ESPS? +esps_min_size = 333; esps_magic = 27162; +if i.nbytes > esps_min_size + fseek(i.fd, 16, -1); + magic = fread(i.fd, 1, 'uint'); + if magic == esps_magic + i.format = 'ESPS'; + return + end +end + + +% Magic didn't succeed, try to determine type based on file name +a = findstr('.', i.fname); +if a + suff = i.fname( a(size(a,2))+1 : size(i.fname,2)); + switch suff + case 's' + i.format = 'short'; + return + case 'f' + i.format = 'float'; + return + % others ? + end +end + + +% Ascii formats +fclose(i.fd); +i.fd = fopen(i.fname, 'rt'); +line = fscanf(i.fd, '%s' , 1); +if strcmp(line, '|WAVE'); + i.format = 'iwave'; % Nottingham IHR |WAVE format + return; +end + +% plain ascii ? +limit = 2048; % test only this many bytes +hunk = char(fread(i.fd, min(limit, i.nbytes))); +chars = unique(hunk); +if isempty(setdiff(chars, [9, 10, 13, ... + double(' abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890+-_.,:;')])) + % ascii + if ~isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.,')])) + % not just numerical: try removing first line + fseek(i.fd, 0, -1); + line = fgets(i.fd); + if -1 == line; return ; end % line too long + i.bytes_to_data = ftell(i.fd); + hunk = char(fread(i.fd, min(limit, i.nbytes-ftell(i.fd)))); + chars = unique(hunk); + warning('skipping first line of what seems like an ascii file'); + end + if isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.')])) + i.format = 'ascii'; + return + end + if isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.,')])) + w.format = 'csv'; + return + end +end + +% Default +if i.nbytes == 2 * round(i.nbytes/2); % even number of bytes + i.format = 'short'; +else + i.format = 'uchar'; +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/sf_info.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,438 @@ +function i = sf_info(i) +% i=sf_info(i) - extract useful info from file + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin ; error('sf_info: no input arguement'); end +if ~isa(i, 'struct') + j.fname=i; + i=j; + i=sf_info(i); + return +end +if ~isfield(i, 'fname'); error('sf_info: no fname field'); end + +% guess format only if unknown +if ~isfield(i, 'format') | isempty(i.format); + i = sf_format(i); + if ~strcmp(i.format, 'matrix') disp(i.format); end +end + +% handle workspace matrices as if they were files +if strcmp(i.format, 'matrix') + [nrows, ncols] = size(i.fname); + i.nchans=ncols; + i.nsamples=nrows; + i.totalsamples = i.nsamples*i.nchans; + i.sr=[]; + return; +end + +% close file if open +if exist('i.fd') & fopen(i.fd) + fclose(i.fd); +end + +% use standard matlab functions for AU and WAV and MACSND +if strcmp(i.format, 'AU') + if isempty(findstr('.', i.fname)) + disp(['WARNING: matlab function AUREAD requires '... + '.au or .snd suffix on file name: ' i.fname]); + end + sz = auread(i.fname, 'size'); + i.nsamples=sz(1); + i.nchans=sz(2); + i.totalsamples = i.nsamples*i.nchans; + [dummy, i.sr, i.samplebits] = auread(i.fname, 1); + return; +end +if strcmp(i.format, 'WAV') + if isempty(findstr('.wav', i.fname)) & isempty(findstr('.WAV', i.fname)) + disp(['WARNING: matlab function WAVREAD requires '... + '.wav suffix on file name: ' i.fname]); + end + sz = wavread(i.fname, 'size'); + i.nsamples=sz(1); + i.nchans=sz(2); + i.totalsamples = i.nsamples*i.nchans; + [dummy, i.sr, i.samplebits] = wavread(i.fname, 1); + return; +end +if strcmp(i.format, 'MACSND') + % must load the data to get info - this is stupid + if ~isempty(findstr(':', i.fname)) + disp(['matlab function READSND cannot handle an '... + 'indirect path: ' i.fname]); + end + if 3==exist('readsnd') + [data, i.sr] = eval('readsnd(i.fname)'); + else + error('cannot read MACSND on this platform'); + end + i.nsamples = size(data,2); + i.nchans = size(data,1); + i.totalsamples = i.nsamples*i.nchans; + return +end + +% reopen file +if strcmp(i.format, 'ascii') ... + | strcmp(i.format, 'csv') | strcmp(i.format, '|WAVE') + [i.fd, msg] = fopen(i.fname, 'rt'); +else + [i.fd, msg] = fopen(i.fname, 'r', 'ieee-be.l64'); +end +if i.fd == -1 + if isempty(msg) + error(['could not open: ', i.fname]); + else + error(msg) + end +end +fd = i.fd; +if ~isfield(i, 'nbytes') + if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; + i.nbytes = ftell(i.fd); +end +fseek(fd, 0, -1); % rewind + + +switch i.format +%%%%%%%%%%%%%%%%%%% AIFF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case {'AIFF','AIFC'} + fseek(fd, 12, 0); % skip container chunk + % skip over spurious chunks + idx=ftell(fd); + while 1 + magic=char(fread(fd,4,'uchar'))'; + if strcmp(magic,'COMM'); break; end; + idx = idx+1; + status = fseek(fd,idx,-1); + if status == -1 + error('expected COMM magic word, found eof'); + end; + end; + + %ckSize=fread(fd,1,'int32'); + %status = fseek(fd, ckSize, 0); % skip to end of chunk + %if status == -1 + % error('unexpected eof'); + %end + + %while (1) + % magic = char(fread(fd, 4, 'char'))'; + % if ~strcmp(magic, 'SSND') + % fseek(fd, -4, 0); % skip back + % break; + % end + % ckSize = fread(fd, 1, 'long'); + % fseek(fd, ckSize, 0); % skip to end of sound chunk + %end + %magic = char(fread (fd, 4, 'char'))'; + %if ~strcmp(magic, 'COMM'); error(['expected COMM, found ', magic]) ; end + commsz = fread(fd, 1, 'int32'); + i.nchans = fread(fd, 1, 'int16'); + i.nsamples = fread(fd, 1, 'uint32'); + i.totalsamples = i.nsamples*i.nchans; + i.samplebits = fread(fd, 1, 'int16'); + switch i.samplebits + case 16 + i.sample_bytes = 2; + i.sample_type = 'int16'; + case 32 + i.sample_bytes = 4; + i.sample_type = 'int32'; % or float? + otherwise + error(['unexpected samplebits: ' num2str(i.samplebits) ]); + end + % read sampling rate using Hideki Kawahara's code: + srex1=fread(fd,1,'uint16'); + srex2=fread(fd,1,'uint64'); + if strcmp(char(i.format),'AIFC') + compress=fread(fd,4,'uchar'); + if ~strcmp(char(compress),'NONE') + error('Compression is not supported.'); + end; + fseek(fd, commsz-22, 0); + end; + i.sr = 2^(srex1-16383)*srex2/hex2dec('8000000000000000'); + %fseek(fd, 12, -1); % skip back to end of container chunk + % skip over eventual common chunk + %while(1) + % magic = char(fread(fd, 4, 'char'))'; + % if ~strcmp(magic, 'COMM') + % fseek(fd, -4, 0); % skip back + % break; + % end + % ckSize = fread(fd, 1, 'long'); + % fseek(fd, ckSize, 0); % skip over chunk + %end + magic=char(fread(fd,4,'uchar'))'; + while ~strcmp(char(magic),'SSND') + [ckSize, count]=fread(fd,1,'int32'); + if ~count; + error('expected chunk size field, found eof'); + return + end + status = fseek(fd, ckSize, 0); % skip to end of chunk + if status == -1 + error('expected SSND magic word, found eof'); + return + end; + magic=char(fread(fd,4,'uchar'))'; + end; + + %magic = char(fread(fd, 4, 'char'))'; + %if ~strcmp(magic, 'SSND') + % error (['expected SSND, found' magic]); + %end + fseek(fd, 12, 0); % skip over ckSize, offset and blocksize fields + i.bytes_to_data = ftell(fd); + if i.totalsamples*i.sample_bytes ~= i.nbytes-i.bytes_to_data + disp(['WARNING: header fields sample_bytes: ' ... + num2str(i.sample_bytes)]); + disp (['and sample and channel count: ' num2str(i.nsamples) ... + ', ' num2str(i.nchans)]); + disp(['are inconsistent with offset to data: ' ... + num2str(i.bytes_to_data) ' and file size: ' ... + num2str(i.nbytes)]); + disp (['(' num2str(i.totalsamples*i.sample_bytes) ... + ' ~= ' num2str(i.nbytes-i.bytes_to_data) ')']); + end + + return +%%%%%%%%%%%%%%%%%%% NIST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'NIST' +% fseek(fd, 0, -1); +% line = fscanf(fd, '%s' , 1); +% if ~strcmp(line, 'NIST_1A'); error(['expected NIST_1A, found ', magic]); end + fseek(fd, 8, 0); % skip over magic string + i.bytes_to_data = fscanf(fd, '%d', 1); + while (1) + key = fscanf(fd, '%s', 1); + if strcmp(key, 'end_head'); break; end; + % read third field according to spec in second field (type) + % this may need refining... + type = fscanf(fd, '%s', 1); + if strcmp(type(1:2), '-s') + bytes_to_read = sscanf(type(3:end), '%d', 1); + fseek(fd, 1, 0); % skip blank + value = char(fread(fd, bytes_to_read, 'char'))'; + else + value = fscanf(fd, '%f', 1); + end + i = setfield(i, key, value); + end + % give standard names to useful fields + if isfield(i, 'channel_count'); i.nchans = i.channel_count; end + if isfield(i, 'sample_count'); i.nsamples = i.sample_count; end + i.totalsamples = i.nsamples*i.nchans; + if isfield(i, 'sample_rate'); + i.sr = i.sample_rate; + i.xunits = 's'; + end + if ~isfield(i, 'sample_coding'); i.sample_coding = 'pcm'; end + i.bytes_to_data = 1024; % needs checking + i.sample_bytes=2; + i.sample_type='int16'; + return +%%%%%%%%%%%%%%%%%%% |WAVE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case '|WAV' + line = fscanf(fd, '%s' , 1); % skip first line + i.nsamples = fscanf(fd, '%d', 1); + i.nchans = fscanf(fd, '%d', 1); + i.totalsamples = i.nsamples*i.nchans; + i.bytes_to_data = ftell(fd); + % channel info handled during data read + i.totalsamples = i.nsamples * i.nchans; + return +%%%%%%%%%%%%%%%%%%% WFF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'wff' + fseek(fd, 4, 0); % skip magic number + i.version = fread(fd, 1, 'long'); + i.type_code = fread(fd, 1, 'long'); + i.nchans = fread(fd, 1, 'long'); + i.info.channel_flags = fread(fd, 1, 'long'); + i.bytes_to_data = fread(fd, 1, 'long'); + fseek(fd, 40, 0); + i.gen_prog_name = fread(fd, 32, 'char'); + i.comment = fread(fd, 32, 'char'); + i.sample_bytes = 2; + i.sample_type = 'int16'; + return; +% channel info handled later during channel read +%%%%%%%%%%%%%%%%%%% ESPS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'ESPS' +% Based on Peter Kabal's afsp package. +% This handles at least one kind of ESPS waveform file. Others? + i.machine_code = fread(fd, 1, 'uint'); + i.version_check_code = fread(fd, 1, 'uint'); + i.bytes_to_data = fread(fd, 1, 'uint'); + i.record_size = fread(fd, 1, 'uint'); + fseek(fd, 20, -1); + i.EDR_ESPS_flag = fread(fd, 1, 'uint'); + i.align_pad_size = fread(fd, 1, 'uint'); + fseek(fd, 32, -1); + i.file_type = fread(fd, 1, 'uint16'); + fseek(fd, 40, -1); + i.file_creation_date_time = char(fread(fd, 26, 'char'))'; + i.header_version = char(fread(fd, 8, 'char'))'; + i.program_name = char(fread(fd, 16, 'char'))'; + i.program_version = char(fread(fd, 8, 'char'))'; + i.compile_date = char(fread(fd, 26, 'char'))'; + i.tag = fread(fd, 1, 'uint'); + fseek(fd, 132, -1); + i.ndoubles = fread(fd, 1, 'uint'); + i.nfloats = fread(fd, 1, 'uint'); + i.nlongs = fread(fd, 1, 'uint'); + i.nshorts = fread(fd, 1, 'uint'); + i.nchars = fread(fd, 1, 'uint'); + i.fixed_header_size = fread(fd, 1, 'uint'); + i.var_header_size = fread(fd, 1, 'uint'); + %w.dunno_what = fread(fd, 1, 'uint'); + fseek(fd, 160, -1); + i.user = char(fread(fd, 8, 'char'))'; + % scan the rest of the header to find sampling rate + a = ftell(fd); + bytes_left_in_header = i.bytes_to_data - ftell(fd); + hunk = char(fread(fd, bytes_left_in_header, 'uchar'))'; + b = findstr(hunk, 'record_freq'); + %fseek(fd, a+b-1, -1); + %w.xxxx = char(fread(fd, 12, 'char'))'; + %w.count = fread(fd, 1, 'uint'); + %w.data_code = fread(fd, 1, 'ushort'); + fseek(fd, a+b-1+12+2+4, -1); + i.sr = fread(fd, 1, 'double'); + i.nchans = i.ndoubles+i.nfloats+i.nlongs ... + +i.nshorts+i.nchars; + recordsize = i.ndoubles*8 + i.nfloats*4 + i.nlongs*4 ... + +i.nshorts*2+i.nchars; + i.nsamples = (i.nbytes - i.bytes_to_data) / recordsize; + i.totalsamples = i.nsamples * i.nchans; + i.sample_bytes = 2; % bug + i.sample_type = 'int16'; % bug + return; +%%%%%%%%%%%%% PLAIN ASCII %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'ascii' + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + line = fgetl(fd); + i.nchans = size(sscanf(line, '%f'), 1); + nlines = 1; + while (1) + line = fgets(i.fd); + if isa(line, 'double') & -1 == line; break; end + nlines = nlines+1; + end + i.nsamples = nlines; + i.totalsamples = i.nsamples*i.nchans; + i.sr=1; + return +%%%%%%%%%%%%% COMMA-SEPARATED VALUES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'csv'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + % todo + return +%%%%%%%%%%%%% Binary %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case {'uchar'}; % should take care of signed/unsigned + i.sample_bytes = 1; + i.sample_type = 'uchar'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'short', 'int16'}; + i.sample_bytes = 2; + i.sample_type = 'int16'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'long', 'int32'}; + i.sample_bytes = 4; + i.sample_type = 'int32'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'float', 'float32'}; + i.sample_bytes = 4; + i.sample_type = 'float32'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'double', 'float64'}; + i.sample_bytes = 8; + i.sample_type = 'float64'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +otherwise + error(['unknown format: >' i.format '<']); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/sf_wave.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,193 @@ +function y = sf_wave(i, samples, chans) +% y=sf_wave(i,samples,chans) - read data from file +% +% y: data read from file (column vector or matrix) +% +% i: structure containing information about file +% samples: range of samples to read ([start stop]) - default ([]) means entire file +% chans: range of channels to read ([lo hi]) - default ([]) means all channels\ +% + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin | ~isfield(i, 'fname') + error('usage: i.fname = "file name", then call y=sf_wave(i)'); +end +if ~isfield(i, 'format') + i = sf_info(i); +end +% defaults +if nargin<2 | isempty(samples); + if (i.nsamples) samples=[1 i.nsamples]; end; +end +if nargin<3 | isempty(chans) + if (i.nchans) chans=[1 i.nchans]; end; +end + +% clip +if samples(1) < 1 + samples(1) = 1; +end +if samples(2) > i.nsamples + samples(2) = i.nsamples; +end +if samples(2) < samples(1) + y=[];return + error(['start sample after stop sample: ' num2str(samples(1)) '-' num2str(samples(2))]); +end +if chans(1) < 1 | chans(2) > i.nchans + error(['requested inexistent channels: ' num2str(chans(1)) '-' num2str(chans(2))]); +end + + +% workspace matrix +if strcmp(i.format, 'matrix') + y = i.fname(samples(1):samples(2), chans(1):chans(2)); + return +end + +% use matlab functions for AU and WAV and MACSND +if strcmp(i.format, 'AU') + y=auread(i.fname, [samples(1) samples(2)]); + y=y(:,chans(1):chans(2)); + return; +end +if strcmp(i.format, 'WAV') + y=wavread(i.fname, [samples(1) samples(2)]); + y=y(:,chans(1):chans(2)); + return; +end +if strcmp(i.format, 'MACSND') + if 3==exist('readsnd') + y = eval('readsnd(i.fname)'); + else + error('cannot read MACSND on this platform'); + end + y = y(samples(1):samples(2),chans(1):chans(2)); + return; +end + +% close if open +% if fopen(i.fd) +% fclose(i.fd); +% end + +if ~isfield(i, 'bytes_to_data') + i.bytes_to_data=0; +end + +% ascii formats +if strcmp(i.format, 'ascii') | strcmp(i.format, 'csv') | strcmp(i.format, 'IWAVE') + i.fd = fopen(i.fname, 'rt'); + fseek(i.fd, i.bytes_to_data, -1); + + switch i.format + case 'ascii' + nsamples = samples(2) - samples(1) + 1; + nchans = chans(2) - chans(1) + 1; + y = zeros(nsamples, nchans); + % skip to start + for j=1:samples(1)-1 + line = fgetl(i.fd); + if isempty(line) | line == -1; error('unexpected eof'); end + end + k=1; + % read + for j=samples(1) : samples(2) + line = fgetl(i.fd); + if isempty(line) | line == -1; error('unexpected eof'); end + a = sscanf(line, '%f'); + y(k,:) = a(chans(1):chans(2))'; + k = k+1; + end + case 'cvs' + error('not implemented'); + case 'IWAVE' + error('not implemented'); + end + fclose(i.fd); + return +end + +% binary formats +fr = samples(2) - samples(1) + 1; +skip_samples = i.nchans * (samples(1) - 1); +switch i.format +case 'uchar' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 1, 0); + y = fread(i.fd, [fr, i.nchans], 'uchar'); + fclose(i.fd); +case 'short' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + fclose(i.fd); +case 'long' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 4, 0); + y = fread(i.fd, [fr, i.nchans], 'long'); + fclose(i.fd); +case 'float' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 4, 0); + y = fread(i.fd, [fr, i.nchans], 'float'); + fclose(i.fd); +case 'double' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 8, 0); + y = fread(i.fd, [fr, i.nchans], 'double'); + fclose(i.fd); +case 'NIST' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + y = zeros(i.nsamples, i.nchans); + switch i.sample_coding + case 'pcm' + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + otherwise + error(['cannot handle NIST sample_coding = ', i.sample_coding]); + end + fclose(i.fd); +case 'ESPS' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + +case 'AIFF' + i.fd = fopen(i.fname, 'r', 's'); + fseek(i.fd, i.bytes_to_data, -1); + % should check sample size + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + fclose(i.fd); + +otherwise + error(['don''t know how to load format = ', i.format]); +end + +y = y(:,chans(1):chans(2));
Binary file Code/Descriptors/yin/private/src/cumnorm_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
Binary file Code/Descriptors/yin/private/src/cumnorm_inplace.%b5 Data/RESOURCE.FRK/MATLAB~1.TDT has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/cumnorm_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,92 @@ +/* +* cumnorm_inplace.c +* cumulative mean-normalization of diff function +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] + +/* Output Arguments */ + +static void cumsum_inplace( + double *xp, /* matrix to cumulative-mean-normalize (along columns) */ + int m, /* rows */ + int n /* columns */ + ) +{ + int j,k; + double z, sum, mean; + + + for (k=0; k<n; k++) { /* for each column */ + SET(xp[k*m+0])=1; + sum = 0.0; + for (j=1;j<m;j++) { /* for each row */ + z = GET(xp[k*m+j]); + z = (z > 0) ? z:0; /* clip to remove numerical artifacts */ + sum += z; + z = (sum>0) ? (z / (sum/j)) : 1; /* cumulative-mean-normalize */ + SET(xp[k*m+j])=z; + } + } +/* + for (j=0;j<m;j++) { // for each row (time) + sum=0; + SET(xp[j])=1; + for (k=1; k<n; k++) { // for each column (lag) + z=GET(xp[m*k+j]); + sum+=z; + mean=sum/k; + SET(xp[m*k+j]) = (mean > 0) ? z/mean : 1; + } + } +*/ + + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp; + int nx, mx; + + /* Check for proper number of arguments */ + if (nrhs != 1) { + mexErrMsgTxt("CUMSUM_INPLACE takes 1 input argument"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("CUMSUM_INPLACE: X should be doubles"); + } + mx=mxGetM(X_IN); /* rows */ + nx=mxGetN(X_IN); /* columns */ + if (nx*mx == 0) { + mexErrMsgTxt("CUMSUM_INPLACE: input matrix is empty"); + } + + xp = mxGetPr(X_IN); + checkin_matrix((mxArray *) X_IN); + + + /* Do the actual computations in a subroutine */ + cumsum_inplace(xp,mx,nx); + + checkout_matrix((mxArray *) X_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/cumnorm_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,9 @@ +% cumnorm_inplace(x) - cumulative mean-normalization +% +% x: matrix to normalize +% +% Each column of x is normalized by dividing each sample by the mean of +% samples with indices smaller or equal to it. +% +% Mex function. +% Warning: the input argument is assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/cumnorm_inplace.mexw64.manifest Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +<?xml version='1.0' encoding='UTF-8' standalone='yes'?> +<assembly xmlns='urn:schemas-microsoft-com:asm.v1' manifestVersion='1.0'> + <trustInfo xmlns="urn:schemas-microsoft-com:asm.v3"> + <security> + <requestedPrivileges> + <requestedExecutionLevel level='asInvoker' uiAccess='false' /> + </requestedPrivileges> + </security> + </trustInfo> +</assembly>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/cumnorm_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod copy.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,14 @@ +% prd=dftoperiod(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) +% is the period estimate. +% +% Mex function.
Binary file Code/Descriptors/yin/private/src/dftoperiod.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,144 @@ +/* +* dftoperiod.c +* estimate period from df +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define D_IN prhs[0] +#define B_IN prhs[1] +#define T_IN prhs[2] + +/* output arguments */ +#define PRD_OUT plhs[0] + +static void ddftoperiods( + double *dp, /* input df */ + double *bp, /* bounds */ + double *prdp, /* periods */ + double thresh, /* threshold */ + int m, /* nrows of df */ + int n /* ncols of df */ + ) +{ + int j,k, p, lo, hi, maxj, goodflag; + double z, min, globalmin; + + /* bounds */ + lo= GET(bp[0]); + hi= GET(bp[1]); + if (lo<0 || hi>m ) { + mexPrintf("%d %d\n", lo, hi); + mexErrMsgTxt("DFTOPERIOD: bounds are out of bounds"); + } + + #define MARGIN 1.5 + + for (k=0;k<n;k++) { /* for each col */ + + goodflag=0; + + /* estimate global min + globalmin=GET(dp[k*m+lo]); + for (j=lo;j<hi; j++) { + z = GET(dp[k*m+j]); + if (z<globalmin) { + globalmin = z; + } + } + + thresh += globalmin; */ + + /* estimate period */ + p=lo; + min=GET(dp[k*m+lo]); + maxj=hi; + for (j=lo;j<maxj; j++) { /* for each row */ + + z = GET(dp[k*m+j]); + + if (z<thresh && z <= min) { /* good candidate, restrict search range */ + maxj=j*MARGIN; + if (maxj> hi) { maxj = hi; } + } + + if (z<min) { /* update min & period estimate */ + min = z; + p=j; + } + + } + + /* mexPrintf("%f %d %f\n", thresh, p, min); */ + + SET(prdp[k])=p; + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *dp, *bp, *tp, *prdp, thresh; + int m, n; + + /* Check for proper number of arguments */ + if (nrhs != 3){ + mexErrMsgTxt("DFTOPERIOD requires 3 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(D_IN) || mxIsComplex(D_IN) || + mxIsSparse(D_IN) || !mxIsDouble(D_IN) ) { + mexErrMsgTxt("DFTOPERIOD: D should be doubles"); + } + if (!mxIsNumeric(B_IN) || mxIsComplex(B_IN) || + mxIsSparse(B_IN) || !mxIsDouble(B_IN) ) { + mexErrMsgTxt("DFTOPERIOD: B should be doubles"); + } + if (!mxIsNumeric(T_IN) || mxIsComplex(T_IN) || + mxIsSparse(T_IN) || !mxIsDouble(T_IN) ) { + mexErrMsgTxt("DFTOPERIOD: T should be double"); + } + m=mxGetM(D_IN); /* rows */ + n=mxGetN(D_IN); /* columns */ + + if (mxGetM(B_IN) * mxGetN(B_IN) != 2) { + mexErrMsgTxt("DFTOPERIOD: B should be a size 2 vector"); + } + if (mxGetM(T_IN)*mxGetN(T_IN) != 1) { + mexErrMsgTxt("DFTOPERIOD: T should be scalar"); + } + + /* Create matrix to return */ + PRD_OUT = mxCreateDoubleMatrix(1,n, mxREAL); + + dp = mxGetPr(D_IN); + bp = mxGetPr(B_IN); + tp = mxGetPr(T_IN); thresh=tp[0]; + prdp = mxGetPr(PRD_OUT); + + checkin_matrix((mxArray *) D_IN); + checkin_matrix((mxArray *) B_IN); + checkin_matrix((mxArray *) PRD_OUT); + + /* Do the actual computations in a subroutine */ + ddftoperiods(dp,bp,prdp,thresh,m,n); + + checkout_matrix((mxArray *) D_IN); + checkout_matrix((mxArray *) B_IN); + checkout_matrix((mxArray *) PRD_OUT); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,13 @@ +% prd=dftoperiod(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) is the period estimate. +% +% Mex function.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod.mexw64.manifest Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +<?xml version='1.0' encoding='UTF-8' standalone='yes'?> +<assembly xmlns='urn:schemas-microsoft-com:asm.v1' manifestVersion='1.0'> + <trustInfo xmlns="urn:schemas-microsoft-com:asm.v3"> + <security> + <requestedPrivileges> + <requestedExecutionLevel level='asInvoker' uiAccess='false' /> + </requestedPrivileges> + </security> + </trustInfo> +</assembly>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/private/src/dftoperiod2.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod2.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,144 @@ +/* +* dftoperiods.c +* estimate period from df +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define D_IN prhs[0] +#define B_IN prhs[1] +#define T_IN prhs[2] + +/* output arguments */ +#define PRD_OUT plhs[0] + +static void ddftoperiods( + double *dp, /* input df */ + double *bp, /* bounds */ + double *prdp, /* periods */ + double thresh, /* threshold */ + int m, /* nrows of df */ + int n /* ncols of df */ + ) +{ + int j,k, p, lo, hi, maxj, goodflag; + double z, min, globalmin; + + /* bounds */ + lo= GET(bp[0]); + hi= GET(bp[1]); + if (lo<0 || hi>m ) { + mexPrintf("%d %d\n", lo, hi); + mexErrMsgTxt("DFTOPERIOD: bounds are out of bounds"); + } + + #define MARGIN 1.5 + + for (k=0;k<n;k++) { /* for each col */ + + goodflag=0; + + /* estimate global min */ + globalmin=GET(dp[k*m+lo]); + for (j=lo;j<hi; j++) { /* for each row */ + z = GET(dp[k*m+j]); + if (z<globalmin) { /* update min & period estimate */ + globalmin = z; + } + } + + thresh += globalmin; + + /* estimate period */ + p=lo; + min=GET(dp[k*m+lo]); + maxj=hi; + for (j=lo;j<maxj; j++) { /* for each row */ + + z = GET(dp[k*m+j]); + + if (z<thresh && z <= min) { /* good candidate, restrict search range */ + maxj=j*MARGIN; + if (maxj> hi) { maxj = hi; } + } + + if (z<min) { /* update min & period estimate */ + min = z; + p=j; + } + + } + + /* mexPrintf("%f %d %f\n", thresh, p, min); */ + + SET(prdp[k])=p; + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *dp, *bp, *tp, *prdp, thresh; + int m, n; + + /* Check for proper number of arguments */ + if (nrhs != 3){ + mexErrMsgTxt("DFTOPERIOD requires 3 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(D_IN) || mxIsComplex(D_IN) || + mxIsSparse(D_IN) || !mxIsDouble(D_IN) ) { + mexErrMsgTxt("DFTOPERIOD: D should be doubles"); + } + if (!mxIsNumeric(B_IN) || mxIsComplex(B_IN) || + mxIsSparse(B_IN) || !mxIsDouble(B_IN) ) { + mexErrMsgTxt("DFTOPERIOD: B should be doubles"); + } + if (!mxIsNumeric(T_IN) || mxIsComplex(T_IN) || + mxIsSparse(T_IN) || !mxIsDouble(T_IN) ) { + mexErrMsgTxt("DFTOPERIOD: T should be double"); + } + m=mxGetM(D_IN); /* rows */ + n=mxGetN(D_IN); /* columns */ + + if (mxGetM(B_IN) * mxGetN(B_IN) != 2) { + mexErrMsgTxt("DFTOPERIOD: B should be a size 2 vector"); + } + if (mxGetM(T_IN)*mxGetN(T_IN) != 1) { + mexErrMsgTxt("DFTOPERIOD: T should be scalar"); + } + + /* Create matrix to return */ + PRD_OUT = mxCreateDoubleMatrix(1,n, mxREAL); + + dp = mxGetPr(D_IN); + bp = mxGetPr(B_IN); + tp = mxGetPr(T_IN); thresh=tp[0]; + prdp = mxGetPr(PRD_OUT); + + checkin_matrix((mxArray *) D_IN); + checkin_matrix((mxArray *) B_IN); + checkin_matrix((mxArray *) PRD_OUT); + + /* Do the actual computations in a subroutine */ + ddftoperiods(dp,bp,prdp,thresh,m,n); + + checkout_matrix((mxArray *) D_IN); + checkout_matrix((mxArray *) B_IN); + checkout_matrix((mxArray *) PRD_OUT); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod2.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,17 @@ +% prd=dftoperiod2(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) +% is the period estimate. +% +% This version differs from dftoperiod in that the threshold is +% incremented by the global minimum of the difference function. +% +% Mex function.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod2.mexw64.manifest Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +<?xml version='1.0' encoding='UTF-8' standalone='yes'?> +<assembly xmlns='urn:schemas-microsoft-com:asm.v1' manifestVersion='1.0'> + <trustInfo xmlns="urn:schemas-microsoft-com:asm.v3"> + <security> + <requestedPrivileges> + <requestedExecutionLevel level='asInvoker' uiAccess='false' /> + </requestedPrivileges> + </security> + </trustInfo> +</assembly>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/dftoperiod2.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/private/src/interp_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/interp_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,101 @@ +/* +* interp_inplace.c +* linear interpolation +* +* Alain de Cheveigné, CNRS/Ircam +* (c) 2003 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] +#define IDX_IN prhs[1] + +/* Output Arguments */ + +static void interp_inplace( + double *xp, /* input vector */ + double *idxp, /* index vector */ + int m, /* rows input vector */ + int midx /* rows index vector */ + ) +{ + int j,idxi; + double a,b,idx,idxf; + + + for (j=0;j<midx;j++) { /* for each index */ + idx=GET(idxp[j]); + idxi=floor(idx); + idxf=idx-idxi; + if (idxi<0) { + a=GET(xp[0]); + b=GET(xp[1]); + SET(idxp[j]) = a + idx*(b-a); + } else if (idxi>=m) { + a=GET(xp[m-2]); + b=GET(xp[m-1]); + SET(idxp[j]) = b + (idx-m+1)*(b-a); + } else { + a=GET(xp[idxi]); + b=GET(xp[idxi+1]); + SET(idxp[j]) = a + idxf*(b-a); + } + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *idxp; + int nx, mx, nidx, midx; + + /* Check for proper number of arguments */ + if (nrhs !=2 ) { + mexErrMsgTxt("INTERP_INPLACE takes 2 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("INTERP_INPLACE: X should be doubles"); + } + if (!mxIsNumeric(IDX_IN) || mxIsComplex(IDX_IN) || + mxIsSparse(IDX_IN) || !mxIsDouble(IDX_IN) ) { + mexErrMsgTxt("INTERP_INPLACE: Y should be doubles"); + } + mx=mxGetM(X_IN); /* rows */ + nx=mxGetN(X_IN); /* columns */ + midx=mxGetM(IDX_IN); + nidx=mxGetN(IDX_IN); + + if (nx>1 || nidx>1) { + mexErrMsgTxt("INTERP_INPLACE: X and IDX should be column vectors"); + } + if (mx<1) { + mexErrMsgTxt("INTERP_INPLACE: X should have at least two samples"); + } + + xp = mxGetPr(X_IN); + idxp = mxGetPr(IDX_IN); + + checkin_matrix((mxArray *) IDX_IN); + checkin_matrix((mxArray *) X_IN); + + /* Do the actual computations in a subroutine */ + interp_inplace(xp,idxp,mx,midx); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix((mxArray *) IDX_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/interp_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/private/src/mininrange.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/mininrange.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,119 @@ +/* +* minInRange.c +* for each sample, return index of the smallest value within an interval +* of that sample +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ + +#define X_IN prhs[0] +#define INTERVAL_IN prhs[1] + +/* Output Arguments */ + +#define IDX_OUT plhs[0] + +static void mininrange( + double *xp, + double *intp, + double *idxp, + unsigned int n + ) +{ + int h,k, idx, left, right, interval; + double min; + double max; + + + for (k=0; k<n; k++) { + interval=(int) GET(intp[k]); + + left = k - interval/2; + right = left+interval; + + if (left<0) left=0; /* clip */ + if (right>n) right=n; + + min = GET(xp[k]); + idx = k; + for (h=left;h<right; h++) { + if (GET(xp[h]) < min) { /* update min */ + min = GET(xp[h]); + idx = h; + } + } + SET(idxp[k])=idx+1; + } + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *idxp, *intp; + unsigned int m,n; + + /* Check for proper number of arguments */ + + if (nrhs != 2) { + mexErrMsgTxt("MININRANGE requires two input arguments"); + } else if (nlhs != 1) { + mexErrMsgTxt("MININRANGE requires one output argument"); + } + + /* Check type of input */ + + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN)) { + mexErrMsgTxt("MININRANGE: X should be matrix of doubles"); + } + if (!mxIsNumeric(INTERVAL_IN) || mxIsComplex(INTERVAL_IN) || + mxIsSparse(INTERVAL_IN) || !mxIsDouble(INTERVAL_IN)) { + mexErrMsgTxt("MININRANGE X should be matrix of doubles"); + } + + m = mxGetM(X_IN); /* rows */ + n = mxGetN(X_IN); /* columns */ + if (m>1 || n <=1) { + mexErrMsgTxt("MININRANGE X should be row vector"); + } + if (m != mxGetM(INTERVAL_IN) || n != mxGetN(INTERVAL_IN)) { + mexErrMsgTxt("MININRANGE: INTERVAL should be of same size as X"); + } + + /* Create matrix to return */ + + IDX_OUT = mxCreateDoubleMatrix(1, n, mxREAL); + + /* Assign pointers to the various parameters */ + + xp = mxGetPr(X_IN); + intp = mxGetPr(INTERVAL_IN); + idxp = mxGetPr(IDX_OUT); + + checkin_matrix((mxArray *) X_IN); + checkin_matrix((mxArray *) INTERVAL_IN); + checkin_matrix(IDX_OUT); + + /* Do the actual computations in a subroutine */ + + mininrange(xp,intp,idxp,n); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix((mxArray *) INTERVAL_IN); + checkout_matrix(IDX_OUT); + return; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/mininrange.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/private/src/minparabolic.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/minparabolic.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,121 @@ +/* +* minInterp.c +* determine position of minimum of parabolic fit to local minima +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] + + +/* Output Arguments */ +#define MIN_OUT plhs[0] +#define IND_OUT plhs[1] + + +static void mininterp( + double *xp, + unsigned int m, + unsigned int n, + double *minp, + double *indp + ) +{ + double x1,x2,x3,a,b,min,ind; + int i,j,k; + + /* each column is interpolated independently */ + /* each local minimum is replaced by the value of the interpolated minimum, and the + position of the interplolated minimum is recorded in indp */ + + for (j=0; j<n; j++) { + /* first row: don't touch */ + SET(indp[j*m])=0; + SET(minp[j*m])=GET(xp[j*m]); + /* intermediate: refine */ + for (k=1; k<m-1; k++) { + i=k+j*m; + x1=GET(xp[i-1]); x2=GET(xp[i]); x3=GET(xp[i+1]); + SET(indp[i]) = 0; /* default */ + SET(minp[i]) = x2; + if ((x2<=x1) && (x2<=x3)) { + a = (x1 + x3 - 2*x2)/2; + b = (x3 - x1)/2; + if (a) { + SET(indp[i]) = - b / (2*a); + SET(minp[i]) = x2 + - (b*b) / (4*a); + } + } + } + /* last row: don't touch */ + SET(indp[m-1 + j*m])=0; + SET(minp[m-1 + j*m])=GET(xp[m-1 + j*m]); + } + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *minp, *indp, *xp; + unsigned int m,n; + + /* Check for proper number of arguments */ + + if (nrhs != 1) { + mexErrMsgTxt("minParabolic requires one input argument."); + } else if (nlhs > 2) { + mexErrMsgTxt("minParabolic requires one or two output arguments"); + } + + + /* Check type of input */ + + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN)) { + mexErrMsgTxt("minParabolic requires that X be a matrix of doubles"); + } + + m = mxGetM(X_IN); /* rows */ + n = mxGetN(X_IN); /* columns */ + if (m<3) { mexErrMsgTxt("number of rows should be 3 or more"); } + + + /* Create two matrices to return */ + + MIN_OUT = mxCreateDoubleMatrix(m, n, mxREAL); + IND_OUT = mxCreateDoubleMatrix(m, n, mxREAL); + + + /* Assign pointers to the various parameters */ + + minp = mxGetPr(MIN_OUT); + indp = mxGetPr(IND_OUT); + xp = mxGetPr(X_IN); + + checkin_matrix((mxArray *) X_IN); + checkin_matrix(MIN_OUT); + checkin_matrix(IND_OUT); + + /* Do the actual computations in a subroutine */ + + + mininterp(xp,m,n,minp,indp); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix(MIN_OUT); + checkout_matrix(IND_OUT); + return; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/minparabolic.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,9 @@ +% [y,idx] = minparabolic(x) - locate minima using parabolic interpolation +% +% x: input vector +% y: same as x, with local minima replaced by interpolated minima +% idx: position of interpolated minima (fractionary) +% +% At each local minimum of x, that sample and its neighbors are used to +% determine a parabola. The minimum of the parabola is used as the new +% minimum value, the abscissa of the minimum is used as its position.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/minparabolic.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/mwrap.h Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,268 @@ +/* +19 Dec 94 - Sept 95 +Alain de Cheveigne, LLF, CNRS/Universite Paris 7. +Modified: Sept 2000 + +file mwrap.h + +To use the malloc wrappers: +- include this file in all source files that call malloc routines, +- link with mwrap.o, +- compile all files with -DMWRAP to turn on wrapping and checking. + +*/ + +#ifndef MWRAP_H +#define MWRAP_H + +#include <stdio.h> +#include <stdlib.h> + +/* NOTE: (char *) is used for pointer arithmetic. + Pointer arguments of all macros are cast to (char *). */ + +/* node structure for binary splay tree */ +typedef struct mwrap_node_struct *mwrap_node ; +struct mwrap_node_struct { + char *lo; /* base of block */ + char *hi; /* one beyond last address in block */ + mwrap_node up; + mwrap_node left; + mwrap_node right; +}; + + +/* Abbreviations for macros. Change in case of conflict. + + Assignment macros: + GET(x) equivalent to 'x' on a rhs, but memory address is checked before access + SET(x) equivalent to 'x' on a lhs, but memory address is checked before assignment + CPY(p,q,n) copy n values starting from p to n values starting from q, after checking both + + Test macros: + POK(p) pointer p is within an allocated block. + BOK(b) pointer b is the base of an allocated block. + PQOK(p, q) pointers p and q are in same allocated block. + PBOK(p, b) pointer p is in block of base b. + + MWRAPOK() check internal tree for sanity (slow) + + Query macros: + PBASE(p) base of block containing p, NULL if there is none. + PBUMP(p) one beyond last byte of block. + PSIZE(p) size of block containing p, -1 if none. + + Explicit block registering macros. + CHECKIN(lo, hi) register block of bounds lo, hi + CHECKOUT(lo) unregister block of base lo + +*/ + +#define GET(x) MWRAP_GET(x) +#define SET(x) MWRAP_SET(x) +#define CPY(p,q,n) MWRAP_CPY(p,q,n) +#define POK(p) MWRAP_POK(p) +#define BOK(b) MWRAP_BOK(b) +#define PQOK(p,q) MWRAP_PQOK((p),(q)) +#define PBOK(p,b) MWRAP_PBOK((p),(b)) +#define MWRAPOK() MWRAP_MWRAPOK() +#define PBASE(p) MWRAP_PBASE(p) +#define PSIZE(p) MWRAP_PSIZE(p) +#define CHECKIN(lo, hi) MWRAP_CHECKIN((lo), (hi)) +#define CHECKOUT(lo) MWRAP_CHECKOUT(lo) + + +/* Definitions of macros. + All are conditional on MWRAP being defined. + For speed, macros check for a "lucky" hit (block already + at root) before handing over to functions. */ + +#ifdef MWRAP + +#define MWRAP_GET(x) \ + ( (MWRAP_POK(&x)) ? (x) : (x) ) + +#define MWRAP_SET(x) \ + MWRAP_POK(&x) ; (x) + +#define MWRAP_CPY(p, q, n) \ + do {int i; \ + MWRAP_PQOK((p),(p)+(n-1)); \ + MWRAP_PQOK((q),(q)+(n-1)); \ + for(i=0;i<n;i++) (q)[i]=(p)[i]; \ + } while(0) + +/* #define MWRAP_POK(p) \ + do { if (mwrap_tree \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->hi) break; \ + mwrap_pok((char *) (p), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_POK(p) \ + ((mwrap_tree && (char *) (p) >= mwrap_tree->lo && (char *) (p) <mwrap_tree->hi) ? \ + 1 : mwrap_pok((char *) (p), __FILE__, __LINE__)) + + +/* #define MWRAP_BOK(b) \ + do { if (mwrap_tree \ + && (char *) (b) == mwrap_tree->lo) break; \ + mwrap_bok((char *) (b), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_BOK(b) \ + ((mwrap_tree && (char *) (b) == mwrap_tree->lo) ? 1 : mwrap_bok((char *) (b), __FILE__, __LINE__)) + + +/* #define MWRAP_PQOK(p, q) \ + do { if (mwrap_tree \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->hi \ + && (char *) (q) >= mwrap_tree->lo \ + && (char *) (q) < mwrap_tree->lo) break; \ + mwrap_pqok((char *)(p), (char *)(q), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_PQOK(p,q) \ + ((mwrap_tree \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->hi \ + && (char *) (q) >= mwrap_tree->lo \ + && (char *) (q) < mwrap_tree->lo ) ? \ + 1 : mwrap_pqok((char *)(p), (char *)(q), __FILE__, __LINE__)) + +/* #define MWRAP_PBOK(p, b) \ + do { if (mwrap_tree \ + && (char *) (b) == mwrap_tree->lo \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->lo) break; \ + mwrap_pbok((char *)(p), (char *)(b), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_PBOK(p, b) \ + (( mwrap_tree \ + && (char *) (b) == mwrap_tree->lo \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->lo) ? \ + 1 : mwrap_pbok((char *)(p), (char *)(b), __FILE__, __LINE__)) + +#define MWRAP_MWRAPOK() mwrap_checktree(__FILE__, __LINE__) + +#define MWRAP_PBASE(p) mwrap_pbase((char *) p) +#define MWRAP_PBUMP(p) mwrap_pbump((char *) p) +#define MWRAP_PSIZE(p) mwrap_psize((char *) p) +#define MWRAP_CHECKIN(lo, hi) \ + mwrap_checkin((char *) (lo), (char *) (hi), __FILE__, __LINE__) +#define MWRAP_CHECKOUT(lo) \ + mwrap_checkout( (char *) (lo), __FILE__, __LINE__) + +#else + +#define MWRAP_GET(x) (x) +#define MWRAP_SET(x) (x) +#define MWRAP_CPY(p,q,n) +#define MWRAP_POK(p) +#define MWRAP_BOK(p) +#define MWRAP_PBOK(p,b) +#define MWRAP_PQOK(p,q) +#define MWRAP_MWRAPOK() +#define MWRAP_PBASE(p) (void *) mwrap_squeal(__FILE__, __LINE__) +#define MWRAP_PBUMP(p) (void *) mwrap_squeal(__FILE__, __LINE__) +#define MWRAP_PSIZE(p) (long) mwrap_squeal(__FILE__, __LINE__) +#define MWRAP_CHECKIN(lo, hi) +#define MWRAP_CHECKOUT(lo) +#endif + + +/* Wrapping macros. + + If MWRAP is defined, all malloc routines are wrapped and the + original routines are available as "nowrap_malloc()", etc.. + + If MWRAP is not defined, malloc routines are not wrapped + and "nowrap_malloc()", etc. are #defined to "malloc()", etc. +*/ + +#ifdef MWRAP + +#ifndef MWRAP_C +#define malloc(size) mwrap_malloc((size), __FILE__, __LINE__) +#define free(p) mwrap_free((p), __FILE__, __LINE__) +#define realloc(p, size) mwrap_realloc((p), (size), __FILE__, __LINE__) +#define calloc(n, size) mwrap_calloc((n), (size), __FILE__, __LINE__) +#define cfree(p) mwrap_cfree((p), __FILE__, __LINE__) +#define valloc(size) mwrap_valloc((size), __FILE__, __LINE__) +#define vfree(p) mwrap_vfree((p), __FILE__, __LINE__) +#define memalign(a, p) mwrap_memalign((a), (p), __FILE__, __LINE__) +#endif /* MWRAP_C */ + +#else + +#ifndef MWRAP_C +#define nowrap_malloc(size) malloc(size) +#define nowrap_free(p) free(p) +#define nowrap_realloc(p, size) realloc(p) +#define nowrap_calloc(n, size) calloc((n), (size)) +#define nowrap_cfree(p) cfree(p) +#define nowrap_valloc(size) valloc(size) +#define nowrap_vfree(p) vfree(p) +#define nowrap_memalign(a, p) memalign((a), (p)) +#endif /* MWRAP_C */ + +#endif /* MWRAP */ + + +/* three global variables */ + +#ifndef MWRAP_C +extern mwrap_node mwrap_tree; /* the tree */ +extern char *mwrap_file; /* caller source file */ +extern long mwrap_line; /* line in caller source file */ +#endif + + +/* Data types used by malloc library. Modify if necessary */ +typedef void DATATYPE; /* returned by malloc(), etc. */ +typedef size_t SIZETYPE; /* fed to malloc(), etc. */ +typedef void VOIDTYPE; /* returned by free(), etc. */ + + +/* prototypes of routines made public */ + +/* malloc wrappers */ +DATATYPE *mwrap_malloc(SIZETYPE size, char *file, long line); +VOIDTYPE mwrap_free(DATATYPE *p, char *file, long line); +DATATYPE *mwrap_realloc(DATATYPE *p, SIZETYPE size, char *file, long line); +DATATYPE *mwrap_calloc(SIZETYPE n, SIZETYPE size, char *file, long line); +VOIDTYPE mwrap_cfree(DATATYPE *p, char *file, long line); +DATATYPE *mwrap_valloc(SIZETYPE size, char *file, long line); +VOIDTYPE mwrap_vfree(DATATYPE *p, char *file, long line); +DATATYPE *mwrap_memalign(SIZETYPE alignment, SIZETYPE size, char *file, long line); + +/* unwrapped versions */ +#ifdef MWRAP +DATATYPE *nowrap_malloc(SIZETYPE size); +VOIDTYPE nowrap_free(DATATYPE *p); +DATATYPE *nowrap_realloc(DATATYPE *p, SIZETYPE size); +DATATYPE *nowrap_calloc(SIZETYPE n, SIZETYPE size); +VOIDTYPE nowrap_cfree(DATATYPE *p); +DATATYPE *nowrap_valloc(SIZETYPE size); +VOIDTYPE nowrap_vfree(DATATYPE *p); +DATATYPE *nowrap_memalign(SIZETYPE alignment, SIZETYPE size); +#endif + +/* explicit block registering */ +void mwrap_checkin(char *lo, char *hi, char *file, long line); +void mwrap_checkout(char *lo, char *file, long line); + +/* tests & queries */ +int mwrap_pok(char *p, char *file, long line); +int mwrap_bok(char *b, char *file, long line); +int mwrap_pqok(char *p, char *q, char *file, long line); +int mwrap_pbok(char *p, char *b, char *file, long line); +char *mwrap_pbase(char *p); +char *mwrap_pbump(char *p); +long mwrap_psize(char *p); +long mwrap_squeal(char *file, long line); + +/* change default functions */ +void mwrap_set_errorfunction(void (*f)(int level)); +void mwrap_set_messagefunction(void (*f)(int level, char * message)); + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/mwrap_check.h Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,55 @@ +/* +* stuff to allow checking with mwrap +* +*/ + +#include "mwrap.h" +#include "string.h" + +#define MACINTOSH /* compiler (?) converts \n to \r: convert them back */ + +void checkin_matrix(mxArray *m); +void checkout_matrix(mxArray *m); +void mex_messagefunction(int level, char* message); +void mex_errorfunction(int level); + +/* check matrix into mwrap's tree */ +void checkin_matrix(mxArray *m) +{ + char *base, *top; + + base = (char *) mxGetPr(m); + top = base + mxGetN(m) * mxGetM(m) * sizeof(double); + /* mexPrintf("%d %d\n", base, top); */ + CHECKIN(base, top); +} + +/* check matrix out of mwrap's tree */ +void checkout_matrix(mxArray *m) +{ + char *base; + + base = (char *) mxGetPr(m); + CHECKOUT(base); +} + +/* message function to give to mwrap */ +void mex_messagefunction(int level, char* message) { +#ifdef MACINTOSH + {char *c; + c = strchr(message, (int) '\r'); + + while(c) { + *c = '\n'; + c = strchr(c, (int) '\r'); + } + } +#endif + mexPrintf(message); +} + +/* how to exit mex function */ +void mex_errorfunction(int level) { + /* mexPrintf("\n"); */ + mexErrMsgTxt(" "); +}
Binary file Code/Descriptors/yin/private/src/rdiff_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rdiff_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,170 @@ +/* +* rdiff_inplace.c +* running difference function +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] +#define Y_IN prhs[1] +#define R_IN prhs[2] +#define L_IN prhs[3] +#define N_IN prhs[4] + +/* Output Arguments */ + +static void rdiff_inplace( + double *xp, /* input vector 1 */ + double *yp, /* input vector 2 */ + double *rp, /* df matrix (time X lag) */ + double *lp, /* lag matrix */ + int N, /* size of summation window */ + int m, /* rows in corr matrix */ + int n /* columns in corr matrix */ + ) +{ + int j,k,i, xlag, ylag; + double z, d, *r, *x, *y, *bump; + + + for (j=0;j<n;j++) { /* for each column (lag pair) */ + xlag = GET(lp[2*j]); + ylag = GET(lp[2*j+1]); + if (N==1) { + r=&rp[j*m]; + x=&xp[xlag]; + y=&yp[ylag]; + bump=r+m; + while(r<bump) { /* for each row */ + d = GET(*x)-GET(*y); + SET(*r) = d*d; + r++; x++; y++; + } + } else { + for (k=0; k<m; k++) { /* for each row */ + z=0; + x=&xp[xlag+k*N]; + y=&yp[ylag+k*N]; + bump=x+N; + while(x<bump) { + d = GET(*x)-GET(*y); + z+=d*d; + x++; y++; + } + SET(rp[j*m+k]) = z; + } + } + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *yp, *rp, *lp, *up, *np; + int nx, mx, ny, my, nr, mr, nl, ml, xmax, ymax, j, N, xlag, ylag; + + /* Check for proper number of arguments */ + if (nrhs < 4 || nrhs > 5) { + mexErrMsgTxt("RDIFF_INPLACE takes 4 or 5 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: X should be doubles"); + } + if (!mxIsNumeric(Y_IN) || mxIsComplex(Y_IN) || + mxIsSparse(Y_IN) || !mxIsDouble(Y_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: Y should be doubles"); + } + if (!mxIsNumeric(R_IN) || mxIsComplex(R_IN) || + mxIsSparse(R_IN) || !mxIsDouble(R_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: R should be doubles"); + } + if (nrhs==5) { + if (!mxIsNumeric(N_IN) || mxIsComplex(N_IN) || + mxIsSparse(N_IN) || !mxIsDouble(N_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: N should double"); + } + np = mxGetPr(N_IN); + N=*np; + } else { + N=1; + } + mx=mxGetM(X_IN); /* rows */ + nx=mxGetN(X_IN); /* columns */ + my=mxGetM(Y_IN); + ny=mxGetN(Y_IN); + mr=mxGetM(R_IN); + nr=mxGetN(R_IN); + ml=mxGetM(L_IN); + nl=mxGetN(L_IN); + + if (nx>1 || ny>1) { + mexErrMsgTxt("RDIFF_INPLACE: X and Y should be column vectors"); + } + if (nr !=nl) { + mexErrMsgTxt("RDIFF_INPLACE: result and lag matrices should have same number of columns"); + } + if ( ml != 2) { + mexErrMsgTxt("RDIFF_INPLACE: lags should be a two row matrix"); + } + if (!nx || !mx || !ny || !my || !nr || !mr || !nl || !ml) { + mexErrMsgTxt("RDIFF_INPLACE: one of the input matrices is empty"); + } + + xp = mxGetPr(X_IN); + yp = mxGetPr(Y_IN); + rp = mxGetPr(R_IN); + lp = mxGetPr(L_IN); + + checkin_matrix((mxArray *) L_IN); + + /* find maximum lag for x and y and check that no lag is negative */ + xmax=GET(lp[0]); + ymax=GET(lp[1]); + for (j=1;j<nl;j++) { + xlag = GET(lp[2*j]); + ylag = GET(lp[2*j+1]); + if ( xlag<0 || ylag<0 ) { + mexErrMsgTxt("RDIFF_INPLACE: LAGS should be non-negative"); + } + if (xlag>xmax ) xmax=xlag; + if (ylag>ymax ) ymax=ylag; + } + + if ( (N*mr+ymax)> my) { + mexPrintf("mr*N: %d, ymax: %d, my: %d\n", mr*N, ymax, my); + mexErrMsgTxt("RDIFF_INPLACE: Y data too short"); + } + if ( (N*mr+xmax)> mx) { + mexPrintf("mr*N: %d, xmax: %d, mx: %d\n", mr*N, xmax, mx); + mexErrMsgTxt("RDIFF_INPLACE: X data too short"); + } + + checkin_matrix((mxArray *) X_IN); + if (yp!=xp) checkin_matrix((mxArray *) Y_IN); + checkin_matrix((mxArray *) R_IN); + + /* Do the actual computations in a subroutine */ + rdiff_inplace(xp,yp,rp,lp,N,mr,nr); + + checkout_matrix((mxArray *) X_IN); + if (yp!=xp) checkout_matrix((mxArray *) Y_IN); + checkout_matrix((mxArray *) R_IN); + checkout_matrix((mxArray *) L_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rdiff_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,19 @@ +% rdiff_inplace(x,y,d,lags,n) - in place running cross-difference function +% +% x: column vector +% y: column vector +% r: result matrix (time X lag) +% lags: 2-column matrix of lags +% n: (samples) frame-rate & window size (default=1) +% +% Vectors x and y are each delayed by amounts specified by lags and subtracted +% sample to sample. The difference is squared and added over a time window +% of n samples. This processing is repeated every n samples, as many times +% as there are columns in r. +% +% A positive lag applied to x causes it to be delayed with respect to y. +% A positive lag applied to y causes it to be delayed with respect to x. +% The number of rows of r must match that of lags. +% +% Mex function. +% Beware: input arguments are assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rdiff_inplace.mexw64.manifest Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +<?xml version='1.0' encoding='UTF-8' standalone='yes'?> +<assembly xmlns='urn:schemas-microsoft-com:asm.v1' manifestVersion='1.0'> + <trustInfo xmlns="urn:schemas-microsoft-com:asm.v3"> + <security> + <requestedPrivileges> + <requestedExecutionLevel level='asInvoker' uiAccess='false' /> + </requestedPrivileges> + </security> + </trustInfo> +</assembly>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rdiff_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rsmooth.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,153 @@ +/* +* rsmooth.c +* smooth columnwise by convolution by a running square window +* multiple passes allowed (--> triangular, gaussian, etc.) +* +* Alain de Cheveigné, CNRS Ircam, Dec 2001. +* (c) CNRS 2001. +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ + +#define X_IN prhs[0] +#define SMOOTH_IN prhs[1] +#define NPASSES_IN prhs[2] +#define CLIP_IN prhs[3] + +/* Output Arguments */ + +#define Y_OUT plhs[0] + +static void rsmooth( + double *xp, + int smooth, + int npasses, + double *yp, + int m, + int n, + int clip + ) +{ + int h,i,j,k,mm; + double *work1, *work2, *a, *b, *tmp, sum; + + mm = m + smooth + npasses*(smooth-1); /* nrows of work buffer */ + work1 = (double *) mxCalloc(mm*n,sizeof(double)); + work2 = (double *) mxCalloc(mm*n,sizeof(double)); + + CHECKIN(work1, (char *) work1 + sizeof(double)*mm*n); + CHECKIN(work2, (char *) work2 + sizeof(double)*mm*n); + + a = work1; + b = work2; + + /* transfer data to buffer, preceded with leading zeros */ + for (k=0; k<n; k++) { + for (j=0;j<m;j++) { + SET(a[ j+smooth+mm*k ]) = GET(xp[ j+m*k ]); + } + } + + /* smooth repeatedly */ + for (h=0;h<npasses;h++) { + for (k=0; k<n; k++) { /* for each column */ + sum=0; + for (j=smooth;j<mm;j++) { /* for each row (inner loop) */ + sum += - GET(a[j-smooth+mm*k]) + GET(a[j+mm*k]); + SET(b[j+mm*k]) = sum/smooth; + } + } + tmp=a; a=b; b=tmp; /* swap for next round */ + } + + /* transfer to output matrix */ + if (clip) { + for (j=0;j<m;j++) { + for (k=0; k<n; k++) { + SET(yp[j+m*k]) = GET(a[ j+(int)(npasses*(smooth-1)/2)+smooth + mm*k]); + } + } + } else { + for (j=0;j<m+npasses*(smooth-1);j++) { + for (k=0; k<n; k++) { + SET(yp[j+(m+npasses*(smooth-1))*k]) = GET(a[j + smooth + mm*k]); + } + } + } + CHECKOUT(work1); + CHECKOUT(work2); + +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *yp; + int m,n,npasses,smooth,clip; + + /* Check for proper number of arguments */ + if (nrhs > 4 || nrhs < 2) { + mexErrMsgTxt("RSMOOTH takes at least 2 and most 4 input arguments"); + } else if (nlhs > 1) { + mexErrMsgTxt("RSMOOTH allows at most 1 output argument"); + } + + /* Check type of input */ + m = mxGetM(X_IN); /* rows */ + n = mxGetN(X_IN); /* columns */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) || m<=1) { + mexErrMsgTxt("RSMOOTH: X should be a column vector or matrix of doubles"); + } + + if (mxGetM(SMOOTH_IN)*mxGetN(SMOOTH_IN) > 1 || mxGetPr(SMOOTH_IN)[0] <1) { + mexErrMsgTxt("RSMOOTH: SMOOTH should be a positive scalar"); + } + smooth = mxGetPr(SMOOTH_IN)[0]; + + if (nrhs > 2) { + npasses = mxGetPr(NPASSES_IN)[0]; + if (mxGetM(NPASSES_IN)*mxGetN(NPASSES_IN) > 1 || npasses <1) { + mexErrMsgTxt("RSMOOTH: NPASSES should be a positive scalar"); + } + } else { + npasses=1; + } + + if (nrhs > 3) { + clip=1; + } else { + clip=0; + } + + /* Create matrix to return */ + if (clip) { + Y_OUT = mxCreateDoubleMatrix(m, n, mxREAL); + } else { + Y_OUT = mxCreateDoubleMatrix(m+npasses*(smooth-1), n, mxREAL); + } + + /* Assign pointers to the various parameters */ + xp = mxGetPr(X_IN); + yp = mxGetPr(Y_OUT); + + checkin_matrix((mxArray *) X_IN); + checkin_matrix(Y_OUT); + + /* Do the actual computations in a subroutine */ + rsmooth(xp,smooth,npasses,yp, m,n,clip); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix(Y_OUT); + return; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rsmooth.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,18 @@ +%y=rsmooth(x,smooth,npasses,trim) - smooth by running convolution +% +% X: input matrix +% SMOOTH: samples - size of square smoothing window +% NPASSES: number of smoothing passes (default=1) +% TRIM: if true, clip Y to same size as X +% +% Y: output matrix +% +% RSMOOTH smooths each column of matrix X by convolution with a square window +% followed by division by the window size. +% Multiple passes allow smoothing with a triangular window (npasses=2), or +% window shapes that approach a gaussian (npasses large). Convolution is +% implemented as a running sum for speed. +% +% Y has NPASSES*(SMOOTH-1) more rows than X unless TRIM is set. +% +% mex function
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rsmooth.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/private/src/rsum_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rsum_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,134 @@ +/* +* rsum_inplace.c +* in-place running sum +* +* Alain de Cheveigné, CNRS/Ircam Jun 2002 +* (c) 2002 CNRS +*/ + +/* Replaces each sample of the input matrix by the sum of itself and its N-1 +* neighbors. The operation is done in-place on the input argument. The last +* N-1 samples of each row are invalid. */ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] +#define N_IN prhs[1] + +/* Output Arguments */ + +static void rsmth( + double *xp, /* matrix to sum */ + double N, /* window size */ + int m, /* rows */ + int n /* columns */ + ) +{ + int j,k, Ni; + double Nr, tmp, sum, *x, *xx, *bump; + + Ni = (int) floor(N); + Nr = N - (double) Ni; + + + if (Nr == 0.0) { + /* N is integer: simple summation */ + for (j=0;j<n;j++) { /* for each column */ + /* prefill running sum */ + sum=0.0; + x=&xp[j*m]; + bump=x+Ni-1; + while (x<bump) { /* for each row */ + sum += GET(*x); + x++; + } + /* N */ + x=&xp[j*m]; + bump=x+m-(Ni-1); + xx=x+Ni-1; + while(x<bump) { + tmp = GET(*x); + sum += GET(*xx); + SET(*x) = sum; + sum -= tmp; + x++; xx++; + } + } + } else { + /* N is fractionary: interpolate */ + for (j=0;j<n;j++) { /* for each column */ + /* prefill running sum */ + sum=0.0; + x=&xp[j*m]; + bump=x+Ni-1; + while (x<bump) { /* for each row */ + sum += GET(*x); + x++; + } + /* N */ + x=&xp[j*m]; + bump=x+m-(Ni-1); + xx=x+Ni-1; + while(x<bump) { + tmp = GET(*x); + sum += GET(*xx); + SET(*x) = sum + Nr * GET(*xx); + sum -= tmp; + x++; xx++; + } + } + } + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *Np, N; + int n, m; + + /* Check for proper number of arguments */ + if (nrhs != 2) { + mexErrMsgTxt("RSUM_INPLACE takes 2 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("RSUM_INPLACE: X should be doubles"); + } + if (!mxIsNumeric(N_IN) || mxIsComplex(N_IN) || + mxIsSparse(N_IN) || !mxIsDouble(N_IN) ) { + mexErrMsgTxt("RSUM_INPLACE: N should be doubles"); + } + m=mxGetM(X_IN); /* rows */ + n=mxGetN(X_IN); /* columns */ + + if (mxGetM(N_IN)*mxGetN(N_IN) != 1) { + mexErrMsgTxt("RSUM_INPLACE: N should be a scalar"); + } + + xp = mxGetPr(X_IN); + Np = mxGetPr(N_IN); + N = *Np; + + if (m < ceil(N)) { + mexErrMsgTxt("RSUM_INPLACE: X should have at least N rows"); + } + + checkin_matrix((mxArray *) X_IN); + + /* Do the actual computations in a subroutine */ + rsmth(xp, N, m, n); + + checkout_matrix((mxArray *) X_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rsum_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +% rsum_inplace(x,N) - in place running sum +% +% x: matrix to sum +% N: samples - size of running sum window +% +% Matrix x is summed row-wise. The last (N-1) samples are not valid. +% N may be fractionary (handled by linear interpolation). +% +% Mex function. +% Beware: input argument is assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/rsum_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/sf_cleanup.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,29 @@ +function i = sf_cleanup(i) +% i=sf_cleanup(i) - cleanup after use (close i.fd if open) + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin ; error('sf_info: no input arguement'); end +if ~isa(i, 'struct') error('sf_info: expected struct'); end + +% close file if open +if isfield(i, 'fd') & fopen(i.fd) + fclose(i.fd); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/private/src/sf_format.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,213 @@ +function i = sf_format(i) +% i=sf_format(i) - try to guess file format from +% magic words or file suffix + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin; error('no argument'); end +if ~isa(i, 'struct') + i.fname=i; + i=sf_format(i); + return +end +if ~isfield(i, 'fname'); error('no fname field'); end + +if isa(i.fname, 'double') + i.format = 'matrix'; return +end + +if isa(i.fname, 'char') + disp(['guessing format of ', i.fname, '...']); +else + disp('guessing format...'); +end + +% Well-known file extensions +a = findstr('.', i.fname); +if a + suff = i.fname( a(size(a,2))+1 : size(i.fname,2)); + suff = lower(suff); + switch suff + case {'short', 'long', 'float', 'double', 'ascii', 'mat'} + i.format = suff; + return +% the following are commented out to exercise header magic +% case {'aiff'} +% i.format = 'AIFF'; +% return +% case {'aifc'} +% i.format = 'AIFC'; +% return +% case 'wav' +% i.format = 'WAV'; +% return +% case 'au' +% i.format = 'AU'; +% case 'sdif' +% i.format = 'sdif'; + % others ? + end +end + +% close file if open +if exist('i.fd') & fopen(i.fd) + fclose(i.fd); +end + + +% open little-endian, check for .wav magic +[i.fd, msg] = fopen(i.fname, 'r', 'l'); +if i.fd == -1 + error(['could not open: >', i.fname, '<']); + if ~isempty(msg) + error(msg); + end +end + +% file size in bytes +if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; +i.nbytes = ftell(i.fd); +if i.nbytes == 0 & strcmp(computer, 'MAC2') + disp('file is empty: perhaps its a macintosh SND resource?'); + i.format = 'MACSND'; % macintosh sound resource (maybe...) + return +end +if i.nbytes < 4; error('file is less than 4 bytes long'); end; +fseek(i.fd, 0, -1); +magic = char(fread(i.fd, 4, 'uchar'))'; +if strcmp(magic, 'RIFF') & i.nbytes >=12 + dummy = fread(i.fd, 1, 'ulong'); % chunk size + magic = char(fread(i.fd, 4, 'uchar'))'; + if strcmp(magic, 'WAVE'); + i.format='WAV'; + return + end +end + +% reopen for standard binary, look for magic +fclose(i.fd); +i.fd = fopen(i.fname, 'r'); +if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; +i.nbytes = ftell(i.fd); +fseek(i.fd, 0, -1); +% magic strings? +magic = char(fread(i.fd, 4, 'uchar'))'; +switch magic +case '.snd' + i.format = 'AU'; + return; +case 'FORM' + fseek(i.fd, 8, -1); + i.format = char(fread(i.fd, 4, 'uchar'))'; + if ~strcmp(i.format, 'AIFF') & ~strcmp(char(i.format),'AIFC') + error (['expected AIFF or AIFC, found ' i.format]); + end; + return; +case 'NIST' + fseek(i.fd, 0, -1); + line = fscanf(i.fd, '%s' , 1); + if ~strcmp(line, 'NIST_1A'); + error(['expected NIST_1A, found ', magic]); + end; + disp(i.format); + return; +case 'SDIF' + i.format = 'SDIF'; + return; +end +% wff magic number? +wff_magic = 195894762; +fseek(i.fd, 0, -1); +magic = fread(i.fd, 1, 'long'); +if magic == wff_magic + i.format = 'wff'; % Nottingham IHR wff format + return +end +% ESPS? +esps_min_size = 333; esps_magic = 27162; +if i.nbytes > esps_min_size + fseek(i.fd, 16, -1); + magic = fread(i.fd, 1, 'uint'); + if magic == esps_magic + i.format = 'ESPS'; + return + end +end + + +% Magic didn't succeed, try to determine type based on file name +a = findstr('.', i.fname); +if a + suff = i.fname( a(size(a,2))+1 : size(i.fname,2)); + switch suff + case 's' + i.format = 'short'; + return + case 'f' + i.format = 'float'; + return + % others ? + end +end + + +% Ascii formats +fclose(i.fd); +i.fd = fopen(i.fname, 'rt'); +line = fscanf(i.fd, '%s' , 1); +if strcmp(line, '|WAVE'); + i.format = 'iwave'; % Nottingham IHR |WAVE format + return; +end + +% plain ascii ? +limit = 2048; % test only this many bytes +hunk = char(fread(i.fd, min(limit, i.nbytes))); +chars = unique(hunk); +if isempty(setdiff(chars, [9, 10, 13, ... + double(' abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890+-_.,:;')])) + % ascii + if ~isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.,')])) + % not just numerical: try removing first line + fseek(i.fd, 0, -1); + line = fgets(i.fd); + if -1 == line; return ; end % line too long + i.bytes_to_data = ftell(i.fd); + hunk = char(fread(i.fd, min(limit, i.nbytes-ftell(i.fd)))); + chars = unique(hunk); + warning('skipping first line of what seems like an ascii file'); + end + if isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.')])) + i.format = 'ascii'; + return + end + if isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.,')])) + w.format = 'csv'; + return + end +end + +% Default +if i.nbytes == 2 * round(i.nbytes/2); % even number of bytes + i.format = 'short'; +else + i.format = 'uchar'; +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/testyin.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,6 @@ +% test YIN + +%yin 'female04_16K.wav' +clear; +clc; +YINPitchDetection( '1.wav'); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +function [prd, fileinfo] = yin(x,p) % YIN - Fundamental frequency (F0) of file or vector % % YIN('FILE') estimates and plots the F0 of signal in FILE. % F0 is plotted in octaves re: 440 Hz, together with residual % (aperiodicity) and power measures. A 'best' estimate is printed. % YIN(X,SR) estimates and plots the F0 of array X assuming sampling rate SR (Hz). % % R=YIN(X) produces no plot but returns result in structure R: % R.f0: fundamental frequency in octaves re: 440 Hz % R.ap: aperiodicity measure (ratio of aperiodic to total power) % R.pwr: period-smoothed instantaneous power % (R also lists the parameters used by YIN) % % YIN(NAME,P) uses parameters stored in P: % P.minf0: Hz - minimum expected F0 (default: 30 Hz) % P.maxf0: Hz - maximum expected F0 (default: SR/(4*dsratio)) % P.thresh: threshold (default: 0.1) % P.relfag: if ~0, thresh is relative to min of difference function (default: 1) % P.hop: samples - interval between estimates (default: 32) % P.range: samples - range of samples ([start stop]) to process % P.bufsize: samples - size of computation buffer (default: 10000) % P.sr: Hz - sampling rate (usually taken from file header) % P.wsize: samples - integration window size (defaut: SR/minf0) % P.lpf: Hz - intial low-pass filtering (default: SR/4) % P.shift 0: shift symmetric, 1: shift right, -1: shift left (default: 0) % % See 'yin.html' for more info. % Version 28 July 2003. % Alain de Cheveign? CNRS/Ircam, 2002. % Copyright (c) 2002 Centre National de la Recherche Scientifique. % % Permission to use, copy, modify, and distribute this software without % fee is hereby granted FOR RESEARCH PURPOSES only, provided that this % copyright notice appears in all copies and in all supporting % documentation, and that the software is not redistributed for any % fee (except for a nominal shipping charge). % % For any other uses of this software, in original or modified form, % including but not limited to consulting, production or distribution % in whole or in part, specific prior permission must be obtained from CNRS. % Algorithms implemented by this software may be claimed by patents owned % by CNRS, France Telecom, Ircam or others. % % The CNRS makes no representations about the suitability of this % software for any purpose. It is provided "as is" without express % or implied warranty. % Hidden parameters are integration window size (set equal to sr/minf0), search % range for 'best neighboring estimate' (set equal to +/- sr/(2*minf0)), maximum % expected width of period dip (stop/start ratio == 1.85), margin for "beam search" % of final estimate (+1.8/-0.6 times the initial estimate). % default parameter values ([]: to be determined) minf0 = 30; % Hz - minimum frequency maxf0 = []; % Hz - maximum frequency wsize = []; % s - integration window size lpf = []; % Hz - lowpass prefiltering cutoff thresh = 0.1; % difference function threshold relflag = 1; % if true threshold is relative to global min of difference function bufsize=10000; % computation buffer size hop = 32; % samples - interval between estimates range=[]; % range of file samples to process sr=[]; % sampling rate shift=0; % flag to control the temporal shift of analysis windows (left/sym/right) plotthreshold=0.2; % aperiodicity above which plot is green or yellow % if 2~=exist('allread') % error('sf routines missing: put them in your path & try again'); % end % handle parameters if nargin<1; help yin; return; end if nargin<2; p=[]; end fileinfo=sf_info(x); if ~isempty(fileinfo.sr) p.sr=fileinfo.sr; end % get sr from file if fileinfo.nchans > 1 disp(['warning: using column 1 of ', num2str(fileinfo.nchans), '-column data']); end if isa(p, 'double') p.sr=p; end if ~isfield(p, 'sr'); p.sr=sr; end if isempty(p.sr); error('YIN2: must specify SR'); end if ~isfield(p, 'range') | isempty(p.range); p.range=[1 fileinfo.nsamples]; end if ~isfield(p, 'minf0'); p.minf0=minf0; end if ~isfield(p, 'thresh'); p.thresh=thresh; end if ~isfield(p, 'relflag'); p.relflag=relflag; end if ~isfield(p, 'bufsize'); p.bufsize=bufsize; end if ~isfield(p, 'hop'); p.hop=hop; end if ~isfield(p, 'maxf0'); p.maxf0=floor(p.sr/4); end % default if ~isfield(p, 'wsize'); p.wsize=ceil(p.sr/p.minf0); end % default if ~isfield(p, 'lpf'); p.lpf=p.sr/4; end % default if mod(p.hop,1); error('hop should be integer'); end if ~isfield(p, 'shift'); p.shift=shift; end % default if ~isfield(p, 'plotthreshold'); p.plotthreshold=plotthreshold; end % default % estimate period r=yink(p,fileinfo); prd=r.r1; % period in samples ap0=r.r2; % gross aperiodicity measure ap= r.r3; % fine aperiodicity measure pwr=r.r4; % period-smoothed instantaneous power f0 = log2(p.sr ./ prd) - log2(440); % convert to octaves re: 440 Hz %figure(); %plot(prd); %dlmwrite('pitch.txt',prd,'delimiter','\n'); %display(length(prd)); % load estimates and major parameters in result structure clear r; r.f0 = f0; r.ap0 = ap0; r.ap = ap; r.pwr = pwr; r.sr = p.sr; r.range=p.range; r.minf0 = p.minf0; r.maxf0 = p.maxf0; r.thresh=p.thresh; r.relflag=p.relflag; r.hop = p.hop; r.bufsize = p.bufsize; r.wsize = p.wsize; r.lpf = p.lpf; r.shift = p.shift; r.plotthreshold=p.plotthreshold; % plot estimates (if nargout == 0) if nargout<1 if isnan(f0) display('no estimates: signal too short or buffersize too small'); return; end % choose sample to report as "the" f0 of the entire signal [mn, idx] = min(ap0); best=f0(idx); disp(['best: ', num2str(2^best*440), 'Hz (', note(best),... ') at ', num2str(idx/(p.sr/p.hop)), 's (aperiodic/total power: ', num2str(mn), ')']); % plot f0 in 3 colors according to periodicity good = f0; good(find(ap0>p.plotthreshold*2)) = nan; best = f0; best(find(ap0>p.plotthreshold)) = nan; subplot(211); fsr=p.sr/p.hop; nframes=size(prd,2); if nframes <2; error('F0 track is too short to plot'); end plot((1:nframes)/fsr, f0, 'y', (1:nframes)/fsr, good, 'g', (1:nframes)/fsr, best, 'b'); lo = max(min(f0),min(good)); hi=min(max(f0),max(good)); set(gca, 'ylim', [lo-0.5; hi+0.5]); set(gca, 'xlim', [1,nframes]/fsr); set(get(gca,'ylabel'), 'string', 'Oct. re: 440 Hz'); set(get(gca,'xlabel'), 'string', 's'); % plot periodicity ap0=max(0,ap0); ap=max(0,ap); ap1=ap; ap(find(ap0>plotthreshold)) = nan; subplot(413); plot((1:nframes), ap0.^0.5, 'b'); %subplot(413); plot((1:nframes), ap0.^0.5, 'c', (1:nframes), ap1.^0.5, 'y', (1:nframes), ap.^0.5, 'b'); set(gca, 'ylim', [0 1]); set(gca, 'xlim', [1, nframes]); set(get(gca,'ylabel'), 'string', 'sqrt(apty)'); % plot power subplot(414); plot((1:nframes), sqrt(pwr), 'b'); set(gca, 'xlim', [1, nframes]); set(get(gca,'ylabel'), 'string', 'sqrt(power)'); set(get(gca,'xlabel'), 'string', 'frames'); if isa(x, 'double') set(gcf, 'Name', 'workspace matrix'); else set (gcf, 'Name', x); end end % convert octave re 440 to note: function s=note(o) n=round(12*o); cents = 100*(12*o-n); oct=floor((n-3)/12)+5; chroma=mod(n,12); chromalist = {'A'; 'A#'; 'B'; 'C'; 'C#'; 'D'; 'D#'; 'E'; 'F'; 'F#';... 'G'; 'G#'}; cents = sprintf('%+.0f',cents); s=[char(chromalist(chroma+1)),num2str(oct),' ',num2str(cents), ' cents']; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/junk/yin2.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,103 @@ +function r=yin2(p,fileinfo) +% YIN2 - fundamental frequency estimator +% new version (feb 2003) +% +% + + +% process signal a chunk at a time +idx=0; +totalhops=round(fileinfo.nsamples / p.hop); +r1=nan*zeros(1,totalhops);r2=nan*zeros(1,totalhops); +r3=nan*zeros(1,totalhops);r4=nan*zeros(1,totalhops); +idx2=0+round(p.wsize/2/p.hop); +while (1) + xx=sf_wave(fileinfo, [idx+1, idx+p.bufsize], []); + xx=xx(:,1); % first channel if multichannel + [prd,ap0,ap,pwr]=yin_helper(xx,p); + n=size(prd ,2); + if (~n) break; end; + idx=idx+n*p.hop; + + r1(idx2+1:idx2+n)= prd; + r2(idx2+1:idx2+n)= ap0; + r3(idx2+1:idx2+n)= ap; + r4(idx2+1:idx2+n)= pwr; + idx2=idx2+n; +end +size(r1) +r.r1=r1; % period estimate +r.r2=r2; % gross aperiodicity measure +r.r3=r3; % fine aperiodicity measure +r.r4=r4; % power +sf_cleanup(fileinfo); +% end of program + + + +% Estimate F0 of a chunk of signal +function [prd,ap0,ap,pwr]=yin_helper(x,p,dd) +[m,n]=size(x); +smooth=ceil(p.sr/p.lpf); +x=rsmooth(x,smooth); % light low-pass smoothing +x=x(smooth+1:end); +maxlag = ceil(p.sr/p.minf0); +minlag = floor(p.sr/p.maxf0); + +hops=floor((m-maxlag-p.wsize)/p.hop); +prd=zeros(1,hops); +ap0=zeros(1,hops); +ap=zeros(1,hops); +pwr=zeros(1,hops); +if hops<1; return; end + +% difference function matrix +dd=zeros(ceil((m-maxlag-2)/p.hop),maxlag+2); % +2 to improve interp near maxlag +lags=[zeros(1,maxlag+2); 1:maxlag+2]; +rdiff_inplace(x,x,dd,lags,p.hop); +rsum_inplace(dd,round(p.wsize/p.hop)); +dd=dd'; + +% parabolic interpolation near min, then cumulative mean normalization +[dd,ddx]=minparabolic(dd); +cumnorm_inplace(dd); + +for j=0:hops-1 + idx=j*p.hop; + d=dd(:,j+1); + dx=ddx(:,j+1); + + % estimate period + pd=dftoperiod(d,[minlag,maxlag],p.thresh*2); + + % gross aperiodicity is value of cumnormed df at minimum: + ap0(j+1)=d(pd)/2; + + % fine tune period based on parabolic interpolation + pd=pd+dx(pd+1)+1; + + % power estimates + k=(1:ceil(pd))'; + x1=x(k+idx); + x2=k+idx+pd-1; + interp_inplace(x,x2); + x3=x2-x1; + x4=x2+x1; + x1=x1.^2; rsum_inplace(x1,pd); + x3=x3.^2; rsum_inplace(x3,pd); + x4=x4.^2; rsum_inplace(x4,pd); + + x1=x1(1)/pd; + x2=x2(1)/pd; + x3=x3(1)/pd; + x4=x4(1)/pd; + % total power + pwr(j+1)=x1; + + % fine aperiodicity + ap(j+1)=(eps+x3)/(eps+(x3+x4)); % accurate, only for valid min + + prd(j+1)=pd; +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/cumnorm_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,9 @@ +% cumnorm_inplace(x) - cumulative mean-normalization +% +% x: matrix to normalize +% +% Each column of x is normalized by dividing each sample by the mean of +% samples with indices smaller or equal to it. +% +% Mex function. +% Warning: the input argument is assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/dftoperiod copy.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,14 @@ +% prd=dftoperiod(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) +% is the period estimate. +% +% Mex function.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/dftoperiod.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,13 @@ +% prd=dftoperiod(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) is the period estimate. +% +% Mex function.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/dftoperiod2.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,17 @@ +% prd=dftoperiod2(d,b,t) - estimate period from difference function +% +% prd: row matrix of period estimates +% +% d: column vector or matrix of difference functions +% b: bounds matrix ([lo, hi]) +% t: threshold +% +% Difference functions are supposed to be cumulative-mean-normalized. +% For each column of d, search for the first minimum between lo and hi that +% falls below threshold. The index of this minimum (re 0) +% is the period estimate. +% +% This version differs from dftoperiod in that the threshold is +% incremented by the global minimum of the difference function. +% +% Mex function.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/minparabolic.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,9 @@ +% [y,idx] = minparabolic(x) - locate minima using parabolic interpolation +% +% x: input vector +% y: same as x, with local minima replaced by interpolated minima +% idx: position of interpolated minima (fractionary) +% +% At each local minimum of x, that sample and its neighbors are used to +% determine a parabola. The minimum of the parabola is used as the new +% minimum value, the abscissa of the minimum is used as its position.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/rdiff_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,19 @@ +% rdiff_inplace(x,y,d,lags,n) - in place running cross-difference function +% +% x: column vector +% y: column vector +% r: result matrix (time X lag) +% lags: 2-column matrix of lags +% n: (samples) frame-rate & window size (default=1) +% +% Vectors x and y are each delayed by amounts specified by lags and subtracted +% sample to sample. The difference is squared and added over a time window +% of n samples. This processing is repeated every n samples, as many times +% as there are columns in r. +% +% A positive lag applied to x causes it to be delayed with respect to y. +% A positive lag applied to y causes it to be delayed with respect to x. +% The number of rows of r must match that of lags. +% +% Mex function. +% Beware: input arguments are assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/rsmooth.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,18 @@ +%y=rsmooth(x,smooth,npasses,trim) - smooth by running convolution +% +% X: input matrix +% SMOOTH: samples - size of square smoothing window +% NPASSES: number of smoothing passes (default=1) +% TRIM: if true, clip Y to same size as X +% +% Y: output matrix +% +% RSMOOTH smooths each column of matrix X by convolution with a square window +% followed by division by the window size. +% Multiple passes allow smoothing with a triangular window (npasses=2), or +% window shapes that approach a gaussian (npasses large). Convolution is +% implemented as a running sum for speed. +% +% Y has NPASSES*(SMOOTH-1) more rows than X unless TRIM is set. +% +% mex function
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/rsum_inplace.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +% rsum_inplace(x,N) - in place running sum +% +% x: matrix to sum +% N: samples - size of running sum window +% +% Matrix x is summed row-wise. The last (N-1) samples are not valid. +% N may be fractionary (handled by linear interpolation). +% +% Mex function. +% Beware: input argument is assigned to. This is not matlab-kosher!
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/sf_cleanup.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,29 @@ +function i = sf_cleanup(i) +% i=sf_cleanup(i) - cleanup after use (close i.fd if open) + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin ; error('sf_info: no input arguement'); end +if ~isa(i, 'struct') error('sf_info: expected struct'); end + +% close file if open +if isfield(i, 'fd') & fopen(i.fd) + fclose(i.fd); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/sf_format.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,213 @@ +function i = sf_format(i) +% i=sf_format(i) - try to guess file format from +% magic words or file suffix + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin; error('no argument'); end +if ~isa(i, 'struct') + i.fname=i; + i=sf_format(i); + return +end +if ~isfield(i, 'fname'); error('no fname field'); end + +if isa(i.fname, 'double') + i.format = 'matrix'; return +end + +if isa(i.fname, 'char') + disp(['guessing format of ', i.fname, '...']); +else + disp('guessing format...'); +end + +% Well-known file extensions +a = findstr('.', i.fname); +if a + suff = i.fname( a(size(a,2))+1 : size(i.fname,2)); + suff = lower(suff); + switch suff + case {'short', 'long', 'float', 'double', 'ascii', 'mat'} + i.format = suff; + return +% the following are commented out to exercise header magic +% case {'aiff'} +% i.format = 'AIFF'; +% return +% case {'aifc'} +% i.format = 'AIFC'; +% return +% case 'wav' +% i.format = 'WAV'; +% return +% case 'au' +% i.format = 'AU'; +% case 'sdif' +% i.format = 'sdif'; + % others ? + end +end + +% close file if open +if exist('i.fd') & fopen(i.fd) + fclose(i.fd); +end + + +% open little-endian, check for .wav magic +[i.fd, msg] = fopen(i.fname, 'r', 'l'); +if i.fd == -1 + error(['could not open: >', i.fname, '<']); + if ~isempty(msg) + error(msg); + end +end + +% file size in bytes +if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; +i.nbytes = ftell(i.fd); +if i.nbytes == 0 & strcmp(computer, 'MAC2') + disp('file is empty: perhaps its a macintosh SND resource?'); + i.format = 'MACSND'; % macintosh sound resource (maybe...) + return +end +if i.nbytes < 4; error('file is less than 4 bytes long'); end; +fseek(i.fd, 0, -1); +magic = char(fread(i.fd, 4, 'uchar'))'; +if strcmp(magic, 'RIFF') & i.nbytes >=12 + dummy = fread(i.fd, 1, 'ulong'); % chunk size + magic = char(fread(i.fd, 4, 'uchar'))'; + if strcmp(magic, 'WAVE'); + i.format='WAV'; + return + end +end + +% reopen for standard binary, look for magic +fclose(i.fd); +i.fd = fopen(i.fname, 'r'); +if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; +i.nbytes = ftell(i.fd); +fseek(i.fd, 0, -1); +% magic strings? +magic = char(fread(i.fd, 4, 'uchar'))'; +switch magic +case '.snd' + i.format = 'AU'; + return; +case 'FORM' + fseek(i.fd, 8, -1); + i.format = char(fread(i.fd, 4, 'uchar'))'; + if ~strcmp(i.format, 'AIFF') & ~strcmp(char(i.format),'AIFC') + error (['expected AIFF or AIFC, found ' i.format]); + end; + return; +case 'NIST' + fseek(i.fd, 0, -1); + line = fscanf(i.fd, '%s' , 1); + if ~strcmp(line, 'NIST_1A'); + error(['expected NIST_1A, found ', magic]); + end; + disp(i.format); + return; +case 'SDIF' + i.format = 'SDIF'; + return; +end +% wff magic number? +wff_magic = 195894762; +fseek(i.fd, 0, -1); +magic = fread(i.fd, 1, 'long'); +if magic == wff_magic + i.format = 'wff'; % Nottingham IHR wff format + return +end +% ESPS? +esps_min_size = 333; esps_magic = 27162; +if i.nbytes > esps_min_size + fseek(i.fd, 16, -1); + magic = fread(i.fd, 1, 'uint'); + if magic == esps_magic + i.format = 'ESPS'; + return + end +end + + +% Magic didn't succeed, try to determine type based on file name +a = findstr('.', i.fname); +if a + suff = i.fname( a(size(a,2))+1 : size(i.fname,2)); + switch suff + case 's' + i.format = 'short'; + return + case 'f' + i.format = 'float'; + return + % others ? + end +end + + +% Ascii formats +fclose(i.fd); +i.fd = fopen(i.fname, 'rt'); +line = fscanf(i.fd, '%s' , 1); +if strcmp(line, '|WAVE'); + i.format = 'iwave'; % Nottingham IHR |WAVE format + return; +end + +% plain ascii ? +limit = 2048; % test only this many bytes +hunk = char(fread(i.fd, min(limit, i.nbytes))); +chars = unique(hunk); +if isempty(setdiff(chars, [9, 10, 13, ... + double(' abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890+-_.,:;')])) + % ascii + if ~isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.,')])) + % not just numerical: try removing first line + fseek(i.fd, 0, -1); + line = fgets(i.fd); + if -1 == line; return ; end % line too long + i.bytes_to_data = ftell(i.fd); + hunk = char(fread(i.fd, min(limit, i.nbytes-ftell(i.fd)))); + chars = unique(hunk); + warning('skipping first line of what seems like an ascii file'); + end + if isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.')])) + i.format = 'ascii'; + return + end + if isempty(setdiff(chars, [9, 10, 13, double(' 1234567890eE+-.,')])) + w.format = 'csv'; + return + end +end + +% Default +if i.nbytes == 2 * round(i.nbytes/2); % even number of bytes + i.format = 'short'; +else + i.format = 'uchar'; +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/sf_info.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,438 @@ +function i = sf_info(i) +% i=sf_info(i) - extract useful info from file + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin ; error('sf_info: no input arguement'); end +if ~isa(i, 'struct') + j.fname=i; + i=j; + i=sf_info(i); + return +end +if ~isfield(i, 'fname'); error('sf_info: no fname field'); end + +% guess format only if unknown +if ~isfield(i, 'format') | isempty(i.format); + i = sf_format(i); + if ~strcmp(i.format, 'matrix') disp(i.format); end +end + +% handle workspace matrices as if they were files +if strcmp(i.format, 'matrix') + [nrows, ncols] = size(i.fname); + i.nchans=ncols; + i.nsamples=nrows; + i.totalsamples = i.nsamples*i.nchans; + i.sr=[]; + return; +end + +% close file if open +if exist('i.fd') & fopen(i.fd) + fclose(i.fd); +end + +% use standard matlab functions for AU and WAV and MACSND +if strcmp(i.format, 'AU') + if isempty(findstr('.', i.fname)) + disp(['WARNING: matlab function AUREAD requires '... + '.au or .snd suffix on file name: ' i.fname]); + end + sz = auread(i.fname, 'size'); + i.nsamples=sz(1); + i.nchans=sz(2); + i.totalsamples = i.nsamples*i.nchans; + [dummy, i.sr, i.samplebits] = auread(i.fname, 1); + return; +end +if strcmp(i.format, 'WAV') + if isempty(findstr('.wav', i.fname)) & isempty(findstr('.WAV', i.fname)) + disp(['WARNING: matlab function WAVREAD requires '... + '.wav suffix on file name: ' i.fname]); + end + sz = wavread(i.fname, 'size'); + i.nsamples=sz(1); + i.nchans=sz(2); + i.totalsamples = i.nsamples*i.nchans; + [dummy, i.sr, i.samplebits] = wavread(i.fname, 1); + return; +end +if strcmp(i.format, 'MACSND') + % must load the data to get info - this is stupid + if ~isempty(findstr(':', i.fname)) + disp(['matlab function READSND cannot handle an '... + 'indirect path: ' i.fname]); + end + if 3==exist('readsnd') + [data, i.sr] = eval('readsnd(i.fname)'); + else + error('cannot read MACSND on this platform'); + end + i.nsamples = size(data,2); + i.nchans = size(data,1); + i.totalsamples = i.nsamples*i.nchans; + return +end + +% reopen file +if strcmp(i.format, 'ascii') ... + | strcmp(i.format, 'csv') | strcmp(i.format, '|WAVE') + [i.fd, msg] = fopen(i.fname, 'rt'); +else + [i.fd, msg] = fopen(i.fname, 'r', 'ieee-be.l64'); +end +if i.fd == -1 + if isempty(msg) + error(['could not open: ', i.fname]); + else + error(msg) + end +end +fd = i.fd; +if ~isfield(i, 'nbytes') + if (-1 == fseek(i.fd, 0, 1)) ; error ('fseek failed'); end; + i.nbytes = ftell(i.fd); +end +fseek(fd, 0, -1); % rewind + + +switch i.format +%%%%%%%%%%%%%%%%%%% AIFF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case {'AIFF','AIFC'} + fseek(fd, 12, 0); % skip container chunk + % skip over spurious chunks + idx=ftell(fd); + while 1 + magic=char(fread(fd,4,'uchar'))'; + if strcmp(magic,'COMM'); break; end; + idx = idx+1; + status = fseek(fd,idx,-1); + if status == -1 + error('expected COMM magic word, found eof'); + end; + end; + + %ckSize=fread(fd,1,'int32'); + %status = fseek(fd, ckSize, 0); % skip to end of chunk + %if status == -1 + % error('unexpected eof'); + %end + + %while (1) + % magic = char(fread(fd, 4, 'char'))'; + % if ~strcmp(magic, 'SSND') + % fseek(fd, -4, 0); % skip back + % break; + % end + % ckSize = fread(fd, 1, 'long'); + % fseek(fd, ckSize, 0); % skip to end of sound chunk + %end + %magic = char(fread (fd, 4, 'char'))'; + %if ~strcmp(magic, 'COMM'); error(['expected COMM, found ', magic]) ; end + commsz = fread(fd, 1, 'int32'); + i.nchans = fread(fd, 1, 'int16'); + i.nsamples = fread(fd, 1, 'uint32'); + i.totalsamples = i.nsamples*i.nchans; + i.samplebits = fread(fd, 1, 'int16'); + switch i.samplebits + case 16 + i.sample_bytes = 2; + i.sample_type = 'int16'; + case 32 + i.sample_bytes = 4; + i.sample_type = 'int32'; % or float? + otherwise + error(['unexpected samplebits: ' num2str(i.samplebits) ]); + end + % read sampling rate using Hideki Kawahara's code: + srex1=fread(fd,1,'uint16'); + srex2=fread(fd,1,'uint64'); + if strcmp(char(i.format),'AIFC') + compress=fread(fd,4,'uchar'); + if ~strcmp(char(compress),'NONE') + error('Compression is not supported.'); + end; + fseek(fd, commsz-22, 0); + end; + i.sr = 2^(srex1-16383)*srex2/hex2dec('8000000000000000'); + %fseek(fd, 12, -1); % skip back to end of container chunk + % skip over eventual common chunk + %while(1) + % magic = char(fread(fd, 4, 'char'))'; + % if ~strcmp(magic, 'COMM') + % fseek(fd, -4, 0); % skip back + % break; + % end + % ckSize = fread(fd, 1, 'long'); + % fseek(fd, ckSize, 0); % skip over chunk + %end + magic=char(fread(fd,4,'uchar'))'; + while ~strcmp(char(magic),'SSND') + [ckSize, count]=fread(fd,1,'int32'); + if ~count; + error('expected chunk size field, found eof'); + return + end + status = fseek(fd, ckSize, 0); % skip to end of chunk + if status == -1 + error('expected SSND magic word, found eof'); + return + end; + magic=char(fread(fd,4,'uchar'))'; + end; + + %magic = char(fread(fd, 4, 'char'))'; + %if ~strcmp(magic, 'SSND') + % error (['expected SSND, found' magic]); + %end + fseek(fd, 12, 0); % skip over ckSize, offset and blocksize fields + i.bytes_to_data = ftell(fd); + if i.totalsamples*i.sample_bytes ~= i.nbytes-i.bytes_to_data + disp(['WARNING: header fields sample_bytes: ' ... + num2str(i.sample_bytes)]); + disp (['and sample and channel count: ' num2str(i.nsamples) ... + ', ' num2str(i.nchans)]); + disp(['are inconsistent with offset to data: ' ... + num2str(i.bytes_to_data) ' and file size: ' ... + num2str(i.nbytes)]); + disp (['(' num2str(i.totalsamples*i.sample_bytes) ... + ' ~= ' num2str(i.nbytes-i.bytes_to_data) ')']); + end + + return +%%%%%%%%%%%%%%%%%%% NIST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'NIST' +% fseek(fd, 0, -1); +% line = fscanf(fd, '%s' , 1); +% if ~strcmp(line, 'NIST_1A'); error(['expected NIST_1A, found ', magic]); end + fseek(fd, 8, 0); % skip over magic string + i.bytes_to_data = fscanf(fd, '%d', 1); + while (1) + key = fscanf(fd, '%s', 1); + if strcmp(key, 'end_head'); break; end; + % read third field according to spec in second field (type) + % this may need refining... + type = fscanf(fd, '%s', 1); + if strcmp(type(1:2), '-s') + bytes_to_read = sscanf(type(3:end), '%d', 1); + fseek(fd, 1, 0); % skip blank + value = char(fread(fd, bytes_to_read, 'char'))'; + else + value = fscanf(fd, '%f', 1); + end + i = setfield(i, key, value); + end + % give standard names to useful fields + if isfield(i, 'channel_count'); i.nchans = i.channel_count; end + if isfield(i, 'sample_count'); i.nsamples = i.sample_count; end + i.totalsamples = i.nsamples*i.nchans; + if isfield(i, 'sample_rate'); + i.sr = i.sample_rate; + i.xunits = 's'; + end + if ~isfield(i, 'sample_coding'); i.sample_coding = 'pcm'; end + i.bytes_to_data = 1024; % needs checking + i.sample_bytes=2; + i.sample_type='int16'; + return +%%%%%%%%%%%%%%%%%%% |WAVE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case '|WAV' + line = fscanf(fd, '%s' , 1); % skip first line + i.nsamples = fscanf(fd, '%d', 1); + i.nchans = fscanf(fd, '%d', 1); + i.totalsamples = i.nsamples*i.nchans; + i.bytes_to_data = ftell(fd); + % channel info handled during data read + i.totalsamples = i.nsamples * i.nchans; + return +%%%%%%%%%%%%%%%%%%% WFF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'wff' + fseek(fd, 4, 0); % skip magic number + i.version = fread(fd, 1, 'long'); + i.type_code = fread(fd, 1, 'long'); + i.nchans = fread(fd, 1, 'long'); + i.info.channel_flags = fread(fd, 1, 'long'); + i.bytes_to_data = fread(fd, 1, 'long'); + fseek(fd, 40, 0); + i.gen_prog_name = fread(fd, 32, 'char'); + i.comment = fread(fd, 32, 'char'); + i.sample_bytes = 2; + i.sample_type = 'int16'; + return; +% channel info handled later during channel read +%%%%%%%%%%%%%%%%%%% ESPS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'ESPS' +% Based on Peter Kabal's afsp package. +% This handles at least one kind of ESPS waveform file. Others? + i.machine_code = fread(fd, 1, 'uint'); + i.version_check_code = fread(fd, 1, 'uint'); + i.bytes_to_data = fread(fd, 1, 'uint'); + i.record_size = fread(fd, 1, 'uint'); + fseek(fd, 20, -1); + i.EDR_ESPS_flag = fread(fd, 1, 'uint'); + i.align_pad_size = fread(fd, 1, 'uint'); + fseek(fd, 32, -1); + i.file_type = fread(fd, 1, 'uint16'); + fseek(fd, 40, -1); + i.file_creation_date_time = char(fread(fd, 26, 'char'))'; + i.header_version = char(fread(fd, 8, 'char'))'; + i.program_name = char(fread(fd, 16, 'char'))'; + i.program_version = char(fread(fd, 8, 'char'))'; + i.compile_date = char(fread(fd, 26, 'char'))'; + i.tag = fread(fd, 1, 'uint'); + fseek(fd, 132, -1); + i.ndoubles = fread(fd, 1, 'uint'); + i.nfloats = fread(fd, 1, 'uint'); + i.nlongs = fread(fd, 1, 'uint'); + i.nshorts = fread(fd, 1, 'uint'); + i.nchars = fread(fd, 1, 'uint'); + i.fixed_header_size = fread(fd, 1, 'uint'); + i.var_header_size = fread(fd, 1, 'uint'); + %w.dunno_what = fread(fd, 1, 'uint'); + fseek(fd, 160, -1); + i.user = char(fread(fd, 8, 'char'))'; + % scan the rest of the header to find sampling rate + a = ftell(fd); + bytes_left_in_header = i.bytes_to_data - ftell(fd); + hunk = char(fread(fd, bytes_left_in_header, 'uchar'))'; + b = findstr(hunk, 'record_freq'); + %fseek(fd, a+b-1, -1); + %w.xxxx = char(fread(fd, 12, 'char'))'; + %w.count = fread(fd, 1, 'uint'); + %w.data_code = fread(fd, 1, 'ushort'); + fseek(fd, a+b-1+12+2+4, -1); + i.sr = fread(fd, 1, 'double'); + i.nchans = i.ndoubles+i.nfloats+i.nlongs ... + +i.nshorts+i.nchars; + recordsize = i.ndoubles*8 + i.nfloats*4 + i.nlongs*4 ... + +i.nshorts*2+i.nchars; + i.nsamples = (i.nbytes - i.bytes_to_data) / recordsize; + i.totalsamples = i.nsamples * i.nchans; + i.sample_bytes = 2; % bug + i.sample_type = 'int16'; % bug + return; +%%%%%%%%%%%%% PLAIN ASCII %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'ascii' + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + line = fgetl(fd); + i.nchans = size(sscanf(line, '%f'), 1); + nlines = 1; + while (1) + line = fgets(i.fd); + if isa(line, 'double') & -1 == line; break; end + nlines = nlines+1; + end + i.nsamples = nlines; + i.totalsamples = i.nsamples*i.nchans; + i.sr=1; + return +%%%%%%%%%%%%% COMMA-SEPARATED VALUES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case 'csv'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + % todo + return +%%%%%%%%%%%%% Binary %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +case {'uchar'}; % should take care of signed/unsigned + i.sample_bytes = 1; + i.sample_type = 'uchar'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'short', 'int16'}; + i.sample_bytes = 2; + i.sample_type = 'int16'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'long', 'int32'}; + i.sample_bytes = 4; + i.sample_type = 'int32'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'float', 'float32'}; + i.sample_bytes = 4; + i.sample_type = 'float32'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +case {'double', 'float64'}; + i.sample_bytes = 8; + i.sample_type = 'float64'; + if isfield(i, 'bytes_to_data') + fseek(fd, i.bytes_to_data, 0); + else + i.bytes_to_data = 0; + end + if ~isfield(i, 'nchans') + i.nchans = 1; + end + i.nsamples = (i.nbytes-i.bytes_to_data)/i.sample_bytes; + i.totalsamples = i.nsamples * i.nchans; + i.sr=1; + return +otherwise + error(['unknown format: >' i.format '<']); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/sf_wave.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,193 @@ +function y = sf_wave(i, samples, chans) +% y=sf_wave(i,samples,chans) - read data from file +% +% y: data read from file (column vector or matrix) +% +% i: structure containing information about file +% samples: range of samples to read ([start stop]) - default ([]) means entire file +% chans: range of channels to read ([lo hi]) - default ([]) means all channels\ +% + +% Alain de Cheveigné, CNRS/Ircam, 2002. +% Copyright (c) 2002 Centre National de la Recherche Scientifique. +% +% Permission to use, copy, modify, and distribute this software without +% fee is hereby granted FOR RESEARCH PURPOSES only, provided that this +% copyright notice appears in all copies and in all supporting +% documentation, and that the software is not redistributed for any +% fee (except for a nominal shipping charge). +% +% For any other uses of this software, in original or modified form, +% including but not limited to consulting, production or distribution +% in whole or in part, specific prior permission must be obtained from CNRS. +% Algorithms implemented by this software may be claimed by patents owned +% by CNRS, France Telecom, Ircam or others. +% +% The CNRS makes no representations about the suitability of this +% software for any purpose. It is provided "as is" without express +% or implied warranty. Beware of the bugs. + +if ~nargin | ~isfield(i, 'fname') + error('usage: i.fname = "file name", then call y=sf_wave(i)'); +end +if ~isfield(i, 'format') + i = sf_info(i); +end +% defaults +if nargin<2 | isempty(samples); + if (i.nsamples) samples=[1 i.nsamples]; end; +end +if nargin<3 | isempty(chans) + if (i.nchans) chans=[1 i.nchans]; end; +end + +% clip +if samples(1) < 1 + samples(1) = 1; +end +if samples(2) > i.nsamples + samples(2) = i.nsamples; +end +if samples(2) < samples(1) + y=[];return + error(['start sample after stop sample: ' num2str(samples(1)) '-' num2str(samples(2))]); +end +if chans(1) < 1 | chans(2) > i.nchans + error(['requested inexistent channels: ' num2str(chans(1)) '-' num2str(chans(2))]); +end + + +% workspace matrix +if strcmp(i.format, 'matrix') + y = i.fname(samples(1):samples(2), chans(1):chans(2)); + return +end + +% use matlab functions for AU and WAV and MACSND +if strcmp(i.format, 'AU') + y=auread(i.fname, [samples(1) samples(2)]); + y=y(:,chans(1):chans(2)); + return; +end +if strcmp(i.format, 'WAV') + y=wavread(i.fname, [samples(1) samples(2)]); + y=y(:,chans(1):chans(2)); + return; +end +if strcmp(i.format, 'MACSND') + if 3==exist('readsnd') + y = eval('readsnd(i.fname)'); + else + error('cannot read MACSND on this platform'); + end + y = y(samples(1):samples(2),chans(1):chans(2)); + return; +end + +% close if open +% if fopen(i.fd) +% fclose(i.fd); +% end + +if ~isfield(i, 'bytes_to_data') + i.bytes_to_data=0; +end + +% ascii formats +if strcmp(i.format, 'ascii') | strcmp(i.format, 'csv') | strcmp(i.format, 'IWAVE') + i.fd = fopen(i.fname, 'rt'); + fseek(i.fd, i.bytes_to_data, -1); + + switch i.format + case 'ascii' + nsamples = samples(2) - samples(1) + 1; + nchans = chans(2) - chans(1) + 1; + y = zeros(nsamples, nchans); + % skip to start + for j=1:samples(1)-1 + line = fgetl(i.fd); + if isempty(line) | line == -1; error('unexpected eof'); end + end + k=1; + % read + for j=samples(1) : samples(2) + line = fgetl(i.fd); + if isempty(line) | line == -1; error('unexpected eof'); end + a = sscanf(line, '%f'); + y(k,:) = a(chans(1):chans(2))'; + k = k+1; + end + case 'cvs' + error('not implemented'); + case 'IWAVE' + error('not implemented'); + end + fclose(i.fd); + return +end + +% binary formats +fr = samples(2) - samples(1) + 1; +skip_samples = i.nchans * (samples(1) - 1); +switch i.format +case 'uchar' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 1, 0); + y = fread(i.fd, [fr, i.nchans], 'uchar'); + fclose(i.fd); +case 'short' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + fclose(i.fd); +case 'long' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 4, 0); + y = fread(i.fd, [fr, i.nchans], 'long'); + fclose(i.fd); +case 'float' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 4, 0); + y = fread(i.fd, [fr, i.nchans], 'float'); + fclose(i.fd); +case 'double' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 8, 0); + y = fread(i.fd, [fr, i.nchans], 'double'); + fclose(i.fd); +case 'NIST' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + y = zeros(i.nsamples, i.nchans); + switch i.sample_coding + case 'pcm' + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + otherwise + error(['cannot handle NIST sample_coding = ', i.sample_coding]); + end + fclose(i.fd); +case 'ESPS' + i.fd = fopen(i.fname, 'r'); + fseek(i.fd, i.bytes_to_data, -1); + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + +case 'AIFF' + i.fd = fopen(i.fname, 'r', 's'); + fseek(i.fd, i.bytes_to_data, -1); + % should check sample size + fseek(i.fd, skip_samples * 2, 0); + y = fread(i.fd, [fr, i.nchans], 'short'); + fclose(i.fd); + +otherwise + error(['don''t know how to load format = ', i.format]); +end + +y = y(:,chans(1):chans(2));
Binary file Code/Descriptors/yin/yin/private/src/cumnorm_inplace.%b5 Data/CW Settings.stg has changed
Binary file Code/Descriptors/yin/yin/private/src/cumnorm_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
Binary file Code/Descriptors/yin/yin/private/src/cumnorm_inplace.%b5 Data/RESOURCE.FRK/MATLAB~1.TDT has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/cumnorm_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,92 @@ +/* +* cumnorm_inplace.c +* cumulative mean-normalization of diff function +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] + +/* Output Arguments */ + +static void cumsum_inplace( + double *xp, /* matrix to cumulative-mean-normalize (along columns) */ + int m, /* rows */ + int n /* columns */ + ) +{ + int j,k; + double z, sum, mean; + + + for (k=0; k<n; k++) { /* for each column */ + SET(xp[k*m+0])=1; + sum = 0.0; + for (j=1;j<m;j++) { /* for each row */ + z = GET(xp[k*m+j]); + z = (z > 0) ? z:0; /* clip to remove numerical artifacts */ + sum += z; + z = (sum>0) ? (z / (sum/j)) : 1; /* cumulative-mean-normalize */ + SET(xp[k*m+j])=z; + } + } +/* + for (j=0;j<m;j++) { // for each row (time) + sum=0; + SET(xp[j])=1; + for (k=1; k<n; k++) { // for each column (lag) + z=GET(xp[m*k+j]); + sum+=z; + mean=sum/k; + SET(xp[m*k+j]) = (mean > 0) ? z/mean : 1; + } + } +*/ + + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp; + int nx, mx; + + /* Check for proper number of arguments */ + if (nrhs != 1) { + mexErrMsgTxt("CUMSUM_INPLACE takes 1 input argument"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("CUMSUM_INPLACE: X should be doubles"); + } + mx=mxGetM(X_IN); /* rows */ + nx=mxGetN(X_IN); /* columns */ + if (nx*mx == 0) { + mexErrMsgTxt("CUMSUM_INPLACE: input matrix is empty"); + } + + xp = mxGetPr(X_IN); + checkin_matrix((mxArray *) X_IN); + + + /* Do the actual computations in a subroutine */ + cumsum_inplace(xp,mx,nx); + + checkout_matrix((mxArray *) X_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/cumnorm_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/yin/private/src/dftoperiod.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/dftoperiod.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,144 @@ +/* +* dftoperiod.c +* estimate period from df +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define D_IN prhs[0] +#define B_IN prhs[1] +#define T_IN prhs[2] + +/* output arguments */ +#define PRD_OUT plhs[0] + +static void ddftoperiods( + double *dp, /* input df */ + double *bp, /* bounds */ + double *prdp, /* periods */ + double thresh, /* threshold */ + int m, /* nrows of df */ + int n /* ncols of df */ + ) +{ + int j,k, p, lo, hi, maxj, goodflag; + double z, min, globalmin; + + /* bounds */ + lo= GET(bp[0]); + hi= GET(bp[1]); + if (lo<0 || hi>m ) { + mexPrintf("%d %d\n", lo, hi); + mexErrMsgTxt("DFTOPERIOD: bounds are out of bounds"); + } + + #define MARGIN 1.5 + + for (k=0;k<n;k++) { /* for each col */ + + goodflag=0; + + /* estimate global min + globalmin=GET(dp[k*m+lo]); + for (j=lo;j<hi; j++) { + z = GET(dp[k*m+j]); + if (z<globalmin) { + globalmin = z; + } + } + + thresh += globalmin; */ + + /* estimate period */ + p=lo; + min=GET(dp[k*m+lo]); + maxj=hi; + for (j=lo;j<maxj; j++) { /* for each row */ + + z = GET(dp[k*m+j]); + + if (z<thresh && z <= min) { /* good candidate, restrict search range */ + maxj=j*MARGIN; + if (maxj> hi) { maxj = hi; } + } + + if (z<min) { /* update min & period estimate */ + min = z; + p=j; + } + + } + + /* mexPrintf("%f %d %f\n", thresh, p, min); */ + + SET(prdp[k])=p; + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *dp, *bp, *tp, *prdp, thresh; + int m, n; + + /* Check for proper number of arguments */ + if (nrhs != 3){ + mexErrMsgTxt("DFTOPERIOD requires 3 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(D_IN) || mxIsComplex(D_IN) || + mxIsSparse(D_IN) || !mxIsDouble(D_IN) ) { + mexErrMsgTxt("DFTOPERIOD: D should be doubles"); + } + if (!mxIsNumeric(B_IN) || mxIsComplex(B_IN) || + mxIsSparse(B_IN) || !mxIsDouble(B_IN) ) { + mexErrMsgTxt("DFTOPERIOD: B should be doubles"); + } + if (!mxIsNumeric(T_IN) || mxIsComplex(T_IN) || + mxIsSparse(T_IN) || !mxIsDouble(T_IN) ) { + mexErrMsgTxt("DFTOPERIOD: T should be double"); + } + m=mxGetM(D_IN); /* rows */ + n=mxGetN(D_IN); /* columns */ + + if (mxGetM(B_IN) * mxGetN(B_IN) != 2) { + mexErrMsgTxt("DFTOPERIOD: B should be a size 2 vector"); + } + if (mxGetM(T_IN)*mxGetN(T_IN) != 1) { + mexErrMsgTxt("DFTOPERIOD: T should be scalar"); + } + + /* Create matrix to return */ + PRD_OUT = mxCreateDoubleMatrix(1,n, mxREAL); + + dp = mxGetPr(D_IN); + bp = mxGetPr(B_IN); + tp = mxGetPr(T_IN); thresh=tp[0]; + prdp = mxGetPr(PRD_OUT); + + checkin_matrix((mxArray *) D_IN); + checkin_matrix((mxArray *) B_IN); + checkin_matrix((mxArray *) PRD_OUT); + + /* Do the actual computations in a subroutine */ + ddftoperiods(dp,bp,prdp,thresh,m,n); + + checkout_matrix((mxArray *) D_IN); + checkout_matrix((mxArray *) B_IN); + checkout_matrix((mxArray *) PRD_OUT); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/dftoperiod.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/yin/private/src/dftoperiod2.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/dftoperiod2.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,144 @@ +/* +* dftoperiods.c +* estimate period from df +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define D_IN prhs[0] +#define B_IN prhs[1] +#define T_IN prhs[2] + +/* output arguments */ +#define PRD_OUT plhs[0] + +static void ddftoperiods( + double *dp, /* input df */ + double *bp, /* bounds */ + double *prdp, /* periods */ + double thresh, /* threshold */ + int m, /* nrows of df */ + int n /* ncols of df */ + ) +{ + int j,k, p, lo, hi, maxj, goodflag; + double z, min, globalmin; + + /* bounds */ + lo= GET(bp[0]); + hi= GET(bp[1]); + if (lo<0 || hi>m ) { + mexPrintf("%d %d\n", lo, hi); + mexErrMsgTxt("DFTOPERIOD: bounds are out of bounds"); + } + + #define MARGIN 1.5 + + for (k=0;k<n;k++) { /* for each col */ + + goodflag=0; + + /* estimate global min */ + globalmin=GET(dp[k*m+lo]); + for (j=lo;j<hi; j++) { /* for each row */ + z = GET(dp[k*m+j]); + if (z<globalmin) { /* update min & period estimate */ + globalmin = z; + } + } + + thresh += globalmin; + + /* estimate period */ + p=lo; + min=GET(dp[k*m+lo]); + maxj=hi; + for (j=lo;j<maxj; j++) { /* for each row */ + + z = GET(dp[k*m+j]); + + if (z<thresh && z <= min) { /* good candidate, restrict search range */ + maxj=j*MARGIN; + if (maxj> hi) { maxj = hi; } + } + + if (z<min) { /* update min & period estimate */ + min = z; + p=j; + } + + } + + /* mexPrintf("%f %d %f\n", thresh, p, min); */ + + SET(prdp[k])=p; + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *dp, *bp, *tp, *prdp, thresh; + int m, n; + + /* Check for proper number of arguments */ + if (nrhs != 3){ + mexErrMsgTxt("DFTOPERIOD requires 3 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(D_IN) || mxIsComplex(D_IN) || + mxIsSparse(D_IN) || !mxIsDouble(D_IN) ) { + mexErrMsgTxt("DFTOPERIOD: D should be doubles"); + } + if (!mxIsNumeric(B_IN) || mxIsComplex(B_IN) || + mxIsSparse(B_IN) || !mxIsDouble(B_IN) ) { + mexErrMsgTxt("DFTOPERIOD: B should be doubles"); + } + if (!mxIsNumeric(T_IN) || mxIsComplex(T_IN) || + mxIsSparse(T_IN) || !mxIsDouble(T_IN) ) { + mexErrMsgTxt("DFTOPERIOD: T should be double"); + } + m=mxGetM(D_IN); /* rows */ + n=mxGetN(D_IN); /* columns */ + + if (mxGetM(B_IN) * mxGetN(B_IN) != 2) { + mexErrMsgTxt("DFTOPERIOD: B should be a size 2 vector"); + } + if (mxGetM(T_IN)*mxGetN(T_IN) != 1) { + mexErrMsgTxt("DFTOPERIOD: T should be scalar"); + } + + /* Create matrix to return */ + PRD_OUT = mxCreateDoubleMatrix(1,n, mxREAL); + + dp = mxGetPr(D_IN); + bp = mxGetPr(B_IN); + tp = mxGetPr(T_IN); thresh=tp[0]; + prdp = mxGetPr(PRD_OUT); + + checkin_matrix((mxArray *) D_IN); + checkin_matrix((mxArray *) B_IN); + checkin_matrix((mxArray *) PRD_OUT); + + /* Do the actual computations in a subroutine */ + ddftoperiods(dp,bp,prdp,thresh,m,n); + + checkout_matrix((mxArray *) D_IN); + checkout_matrix((mxArray *) B_IN); + checkout_matrix((mxArray *) PRD_OUT); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/dftoperiod2.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/yin/private/src/interp_inplace.%b5 Data/CW Settings.stg has changed
Binary file Code/Descriptors/yin/yin/private/src/interp_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/interp_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,101 @@ +/* +* interp_inplace.c +* linear interpolation +* +* Alain de Cheveigné, CNRS/Ircam +* (c) 2003 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] +#define IDX_IN prhs[1] + +/* Output Arguments */ + +static void interp_inplace( + double *xp, /* input vector */ + double *idxp, /* index vector */ + int m, /* rows input vector */ + int midx /* rows index vector */ + ) +{ + int j,idxi; + double a,b,idx,idxf; + + + for (j=0;j<midx;j++) { /* for each index */ + idx=GET(idxp[j]); + idxi=floor(idx); + idxf=idx-idxi; + if (idxi<0) { + a=GET(xp[0]); + b=GET(xp[1]); + SET(idxp[j]) = a + idx*(b-a); + } else if (idxi>=m) { + a=GET(xp[m-2]); + b=GET(xp[m-1]); + SET(idxp[j]) = b + (idx-m+1)*(b-a); + } else { + a=GET(xp[idxi]); + b=GET(xp[idxi+1]); + SET(idxp[j]) = a + idxf*(b-a); + } + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *idxp; + int nx, mx, nidx, midx; + + /* Check for proper number of arguments */ + if (nrhs !=2 ) { + mexErrMsgTxt("INTERP_INPLACE takes 2 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("INTERP_INPLACE: X should be doubles"); + } + if (!mxIsNumeric(IDX_IN) || mxIsComplex(IDX_IN) || + mxIsSparse(IDX_IN) || !mxIsDouble(IDX_IN) ) { + mexErrMsgTxt("INTERP_INPLACE: Y should be doubles"); + } + mx=mxGetM(X_IN); /* rows */ + nx=mxGetN(X_IN); /* columns */ + midx=mxGetM(IDX_IN); + nidx=mxGetN(IDX_IN); + + if (nx>1 || nidx>1) { + mexErrMsgTxt("INTERP_INPLACE: X and IDX should be column vectors"); + } + if (mx<1) { + mexErrMsgTxt("INTERP_INPLACE: X should have at least two samples"); + } + + xp = mxGetPr(X_IN); + idxp = mxGetPr(IDX_IN); + + checkin_matrix((mxArray *) IDX_IN); + checkin_matrix((mxArray *) X_IN); + + /* Do the actual computations in a subroutine */ + interp_inplace(xp,idxp,mx,midx); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix((mxArray *) IDX_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/interp_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/yin/private/src/mininrange.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/mininrange.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,119 @@ +/* +* minInRange.c +* for each sample, return index of the smallest value within an interval +* of that sample +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ + +#define X_IN prhs[0] +#define INTERVAL_IN prhs[1] + +/* Output Arguments */ + +#define IDX_OUT plhs[0] + +static void mininrange( + double *xp, + double *intp, + double *idxp, + unsigned int n + ) +{ + int h,k, idx, left, right, interval; + double min; + double max; + + + for (k=0; k<n; k++) { + interval=(int) GET(intp[k]); + + left = k - interval/2; + right = left+interval; + + if (left<0) left=0; /* clip */ + if (right>n) right=n; + + min = GET(xp[k]); + idx = k; + for (h=left;h<right; h++) { + if (GET(xp[h]) < min) { /* update min */ + min = GET(xp[h]); + idx = h; + } + } + SET(idxp[k])=idx+1; + } + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *idxp, *intp; + unsigned int m,n; + + /* Check for proper number of arguments */ + + if (nrhs != 2) { + mexErrMsgTxt("MININRANGE requires two input arguments"); + } else if (nlhs != 1) { + mexErrMsgTxt("MININRANGE requires one output argument"); + } + + /* Check type of input */ + + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN)) { + mexErrMsgTxt("MININRANGE: X should be matrix of doubles"); + } + if (!mxIsNumeric(INTERVAL_IN) || mxIsComplex(INTERVAL_IN) || + mxIsSparse(INTERVAL_IN) || !mxIsDouble(INTERVAL_IN)) { + mexErrMsgTxt("MININRANGE X should be matrix of doubles"); + } + + m = mxGetM(X_IN); /* rows */ + n = mxGetN(X_IN); /* columns */ + if (m>1 || n <=1) { + mexErrMsgTxt("MININRANGE X should be row vector"); + } + if (m != mxGetM(INTERVAL_IN) || n != mxGetN(INTERVAL_IN)) { + mexErrMsgTxt("MININRANGE: INTERVAL should be of same size as X"); + } + + /* Create matrix to return */ + + IDX_OUT = mxCreateDoubleMatrix(1, n, mxREAL); + + /* Assign pointers to the various parameters */ + + xp = mxGetPr(X_IN); + intp = mxGetPr(INTERVAL_IN); + idxp = mxGetPr(IDX_OUT); + + checkin_matrix((mxArray *) X_IN); + checkin_matrix((mxArray *) INTERVAL_IN); + checkin_matrix(IDX_OUT); + + /* Do the actual computations in a subroutine */ + + mininrange(xp,intp,idxp,n); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix((mxArray *) INTERVAL_IN); + checkout_matrix(IDX_OUT); + return; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/mininrange.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/yin/private/src/minparabolic.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/minparabolic.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,121 @@ +/* +* minInterp.c +* determine position of minimum of parabolic fit to local minima +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] + + +/* Output Arguments */ +#define MIN_OUT plhs[0] +#define IND_OUT plhs[1] + + +static void mininterp( + double *xp, + unsigned int m, + unsigned int n, + double *minp, + double *indp + ) +{ + double x1,x2,x3,a,b,min,ind; + int i,j,k; + + /* each column is interpolated independently */ + /* each local minimum is replaced by the value of the interpolated minimum, and the + position of the interplolated minimum is recorded in indp */ + + for (j=0; j<n; j++) { + /* first row: don't touch */ + SET(indp[j*m])=0; + SET(minp[j*m])=GET(xp[j*m]); + /* intermediate: refine */ + for (k=1; k<m-1; k++) { + i=k+j*m; + x1=GET(xp[i-1]); x2=GET(xp[i]); x3=GET(xp[i+1]); + SET(indp[i]) = 0; /* default */ + SET(minp[i]) = x2; + if ((x2<=x1) && (x2<=x3)) { + a = (x1 + x3 - 2*x2)/2; + b = (x3 - x1)/2; + if (a) { + SET(indp[i]) = - b / (2*a); + SET(minp[i]) = x2 + - (b*b) / (4*a); + } + } + } + /* last row: don't touch */ + SET(indp[m-1 + j*m])=0; + SET(minp[m-1 + j*m])=GET(xp[m-1 + j*m]); + } + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *minp, *indp, *xp; + unsigned int m,n; + + /* Check for proper number of arguments */ + + if (nrhs != 1) { + mexErrMsgTxt("minParabolic requires one input argument."); + } else if (nlhs > 2) { + mexErrMsgTxt("minParabolic requires one or two output arguments"); + } + + + /* Check type of input */ + + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN)) { + mexErrMsgTxt("minParabolic requires that X be a matrix of doubles"); + } + + m = mxGetM(X_IN); /* rows */ + n = mxGetN(X_IN); /* columns */ + if (m<3) { mexErrMsgTxt("number of rows should be 3 or more"); } + + + /* Create two matrices to return */ + + MIN_OUT = mxCreateDoubleMatrix(m, n, mxREAL); + IND_OUT = mxCreateDoubleMatrix(m, n, mxREAL); + + + /* Assign pointers to the various parameters */ + + minp = mxGetPr(MIN_OUT); + indp = mxGetPr(IND_OUT); + xp = mxGetPr(X_IN); + + checkin_matrix((mxArray *) X_IN); + checkin_matrix(MIN_OUT); + checkin_matrix(IND_OUT); + + /* Do the actual computations in a subroutine */ + + + mininterp(xp,m,n,minp,indp); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix(MIN_OUT); + checkout_matrix(IND_OUT); + return; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/minparabolic.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/mwrap.h Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,268 @@ +/* +19 Dec 94 - Sept 95 +Alain de Cheveigne, LLF, CNRS/Universite Paris 7. +Modified: Sept 2000 + +file mwrap.h + +To use the malloc wrappers: +- include this file in all source files that call malloc routines, +- link with mwrap.o, +- compile all files with -DMWRAP to turn on wrapping and checking. + +*/ + +#ifndef MWRAP_H +#define MWRAP_H + +#include <stdio.h> +#include <stdlib.h> + +/* NOTE: (char *) is used for pointer arithmetic. + Pointer arguments of all macros are cast to (char *). */ + +/* node structure for binary splay tree */ +typedef struct mwrap_node_struct *mwrap_node ; +struct mwrap_node_struct { + char *lo; /* base of block */ + char *hi; /* one beyond last address in block */ + mwrap_node up; + mwrap_node left; + mwrap_node right; +}; + + +/* Abbreviations for macros. Change in case of conflict. + + Assignment macros: + GET(x) equivalent to 'x' on a rhs, but memory address is checked before access + SET(x) equivalent to 'x' on a lhs, but memory address is checked before assignment + CPY(p,q,n) copy n values starting from p to n values starting from q, after checking both + + Test macros: + POK(p) pointer p is within an allocated block. + BOK(b) pointer b is the base of an allocated block. + PQOK(p, q) pointers p and q are in same allocated block. + PBOK(p, b) pointer p is in block of base b. + + MWRAPOK() check internal tree for sanity (slow) + + Query macros: + PBASE(p) base of block containing p, NULL if there is none. + PBUMP(p) one beyond last byte of block. + PSIZE(p) size of block containing p, -1 if none. + + Explicit block registering macros. + CHECKIN(lo, hi) register block of bounds lo, hi + CHECKOUT(lo) unregister block of base lo + +*/ + +#define GET(x) MWRAP_GET(x) +#define SET(x) MWRAP_SET(x) +#define CPY(p,q,n) MWRAP_CPY(p,q,n) +#define POK(p) MWRAP_POK(p) +#define BOK(b) MWRAP_BOK(b) +#define PQOK(p,q) MWRAP_PQOK((p),(q)) +#define PBOK(p,b) MWRAP_PBOK((p),(b)) +#define MWRAPOK() MWRAP_MWRAPOK() +#define PBASE(p) MWRAP_PBASE(p) +#define PSIZE(p) MWRAP_PSIZE(p) +#define CHECKIN(lo, hi) MWRAP_CHECKIN((lo), (hi)) +#define CHECKOUT(lo) MWRAP_CHECKOUT(lo) + + +/* Definitions of macros. + All are conditional on MWRAP being defined. + For speed, macros check for a "lucky" hit (block already + at root) before handing over to functions. */ + +#ifdef MWRAP + +#define MWRAP_GET(x) \ + ( (MWRAP_POK(&x)) ? (x) : (x) ) + +#define MWRAP_SET(x) \ + MWRAP_POK(&x) ; (x) + +#define MWRAP_CPY(p, q, n) \ + do {int i; \ + MWRAP_PQOK((p),(p)+(n-1)); \ + MWRAP_PQOK((q),(q)+(n-1)); \ + for(i=0;i<n;i++) (q)[i]=(p)[i]; \ + } while(0) + +/* #define MWRAP_POK(p) \ + do { if (mwrap_tree \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->hi) break; \ + mwrap_pok((char *) (p), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_POK(p) \ + ((mwrap_tree && (char *) (p) >= mwrap_tree->lo && (char *) (p) <mwrap_tree->hi) ? \ + 1 : mwrap_pok((char *) (p), __FILE__, __LINE__)) + + +/* #define MWRAP_BOK(b) \ + do { if (mwrap_tree \ + && (char *) (b) == mwrap_tree->lo) break; \ + mwrap_bok((char *) (b), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_BOK(b) \ + ((mwrap_tree && (char *) (b) == mwrap_tree->lo) ? 1 : mwrap_bok((char *) (b), __FILE__, __LINE__)) + + +/* #define MWRAP_PQOK(p, q) \ + do { if (mwrap_tree \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->hi \ + && (char *) (q) >= mwrap_tree->lo \ + && (char *) (q) < mwrap_tree->lo) break; \ + mwrap_pqok((char *)(p), (char *)(q), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_PQOK(p,q) \ + ((mwrap_tree \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->hi \ + && (char *) (q) >= mwrap_tree->lo \ + && (char *) (q) < mwrap_tree->lo ) ? \ + 1 : mwrap_pqok((char *)(p), (char *)(q), __FILE__, __LINE__)) + +/* #define MWRAP_PBOK(p, b) \ + do { if (mwrap_tree \ + && (char *) (b) == mwrap_tree->lo \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->lo) break; \ + mwrap_pbok((char *)(p), (char *)(b), __FILE__, __LINE__); } while (0) */ + +#define MWRAP_PBOK(p, b) \ + (( mwrap_tree \ + && (char *) (b) == mwrap_tree->lo \ + && (char *) (p) >= mwrap_tree->lo \ + && (char *) (p) < mwrap_tree->lo) ? \ + 1 : mwrap_pbok((char *)(p), (char *)(b), __FILE__, __LINE__)) + +#define MWRAP_MWRAPOK() mwrap_checktree(__FILE__, __LINE__) + +#define MWRAP_PBASE(p) mwrap_pbase((char *) p) +#define MWRAP_PBUMP(p) mwrap_pbump((char *) p) +#define MWRAP_PSIZE(p) mwrap_psize((char *) p) +#define MWRAP_CHECKIN(lo, hi) \ + mwrap_checkin((char *) (lo), (char *) (hi), __FILE__, __LINE__) +#define MWRAP_CHECKOUT(lo) \ + mwrap_checkout( (char *) (lo), __FILE__, __LINE__) + +#else + +#define MWRAP_GET(x) (x) +#define MWRAP_SET(x) (x) +#define MWRAP_CPY(p,q,n) +#define MWRAP_POK(p) +#define MWRAP_BOK(p) +#define MWRAP_PBOK(p,b) +#define MWRAP_PQOK(p,q) +#define MWRAP_MWRAPOK() +#define MWRAP_PBASE(p) (void *) mwrap_squeal(__FILE__, __LINE__) +#define MWRAP_PBUMP(p) (void *) mwrap_squeal(__FILE__, __LINE__) +#define MWRAP_PSIZE(p) (long) mwrap_squeal(__FILE__, __LINE__) +#define MWRAP_CHECKIN(lo, hi) +#define MWRAP_CHECKOUT(lo) +#endif + + +/* Wrapping macros. + + If MWRAP is defined, all malloc routines are wrapped and the + original routines are available as "nowrap_malloc()", etc.. + + If MWRAP is not defined, malloc routines are not wrapped + and "nowrap_malloc()", etc. are #defined to "malloc()", etc. +*/ + +#ifdef MWRAP + +#ifndef MWRAP_C +#define malloc(size) mwrap_malloc((size), __FILE__, __LINE__) +#define free(p) mwrap_free((p), __FILE__, __LINE__) +#define realloc(p, size) mwrap_realloc((p), (size), __FILE__, __LINE__) +#define calloc(n, size) mwrap_calloc((n), (size), __FILE__, __LINE__) +#define cfree(p) mwrap_cfree((p), __FILE__, __LINE__) +#define valloc(size) mwrap_valloc((size), __FILE__, __LINE__) +#define vfree(p) mwrap_vfree((p), __FILE__, __LINE__) +#define memalign(a, p) mwrap_memalign((a), (p), __FILE__, __LINE__) +#endif /* MWRAP_C */ + +#else + +#ifndef MWRAP_C +#define nowrap_malloc(size) malloc(size) +#define nowrap_free(p) free(p) +#define nowrap_realloc(p, size) realloc(p) +#define nowrap_calloc(n, size) calloc((n), (size)) +#define nowrap_cfree(p) cfree(p) +#define nowrap_valloc(size) valloc(size) +#define nowrap_vfree(p) vfree(p) +#define nowrap_memalign(a, p) memalign((a), (p)) +#endif /* MWRAP_C */ + +#endif /* MWRAP */ + + +/* three global variables */ + +#ifndef MWRAP_C +extern mwrap_node mwrap_tree; /* the tree */ +extern char *mwrap_file; /* caller source file */ +extern long mwrap_line; /* line in caller source file */ +#endif + + +/* Data types used by malloc library. Modify if necessary */ +typedef void DATATYPE; /* returned by malloc(), etc. */ +typedef size_t SIZETYPE; /* fed to malloc(), etc. */ +typedef void VOIDTYPE; /* returned by free(), etc. */ + + +/* prototypes of routines made public */ + +/* malloc wrappers */ +DATATYPE *mwrap_malloc(SIZETYPE size, char *file, long line); +VOIDTYPE mwrap_free(DATATYPE *p, char *file, long line); +DATATYPE *mwrap_realloc(DATATYPE *p, SIZETYPE size, char *file, long line); +DATATYPE *mwrap_calloc(SIZETYPE n, SIZETYPE size, char *file, long line); +VOIDTYPE mwrap_cfree(DATATYPE *p, char *file, long line); +DATATYPE *mwrap_valloc(SIZETYPE size, char *file, long line); +VOIDTYPE mwrap_vfree(DATATYPE *p, char *file, long line); +DATATYPE *mwrap_memalign(SIZETYPE alignment, SIZETYPE size, char *file, long line); + +/* unwrapped versions */ +#ifdef MWRAP +DATATYPE *nowrap_malloc(SIZETYPE size); +VOIDTYPE nowrap_free(DATATYPE *p); +DATATYPE *nowrap_realloc(DATATYPE *p, SIZETYPE size); +DATATYPE *nowrap_calloc(SIZETYPE n, SIZETYPE size); +VOIDTYPE nowrap_cfree(DATATYPE *p); +DATATYPE *nowrap_valloc(SIZETYPE size); +VOIDTYPE nowrap_vfree(DATATYPE *p); +DATATYPE *nowrap_memalign(SIZETYPE alignment, SIZETYPE size); +#endif + +/* explicit block registering */ +void mwrap_checkin(char *lo, char *hi, char *file, long line); +void mwrap_checkout(char *lo, char *file, long line); + +/* tests & queries */ +int mwrap_pok(char *p, char *file, long line); +int mwrap_bok(char *b, char *file, long line); +int mwrap_pqok(char *p, char *q, char *file, long line); +int mwrap_pbok(char *p, char *b, char *file, long line); +char *mwrap_pbase(char *p); +char *mwrap_pbump(char *p); +long mwrap_psize(char *p); +long mwrap_squeal(char *file, long line); + +/* change default functions */ +void mwrap_set_errorfunction(void (*f)(int level)); +void mwrap_set_messagefunction(void (*f)(int level, char * message)); + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/mwrap_check.h Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,55 @@ +/* +* stuff to allow checking with mwrap +* +*/ + +#include "mwrap.h" +#include "string.h" + +#define MACINTOSH /* compiler (?) converts \n to \r: convert them back */ + +void checkin_matrix(mxArray *m); +void checkout_matrix(mxArray *m); +void mex_messagefunction(int level, char* message); +void mex_errorfunction(int level); + +/* check matrix into mwrap's tree */ +void checkin_matrix(mxArray *m) +{ + char *base, *top; + + base = (char *) mxGetPr(m); + top = base + mxGetN(m) * mxGetM(m) * sizeof(double); + /* mexPrintf("%d %d\n", base, top); */ + CHECKIN(base, top); +} + +/* check matrix out of mwrap's tree */ +void checkout_matrix(mxArray *m) +{ + char *base; + + base = (char *) mxGetPr(m); + CHECKOUT(base); +} + +/* message function to give to mwrap */ +void mex_messagefunction(int level, char* message) { +#ifdef MACINTOSH + {char *c; + c = strchr(message, (int) '\r'); + + while(c) { + *c = '\n'; + c = strchr(c, (int) '\r'); + } + } +#endif + mexPrintf(message); +} + +/* how to exit mex function */ +void mex_errorfunction(int level) { + /* mexPrintf("\n"); */ + mexErrMsgTxt(" "); +}
Binary file Code/Descriptors/yin/yin/private/src/rdiff_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/rdiff_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,170 @@ +/* +* rdiff_inplace.c +* running difference function +* +* Alain de Cheveigné, CNRS/Ircam Dec 2001 +* (c) 2001 CNRS +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] +#define Y_IN prhs[1] +#define R_IN prhs[2] +#define L_IN prhs[3] +#define N_IN prhs[4] + +/* Output Arguments */ + +static void rdiff_inplace( + double *xp, /* input vector 1 */ + double *yp, /* input vector 2 */ + double *rp, /* df matrix (time X lag) */ + double *lp, /* lag matrix */ + int N, /* size of summation window */ + int m, /* rows in corr matrix */ + int n /* columns in corr matrix */ + ) +{ + int j,k,i, xlag, ylag; + double z, d, *r, *x, *y, *bump; + + + for (j=0;j<n;j++) { /* for each column (lag pair) */ + xlag = GET(lp[2*j]); + ylag = GET(lp[2*j+1]); + if (N==1) { + r=&rp[j*m]; + x=&xp[xlag]; + y=&yp[ylag]; + bump=r+m; + while(r<bump) { /* for each row */ + d = GET(*x)-GET(*y); + SET(*r) = d*d; + r++; x++; y++; + } + } else { + for (k=0; k<m; k++) { /* for each row */ + z=0; + x=&xp[xlag+k*N]; + y=&yp[ylag+k*N]; + bump=x+N; + while(x<bump) { + d = GET(*x)-GET(*y); + z+=d*d; + x++; y++; + } + SET(rp[j*m+k]) = z; + } + } + } + + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *yp, *rp, *lp, *up, *np; + int nx, mx, ny, my, nr, mr, nl, ml, xmax, ymax, j, N, xlag, ylag; + + /* Check for proper number of arguments */ + if (nrhs < 4 || nrhs > 5) { + mexErrMsgTxt("RDIFF_INPLACE takes 4 or 5 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: X should be doubles"); + } + if (!mxIsNumeric(Y_IN) || mxIsComplex(Y_IN) || + mxIsSparse(Y_IN) || !mxIsDouble(Y_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: Y should be doubles"); + } + if (!mxIsNumeric(R_IN) || mxIsComplex(R_IN) || + mxIsSparse(R_IN) || !mxIsDouble(R_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: R should be doubles"); + } + if (nrhs==5) { + if (!mxIsNumeric(N_IN) || mxIsComplex(N_IN) || + mxIsSparse(N_IN) || !mxIsDouble(N_IN) ) { + mexErrMsgTxt("RDIFF_INPLACE: N should double"); + } + np = mxGetPr(N_IN); + N=*np; + } else { + N=1; + } + mx=mxGetM(X_IN); /* rows */ + nx=mxGetN(X_IN); /* columns */ + my=mxGetM(Y_IN); + ny=mxGetN(Y_IN); + mr=mxGetM(R_IN); + nr=mxGetN(R_IN); + ml=mxGetM(L_IN); + nl=mxGetN(L_IN); + + if (nx>1 || ny>1) { + mexErrMsgTxt("RDIFF_INPLACE: X and Y should be column vectors"); + } + if (nr !=nl) { + mexErrMsgTxt("RDIFF_INPLACE: result and lag matrices should have same number of columns"); + } + if ( ml != 2) { + mexErrMsgTxt("RDIFF_INPLACE: lags should be a two row matrix"); + } + if (!nx || !mx || !ny || !my || !nr || !mr || !nl || !ml) { + mexErrMsgTxt("RDIFF_INPLACE: one of the input matrices is empty"); + } + + xp = mxGetPr(X_IN); + yp = mxGetPr(Y_IN); + rp = mxGetPr(R_IN); + lp = mxGetPr(L_IN); + + checkin_matrix((mxArray *) L_IN); + + /* find maximum lag for x and y and check that no lag is negative */ + xmax=GET(lp[0]); + ymax=GET(lp[1]); + for (j=1;j<nl;j++) { + xlag = GET(lp[2*j]); + ylag = GET(lp[2*j+1]); + if ( xlag<0 || ylag<0 ) { + mexErrMsgTxt("RDIFF_INPLACE: LAGS should be non-negative"); + } + if (xlag>xmax ) xmax=xlag; + if (ylag>ymax ) ymax=ylag; + } + + if ( (N*mr+ymax)> my) { + mexPrintf("mr*N: %d, ymax: %d, my: %d\n", mr*N, ymax, my); + mexErrMsgTxt("RDIFF_INPLACE: Y data too short"); + } + if ( (N*mr+xmax)> mx) { + mexPrintf("mr*N: %d, xmax: %d, mx: %d\n", mr*N, xmax, mx); + mexErrMsgTxt("RDIFF_INPLACE: X data too short"); + } + + checkin_matrix((mxArray *) X_IN); + if (yp!=xp) checkin_matrix((mxArray *) Y_IN); + checkin_matrix((mxArray *) R_IN); + + /* Do the actual computations in a subroutine */ + rdiff_inplace(xp,yp,rp,lp,N,mr,nr); + + checkout_matrix((mxArray *) X_IN); + if (yp!=xp) checkout_matrix((mxArray *) Y_IN); + checkout_matrix((mxArray *) R_IN); + checkout_matrix((mxArray *) L_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/rdiff_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/yin/private/src/rsmooth.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/rsmooth.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,153 @@ +/* +* rsmooth.c +* smooth columnwise by convolution by a running square window +* multiple passes allowed (--> triangular, gaussian, etc.) +* +* Alain de Cheveigné, CNRS Ircam, Dec 2001. +* (c) CNRS 2001. +*/ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ + +#define X_IN prhs[0] +#define SMOOTH_IN prhs[1] +#define NPASSES_IN prhs[2] +#define CLIP_IN prhs[3] + +/* Output Arguments */ + +#define Y_OUT plhs[0] + +static void rsmooth( + double *xp, + int smooth, + int npasses, + double *yp, + int m, + int n, + int clip + ) +{ + int h,i,j,k,mm; + double *work1, *work2, *a, *b, *tmp, sum; + + mm = m + smooth + npasses*(smooth-1); /* nrows of work buffer */ + work1 = (double *) mxCalloc(mm*n,sizeof(double)); + work2 = (double *) mxCalloc(mm*n,sizeof(double)); + + CHECKIN(work1, (char *) work1 + sizeof(double)*mm*n); + CHECKIN(work2, (char *) work2 + sizeof(double)*mm*n); + + a = work1; + b = work2; + + /* transfer data to buffer, preceded with leading zeros */ + for (k=0; k<n; k++) { + for (j=0;j<m;j++) { + SET(a[ j+smooth+mm*k ]) = GET(xp[ j+m*k ]); + } + } + + /* smooth repeatedly */ + for (h=0;h<npasses;h++) { + for (k=0; k<n; k++) { /* for each column */ + sum=0; + for (j=smooth;j<mm;j++) { /* for each row (inner loop) */ + sum += - GET(a[j-smooth+mm*k]) + GET(a[j+mm*k]); + SET(b[j+mm*k]) = sum/smooth; + } + } + tmp=a; a=b; b=tmp; /* swap for next round */ + } + + /* transfer to output matrix */ + if (clip) { + for (j=0;j<m;j++) { + for (k=0; k<n; k++) { + SET(yp[j+m*k]) = GET(a[ j+(int)(npasses*(smooth-1)/2)+smooth + mm*k]); + } + } + } else { + for (j=0;j<m+npasses*(smooth-1);j++) { + for (k=0; k<n; k++) { + SET(yp[j+(m+npasses*(smooth-1))*k]) = GET(a[j + smooth + mm*k]); + } + } + } + CHECKOUT(work1); + CHECKOUT(work2); + +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *yp; + int m,n,npasses,smooth,clip; + + /* Check for proper number of arguments */ + if (nrhs > 4 || nrhs < 2) { + mexErrMsgTxt("RSMOOTH takes at least 2 and most 4 input arguments"); + } else if (nlhs > 1) { + mexErrMsgTxt("RSMOOTH allows at most 1 output argument"); + } + + /* Check type of input */ + m = mxGetM(X_IN); /* rows */ + n = mxGetN(X_IN); /* columns */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) || m<=1) { + mexErrMsgTxt("RSMOOTH: X should be a column vector or matrix of doubles"); + } + + if (mxGetM(SMOOTH_IN)*mxGetN(SMOOTH_IN) > 1 || mxGetPr(SMOOTH_IN)[0] <1) { + mexErrMsgTxt("RSMOOTH: SMOOTH should be a positive scalar"); + } + smooth = mxGetPr(SMOOTH_IN)[0]; + + if (nrhs > 2) { + npasses = mxGetPr(NPASSES_IN)[0]; + if (mxGetM(NPASSES_IN)*mxGetN(NPASSES_IN) > 1 || npasses <1) { + mexErrMsgTxt("RSMOOTH: NPASSES should be a positive scalar"); + } + } else { + npasses=1; + } + + if (nrhs > 3) { + clip=1; + } else { + clip=0; + } + + /* Create matrix to return */ + if (clip) { + Y_OUT = mxCreateDoubleMatrix(m, n, mxREAL); + } else { + Y_OUT = mxCreateDoubleMatrix(m+npasses*(smooth-1), n, mxREAL); + } + + /* Assign pointers to the various parameters */ + xp = mxGetPr(X_IN); + yp = mxGetPr(Y_OUT); + + checkin_matrix((mxArray *) X_IN); + checkin_matrix(Y_OUT); + + /* Do the actual computations in a subroutine */ + rsmooth(xp,smooth,npasses,yp, m,n,clip); + + checkout_matrix((mxArray *) X_IN); + checkout_matrix(Y_OUT); + return; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/rsmooth.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
Binary file Code/Descriptors/yin/yin/private/src/rsum_inplace.%b5 Data/MATLAB PPC MEX-file.tdt has changed
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/rsum_inplace.c Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,134 @@ +/* +* rsum_inplace.c +* in-place running sum +* +* Alain de Cheveigné, CNRS/Ircam Jun 2002 +* (c) 2002 CNRS +*/ + +/* Replaces each sample of the input matrix by the sum of itself and its N-1 +* neighbors. The operation is done in-place on the input argument. The last +* N-1 samples of each row are invalid. */ + +#include <math.h> +#include "mex.h" + +/* #define MWRAP */ +#include "mwrap_check.h" + +/* Input Arguments */ +#define X_IN prhs[0] +#define N_IN prhs[1] + +/* Output Arguments */ + +static void rsmth( + double *xp, /* matrix to sum */ + double N, /* window size */ + int m, /* rows */ + int n /* columns */ + ) +{ + int j,k, Ni; + double Nr, tmp, sum, *x, *xx, *bump; + + Ni = (int) floor(N); + Nr = N - (double) Ni; + + + if (Nr == 0.0) { + /* N is integer: simple summation */ + for (j=0;j<n;j++) { /* for each column */ + /* prefill running sum */ + sum=0.0; + x=&xp[j*m]; + bump=x+Ni-1; + while (x<bump) { /* for each row */ + sum += GET(*x); + x++; + } + /* N */ + x=&xp[j*m]; + bump=x+m-(Ni-1); + xx=x+Ni-1; + while(x<bump) { + tmp = GET(*x); + sum += GET(*xx); + SET(*x) = sum; + sum -= tmp; + x++; xx++; + } + } + } else { + /* N is fractionary: interpolate */ + for (j=0;j<n;j++) { /* for each column */ + /* prefill running sum */ + sum=0.0; + x=&xp[j*m]; + bump=x+Ni-1; + while (x<bump) { /* for each row */ + sum += GET(*x); + x++; + } + /* N */ + x=&xp[j*m]; + bump=x+m-(Ni-1); + xx=x+Ni-1; + while(x<bump) { + tmp = GET(*x); + sum += GET(*xx); + SET(*x) = sum + Nr * GET(*xx); + sum -= tmp; + x++; xx++; + } + } + } + return; +} + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[] + ) +{ + double *xp, *Np, N; + int n, m; + + /* Check for proper number of arguments */ + if (nrhs != 2) { + mexErrMsgTxt("RSUM_INPLACE takes 2 input arguments"); + } + + /* Check type of input */ + if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || + mxIsSparse(X_IN) || !mxIsDouble(X_IN) ) { + mexErrMsgTxt("RSUM_INPLACE: X should be doubles"); + } + if (!mxIsNumeric(N_IN) || mxIsComplex(N_IN) || + mxIsSparse(N_IN) || !mxIsDouble(N_IN) ) { + mexErrMsgTxt("RSUM_INPLACE: N should be doubles"); + } + m=mxGetM(X_IN); /* rows */ + n=mxGetN(X_IN); /* columns */ + + if (mxGetM(N_IN)*mxGetN(N_IN) != 1) { + mexErrMsgTxt("RSUM_INPLACE: N should be a scalar"); + } + + xp = mxGetPr(X_IN); + Np = mxGetPr(N_IN); + N = *Np; + + if (m < ceil(N)) { + mexErrMsgTxt("RSUM_INPLACE: X should have at least N rows"); + } + + checkin_matrix((mxArray *) X_IN); + + /* Do the actual computations in a subroutine */ + rsmth(xp, N, m, n); + + checkout_matrix((mxArray *) X_IN); + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/src/rsum_inplace.prefix Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +/* * Prefix file for the Metrowerks CodeWarrior compiler. * There are no command line options. */ /* $Revision: 1.2 $ */ #define MATLAB_MEX_FILE #ifndef rez /* Both Rez and the C Compiler use this file. Put any C specific * code within ifndef rez..endif clauses, such as here. */ #endif #define NDEBUG \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/yink.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,147 @@ +function r=yink(p,fileinfo) +% YINK - fundamental frequency estimator +% new version (feb 2003) +% +% + +%global jj; +%jj=0; +% process signal a chunk at a time +idx=p.range(1)-1; +totalhops=round((p.range(2)-p.range(1)+1) / p.hop); +r1=nan*zeros(1,totalhops);r2=nan*zeros(1,totalhops); +r3=nan*zeros(1,totalhops);r4=nan*zeros(1,totalhops); +idx2=0+round(p.wsize/2/p.hop); +while (1) + start = idx+1; + stop = idx+p.bufsize; + stop=min(stop, p.range(2)); + xx=sf_wave(fileinfo, [start, stop], []); +% if size(xx,1) == 1; xx=xx'; end + xx=xx(:,1); % first channel if multichannel + [prd,ap0,ap,pwr]=yin_helper(xx,p); + n=size(prd ,2); + if (~n) break; end; + idx=idx+n*p.hop; + + r1(idx2+1:idx2+n)= prd; + r2(idx2+1:idx2+n)= ap0; + r3(idx2+1:idx2+n)= ap; + r4(idx2+1:idx2+n)= pwr; + idx2=idx2+n; +end +r.r1=r1; % period estimate +r.r2=r2; % gross aperiodicity measure +r.r3=r3; % fine aperiodicity measure +r.r4=r4; % power +sf_cleanup(fileinfo); +% end of program + + + +% Estimate F0 of a chunk of signal +function [prd,ap0,ap,pwr]=yin_helper(x,p,dd) +smooth=ceil(p.sr/p.lpf); +x=rsmooth(x,smooth); % light low-pass smoothing +x=x(smooth:end-smooth+1); +[m,n]=size(x); +maxlag = ceil(p.sr/p.minf0); +minlag = floor(p.sr/p.maxf0); +mxlg = maxlag+2; % +2 to improve interp near maxlag + + +hops=floor((m-mxlg-p.wsize)/p.hop); +prd=zeros(1,hops); +ap0=zeros(1,hops); +ap=zeros(1,hops); +pwr=zeros(1,hops); +if hops<1; return; end + +% difference function matrix +dd=zeros(floor((m-mxlg-p.hop)/p.hop),mxlg); +if p.shift == 0 % windows shift both ways + lags1=round(mxlg/2) + round((0:mxlg-1)/2); + lags2=round(mxlg/2) - round((1:mxlg)/2); + lags=[lags1; lags2]; +elseif p.shift == 1 % one window fixed, other shifts right + lags=[zeros(1,mxlg); 1:mxlg]; +elseif p.shift == -1 % one window fixed, other shifts right + lags=[mxlg-1:-1:0; mxlg*ones(1,mxlg)]; +else + error (['unexpected shift flag: ', num2str(p.shift)]); +end +rdiff_inplace(x,x,dd,lags,p.hop); +rsum_inplace(dd,round(p.wsize/p.hop)); +dd=dd'; +[dd,ddx]=minparabolic(dd); % parabolic interpolation near min +cumnorm_inplace(dd);; % cumulative mean-normalize + +% first period estimate +%global jj; +for j=1:hops + d=dd(:,j); + if p.relflag + pd=dftoperiod2(d,[minlag,maxlag],p.thresh); + else + pd=dftoperiod(d,[minlag,maxlag],p.thresh); + end + ap0(j)=d(pd+1); + prd(j)=pd; +end + +% replace each estimate by best estimate in range +range = 2*round(maxlag/p.hop); +if hops>1; prd=prd(mininrange(ap0,range*ones(1,hops))); end +%prd=prd(mininrange(ap0,prd)); + + +% refine estimate by constraining search to vicinity of best local estimate +margin1=0.6; +margin2=1.8; +for j=1:hops + d=dd(:,j); + dx=ddx(:,j); + pd=prd(j); + lo=floor(pd*margin1); lo=max(minlag,lo); + hi=ceil(pd*margin2); hi=min(maxlag,hi); + pd=dftoperiod(d,[lo,hi],0); + ap0(j)=d(pd+1); + pd=pd+dx(pd+1)+1; % fine tune based on parabolic interpolation + prd(j)=pd; + + % power estimates + idx=(j-1)*p.hop; + k=(1:ceil(pd))'; + x1=x(k+idx); + x2=k+idx+pd-1; + interp_inplace(x,x2); + x3=x2-x1; + x4=x2+x1; + x1=x1.^2; rsum_inplace(x1,pd); + x3=x3.^2; rsum_inplace(x3,pd); + x4=x4.^2; rsum_inplace(x4,pd); + + x1=x1(1)/pd; + x2=x2(1)/pd; + x3=x3(1)/pd; + x4=x4(1)/pd; + % total power + pwr(j)=x1; + + % fine aperiodicity + ap(j)=(eps+x3)/(eps+(x3+x4)); % accurate, only for valid min + + %ap(j) + %plot(min(1, d)); pause + + prd(j)=pd; +end + + + +%cumulative mean-normalize +function y=cumnorm(x) +[m,n]=size(x); +y = cumsum(x); +y = (y)./ (eps+repmat((1:m)',1,n)); % cumulative mean +y = (eps+x) ./ (eps+y);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/private/yink.m~ Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,147 @@ +function r=yink(p,fileinfo) +% YINK - fundamental frequency estimator +% new version (feb 2003) +% +% + +%global jj; +%jj=0; +% process signal a chunk at a time +idx=p.range(1)-1; +totalhops=round((p.range(2)-p.range(1)+1) / p.hop); +r1=nan*zeros(1,totalhops);r2=nan*zeros(1,totalhops); +r3=nan*zeros(1,totalhops);r4=nan*zeros(1,totalhops); +idx2=0+round(p.wsize/2/p.hop); +while (1) + start = idx+1; + stop = idx+p.bufsize; + stop=min(stop, p.range(2)); + xx=sf_wave(fileinfo, [start, stop], []); +% if size(xx,1) == 1; xx=xx'; end + xx=xx(:,1); % first channel if multichannel + [prd,ap0,ap,pwr]=yin_helper(xx,p); + n=size(prd ,2); + if (~n) break; end; + idx=idx+n*p.hop; + + r1(idx2+1:idx2+n)= prd; + r2(idx2+1:idx2+n)= ap0; + r3(idx2+1:idx2+n)= ap; + r4(idx2+1:idx2+n)= pwr; + idx2=idx2+n; +end +r.r1=r1; % period estimate +r.r2=r2; % gross aperiodicity measure +r.r3=r3; % fine aperiodicity measure +r.r4=r4; % power +sf_cleanup(fileinfo); +% end of program + + + +% Estimate F0 of a chunk of signal +function [prd,ap0,ap,pwr]=yin_helper(x,p,dd) +smooth=ceil(p.sr/p.lpf); +x=rsmooth(x,smooth); % light low-pass smoothing +x=x(smooth:end-smooth+1); +[m,n]=size(x); +maxlag = ceil(p.sr/p.minf0); +minlag = floor(p.sr/p.maxf0); +mxlg = maxlag+2; % +2 to improve interp near maxlag + + +hops=floor((m-mxlg-p.wsize)/p.hop); +prd=zeros(1,hops); +ap0=zeros(1,hops); +ap=zeros(1,hops); +pwr=zeros(1,hops); +if hops<1; return; end + +% difference function matrix +dd=zeros(floor((m-mxlg-p.hop)/p.hop),mxlg); +if p.shift == 0 % windows shift both ways + lags1=round(mxlg/2) + round((0:mxlg-1)/2); + lags2=round(mxlg/2) - round((1:mxlg)/2); + lags=[lags1; lags2]; +elseif p.shift == 1 % one window fixed, other shifts right + lags=[zeros(1,mxlg); 1:mxlg]; +elseif p.shift == -1 % one window fixed, other shifts right + lags=[mxlg-1:-1:0; mxlg*ones(1,mxlg)]; +else + error (['unexpected shift flag: ', num2str(p.shift)]); +end +rdiff_inplace(x,x,dd,lags,p.hop); +rsum_inplace(dd,round(p.wsize/p.hop)); +dd=dd'; +[dd,ddx]=minparabolic(dd); % parabolic interpolation near min +cumnorm_inplace(dd);; % cumulative mean-normalize + +% first period estimate +%global jj; +for j=1:hops + d=dd(:,j); + if p.relflag + pd=dftoperiod(d,[minlag,maxlag],p.thresh); + else + pd=dftoperiod2(d,[minlag,maxlag],p.thresh); + end + ap0(j)=d(pd+1); + prd(j)=pd; +end + +% replace each estimate by best estimate in range +range = 2*round(maxlag/p.hop); +if hops>1; prd=prd(mininrange(ap0,range*ones(1,hops))); end +%prd=prd(mininrange(ap0,prd)); + + +% refine estimate by constraining search to vicinity of best local estimate +margin1=0.6; +margin2=1.8; +for j=1:hops + d=dd(:,j); + dx=ddx(:,j); + pd=prd(j); + lo=floor(pd*margin1); lo=max(minlag,lo); + hi=ceil(pd*margin2); hi=min(maxlag,hi); + pd=dftoperiod(d,[lo,hi],0); + ap0(j)=d(pd+1); + pd=pd+dx(pd+1)+1; % fine tune based on parabolic interpolation + prd(j)=pd; + + % power estimates + idx=(j-1)*p.hop; + k=(1:ceil(pd))'; + x1=x(k+idx); + x2=k+idx+pd-1; + interp_inplace(x,x2); + x3=x2-x1; + x4=x2+x1; + x1=x1.^2; rsum_inplace(x1,pd); + x3=x3.^2; rsum_inplace(x3,pd); + x4=x4.^2; rsum_inplace(x4,pd); + + x1=x1(1)/pd; + x2=x2(1)/pd; + x3=x3(1)/pd; + x4=x4(1)/pd; + % total power + pwr(j)=x1; + + % fine aperiodicity + ap(j)=(eps+x3)/(eps+(x3+x4)); % accurate, only for valid min + + %ap(j) + %plot(min(1, d)); pause + + prd(j)=pd; +end + + + +%cumulative mean-normalize +function y=cumnorm(x) +[m,n]=size(x); +y = cumsum(x); +y = (y)./ (eps+repmat((1:m)',1,n)); % cumulative mean +y = (eps+x) ./ (eps+y);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/testyin.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,3 @@ +% test YIN + +yin 'female04_16K.wav'
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Descriptors/yin/yin/yin.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,1 @@ +function r=yin(x,p) % YIN - Fundamental frequency (F0) of file or vector % % YIN('FILE') estimates and plots the F0 of signal in FILE. % F0 is plotted in octaves re: 440 Hz, together with residual % (aperiodicity) and power measures. A 'best' estimate is printed. % YIN(X,SR) estimates and plots the F0 of array X assuming sampling rate SR (Hz). % % R=YIN(X) produces no plot but returns result in structure R: % R.f0: fundamental frequency in octaves re: 440 Hz % R.ap: aperiodicity measure (ratio of aperiodic to total power) % R.pwr: period-smoothed instantaneous power % (R also lists the parameters used by YIN) % % YIN(NAME,P) uses parameters stored in P: % P.minf0: Hz - minimum expected F0 (default: 30 Hz) % P.maxf0: Hz - maximum expected F0 (default: SR/(4*dsratio)) % P.thresh: threshold (default: 0.1) % P.relfag: if ~0, thresh is relative to min of difference function (default: 1) % P.hop: samples - interval between estimates (default: 32) % P.range: samples - range of samples ([start stop]) to process % P.bufsize: samples - size of computation buffer (default: 10000) % P.sr: Hz - sampling rate (usually taken from file header) % P.wsize: samples - integration window size (defaut: SR/minf0) % P.lpf: Hz - intial low-pass filtering (default: SR/4) % P.shift 0: shift symmetric, 1: shift right, -1: shift left (default: 0) % % See 'yin.html' for more info. % Version 28 July 2003. % Alain de Cheveign? CNRS/Ircam, 2002. % Copyright (c) 2002 Centre National de la Recherche Scientifique. % % Permission to use, copy, modify, and distribute this software without % fee is hereby granted FOR RESEARCH PURPOSES only, provided that this % copyright notice appears in all copies and in all supporting % documentation, and that the software is not redistributed for any % fee (except for a nominal shipping charge). % % For any other uses of this software, in original or modified form, % including but not limited to consulting, production or distribution % in whole or in part, specific prior permission must be obtained from CNRS. % Algorithms implemented by this software may be claimed by patents owned % by CNRS, France Telecom, Ircam or others. % % The CNRS makes no representations about the suitability of this % software for any purpose. It is provided "as is" without express % or implied warranty. % Hidden parameters are integration window size (set equal to sr/minf0), search % range for 'best neighboring estimate' (set equal to +/- sr/(2*minf0)), maximum % expected width of period dip (stop/start ratio == 1.85), margin for "beam search" % of final estimate (+1.8/-0.6 times the initial estimate). % default parameter values ([]: to be determined) minf0 = 30; % Hz - minimum frequency maxf0 = []; % Hz - maximum frequency wsize = []; % s - integration window size lpf = []; % Hz - lowpass prefiltering cutoff thresh = 0.1; % difference function threshold relflag = 1; % if true threshold is relative to global min of difference function bufsize=10000; % computation buffer size hop = 32; % samples - interval between estimates range=[]; % range of file samples to process sr=[]; % sampling rate shift=0; % flag to control the temporal shift of analysis windows (left/sym/right) plotthreshold=0.2; % aperiodicity above which plot is green or yellow % if 2~=exist('allread') % error('sf routines missing: put them in your path & try again'); % end % handle parameters if nargin<1; help yin; return; end if nargin<2; p=[]; end fileinfo=sf_info(x); if ~isempty(fileinfo.sr) p.sr=fileinfo.sr; end % get sr from file if fileinfo.nchans > 1 disp(['warning: using column 1 of ', num2str(fileinfo.nchans), '-column data']); end if isa(p, 'double') p.sr=p; end if ~isfield(p, 'sr'); p.sr=sr; end if isempty(p.sr); error('YIN2: must specify SR'); end if ~isfield(p, 'range') | isempty(p.range); p.range=[1 fileinfo.nsamples]; end if ~isfield(p, 'minf0'); p.minf0=minf0; end if ~isfield(p, 'thresh'); p.thresh=thresh; end if ~isfield(p, 'relflag'); p.relflag=relflag; end if ~isfield(p, 'bufsize'); p.bufsize=bufsize; end if ~isfield(p, 'hop'); p.hop=hop; end if ~isfield(p, 'maxf0'); p.maxf0=floor(p.sr/4); end % default if ~isfield(p, 'wsize'); p.wsize=ceil(p.sr/p.minf0); end % default if ~isfield(p, 'lpf'); p.lpf=p.sr/4; end % default if mod(p.hop,1); error('hop should be integer'); end if ~isfield(p, 'shift'); p.shift=shift; end % default if ~isfield(p, 'plotthreshold'); p.plotthreshold=plotthreshold; end % default % estimate period r=yink(p,fileinfo); prd=r.r1; % period in samples ap0=r.r2; % gross aperiodicity measure ap= r.r3; % fine aperiodicity measure pwr=r.r4; % period-smoothed instantaneous power f0 = log2(p.sr ./ prd) - log2(440); % convert to octaves re: 440 Hz figure(); plot(prd); dlmwrite('pitch.txt',prd,'delimiter','\n'); display(length(prd)); % load estimates and major parameters in result structure clear r; r.f0 = f0; r.ap0 = ap0; r.ap = ap; r.pwr = pwr; r.sr = p.sr; r.range=p.range; r.minf0 = p.minf0; r.maxf0 = p.maxf0; r.thresh=p.thresh; r.relflag=p.relflag; r.hop = p.hop; r.bufsize = p.bufsize; r.wsize = p.wsize; r.lpf = p.lpf; r.shift = p.shift; r.plotthreshold=p.plotthreshold; % plot estimates (if nargout == 0) if nargout<1 if isnan(f0) display('no estimates: signal too short or buffersize too small'); return; end % choose sample to report as "the" f0 of the entire signal [mn, idx] = min(ap0); best=f0(idx); disp(['best: ', num2str(2^best*440), 'Hz (', note(best),... ') at ', num2str(idx/(p.sr/p.hop)), 's (aperiodic/total power: ', num2str(mn), ')']); % plot f0 in 3 colors according to periodicity good = f0; good(find(ap0>p.plotthreshold*2)) = nan; best = f0; best(find(ap0>p.plotthreshold)) = nan; subplot(211); fsr=p.sr/p.hop; nframes=size(prd,2); if nframes <2; error('F0 track is too short to plot'); end plot((1:nframes)/fsr, f0, 'y', (1:nframes)/fsr, good, 'g', (1:nframes)/fsr, best, 'b'); lo = max(min(f0),min(good)); hi=min(max(f0),max(good)); set(gca, 'ylim', [lo-0.5; hi+0.5]); set(gca, 'xlim', [1,nframes]/fsr); set(get(gca,'ylabel'), 'string', 'Oct. re: 440 Hz'); set(get(gca,'xlabel'), 'string', 's'); % plot periodicity ap0=max(0,ap0); ap=max(0,ap); ap1=ap; ap(find(ap0>plotthreshold)) = nan; subplot(413); plot((1:nframes), ap0.^0.5, 'b'); %subplot(413); plot((1:nframes), ap0.^0.5, 'c', (1:nframes), ap1.^0.5, 'y', (1:nframes), ap.^0.5, 'b'); set(gca, 'ylim', [0 1]); set(gca, 'xlim', [1, nframes]); set(get(gca,'ylabel'), 'string', 'sqrt(apty)'); % plot power subplot(414); plot((1:nframes), sqrt(pwr), 'b'); set(gca, 'xlim', [1, nframes]); set(get(gca,'ylabel'), 'string', 'sqrt(power)'); set(get(gca,'xlabel'), 'string', 'frames'); if isa(x, 'double') set(gcf, 'Name', 'workspace matrix'); else set (gcf, 'Name', x); end end % convert octave re 440 to note: function s=note(o) n=round(12*o); cents = 100*(12*o-n); oct=floor((n-3)/12)+5; chroma=mod(n,12); chromalist = {'A'; 'A#'; 'B'; 'C'; 'C#'; 'D'; 'D#'; 'E'; 'F'; 'F#';... 'G'; 'G#'}; cents = sprintf('%+.0f',cents); s=[char(chromalist(chroma+1)),num2str(oct),' ',num2str(cents), ' cents']; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/General/getVariables.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,13 @@ +% speech range = around 80Hz to 260Hz +% singing range extend up to 1050 Hz +function [ out ] = getVariables( name ) + +switch name + case 'getMinFreq' + out = 100; %100Hz is the default minimum fundatmental frequency + + case 'getMaxFreq' + out = 1050; %1050 Hz is the default maximum fundamental frequency + otherwise + out = 'EEE'; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/General/openFile.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,10 @@ +function [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ) + + [ x, fs ] = wavread(['../' currentSampleName]); + frameLength = 2^(ceil(log2(fs/100))); + noOfFrames = floor( length(x) / frameLength); %skip the final samples if necessary + + % Convert stereo to mono + if (size(x, 2)==2) + x = mean(x')'; + end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/ProcessDatabase.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,467 @@ +function [] = ProcessDatabase( databaseName, functionName, sampleStartNum, OVERWRITE ) +% run from D:\Dropbox\BUPTResearch2011\EmotionDetectionCode\Code +machineName = getenv('COMPUTERNAME'); + +switch machineName + case 'LAPTOP' + % for Dawn's laptop + root = 'D:\Dropbox\' + case 'SLATE1' + % for Dawn's Slate + root = 'E:\Dropbox\' + case 'DAWNBLACK-PC' + % for Dawn's work PC + root = 'C:\Users\dawn\Dropbox\' +end + +addpath ([root 'BUPTResearch2011\emotionDetectionCode\Code']) +addpath ([root 'BUPTResearch2011\emotionDetectionCode\Code\Descriptors\Matlab\MPEG7']) +addpath ([root 'BUPTResearch2011\emotionDetectionCode\Code\Descriptors\Matlab\Common']) +addpath ([root 'BUPTResearch2011\emotionDetectionCode\Code\Descriptors\Matlab\MPEG7\FromWeb']) +addpath ([root 'BUPTResearch2011\emotionDetectionCode\Code\Descriptors\PRAAT']) +addpath ([root 'BUPTResearch2011\emotionDetectionCode\Code\Collation']) +addpath ([root 'BUPTResearch2011\emotionDetectionCode\Code\General']) + +switch databaseName + case 'ChineseOpera' + cd([root 'BUPTResearch2011/Data/Opera/TestDatabase']) + startEmotion = 4; + case 'MandarinSpeech' + cd([root 'BUPTResearch2011/Data/Database/EditedRecording']) + startEmotion = 3; + case 'SpeechTestFiles' + cd([root 'BUPTResearch2011\Data\SpeechTestFiles\MyAnnotatedFiles']) + startEmotion = 3; +end + +fileStructure = dir; + +% how many samples do we have? +noOfSamples = size( fileStructure ); +firstfileOpen = 1; + +if( sampleStartNum < 3 ) + sampleStartNum = 3; +end + +for sampleNum = sampleStartNum : noOfSamples(1) + sampleDirName = fileStructure(sampleNum).name + if( fileStructure(sampleNum).isdir == 1 ) % only directories + cd( sampleDirName ); + % how many emotions for that sample? + validEmotionList = dir; + noOfEmotions = size( validEmotionList ); + + for emotionNum = startEmotion : noOfEmotions(1) + emotionName = validEmotionList(emotionNum).name; + cd( emotionName ); + % how many samples for that emotion? + sampleNames = dir; + noOfSamples = size( sampleNames ); + + for sampleNumber = 3 : noOfSamples(1) + if( sampleNames(sampleNumber).isdir == 0 ) % skip directories + currentSampleName = sampleNames(sampleNumber).name + % get the file type + extName = currentSampleName(length( currentSampleName ) - 3: end); + % is it a .wav file? + if strcmp( extName, '.wav' ) + newDirName = currentSampleName( 1 : length( currentSampleName ) - 4 ); + % if none exists, make a new directory for all the results of + % our calculations. + mkdir( [ newDirName '_metrics'] ); + cd ( [ newDirName '_metrics'] ); + + % call your function here + switch functionName + % functions that start with 'get' simply read + % from existing text files + % functions that start with 'detect' check if + % the results file already exists, if they do + % not (or the OVERWRITE flag is true) then the + % results are calculated afresh. + + % case 'getsampleDependantThresholds' + % statsFileName = '../../../../FeatureSets/singingsampleDependantThresholds.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % % is this the first sample in this emotion for + % % this sample? + % if( sampleNumber == 3 ) + % if ( firstfileOpen == 0 ) + % if(strfind( oldDirName , 'fem') > 0 ) + % sampleName = oldDirName( 6 : end ); + % meanET = mean( energyThreshold ); + % meanSCT = mean( spectralCentroidThreshold ); + % fprintf( statsFileID, '%s \t %f \t %f \n', sampleName, meanET, meanSCT ); + % elseif(strfind( oldDirName , 'male') > 0 ) + % sampleName = oldDirName( 6 : end ); + % meanET = mean( energyThreshold ); + % meanSCT = mean( spectralCentroidThreshold ); + % fprintf( statsFileID, '%s \t %f \t %f \n', sampleName, meanET, meanSCT ); + % end + % end + % + % spectralCentroidThreshold = []; + % energyThreshold = []; + % firstfileOpen = 0; + % end + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % [st , et] = getsampleDependantThresholds( newDirName, x, fs, statsFileID, frameLength, noOfFrames ); + % spectralCentroidThreshold = [spectralCentroidThreshold st]; + % energyThreshold = [energyThreshold et]; + % oldDirName = newDirName; + % + % fclose( statsFileID ); + % case 'getPitchStatistics' + % statsHeaderFileName = '../../../../FeatureSets/singingHeaderPitchStats.txt'; + % statsHeaderFileID = fopen( statsHeaderFileName, 'w' ); + % if( firstfileOpen == 1 ) + % fprintf( statsHeaderFileID, 'name \t unvoiced to voiced frame ratio \t mean voiced pitch \t median voiced pitch \t stdev voiced pitch \t variance voiced pitch \t min voiced pitch \t max voiced pitch \t range voiced pitch \t mean voiced pitch grad \t median voiced pitch grad \t stdev pitch grad \t variance voiced pitch grad \t min voiced pitch grad \t max voiced pitch grad \t range voiced pitch grad \n'); + % firstfileOpen = 0; + % end + % fclose(statsHeaderFileID); + % + % statsFileName = '../../../../FeatureSets/singingPitchStats.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getPitchStatistics( newDirName, x, fs, statsFileID, frameLength, noOfFrames ); + % fclose( statsFileID ); + % + % + % + % + % + % %%%%%%%%%%%%%%% careful %%%%%%%%%%%%%%% + % case 'getHNRStatistics' + % statsHeaderFileName = '../../../../FeatureSets/singingHeaderHNRStats.txt'; + % statsHeaderFileID = fopen( statsHeaderFileName, 'w' ); + % if( firstfileOpen == 1 ) + % fprintf( statsHeaderFileID, 'name \t \t \t \t \t \n'); + % firstfileOpen = 0; + % end + % fclose(statsHeaderFileID); + % + % statsFileName = '../../../../FeatureSets/singingHNRStats.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getHNRStatistics( newDirName, x, fs, statsFileID, frameLength, noOfFrames ); + % fclose( statsFileID ); + % + % + % + % + % case 'getEnergyStatistics' + % statsFileName = '../../../../FeatureSets/singingEnergyStats.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % if( firstfileOpen == 1 ) + % fprintf( statsFileID, 'name \t mean energy \t mean voiced energy \t mean unvoiced energy \n'); + % firstfileOpen = 0; + % end + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getEnergyStatistics( newDirName, x, fs, statsFileID, frameLength, noOfFrames ); + % fclose( statsFileID ); + % + % + % % %%%%% for energy-related features %%%%% + % % case 'getOtherEnergyStatistics' + % % statsFileName = '../../../../FeatureSets/singingOtherEnergyStats.txt'; + % % statsFileID = fopen( statsFileName, 'a' ); + % % if( firstfileOpen == 1 ) + % % fprintf( statsFileID, 'name \t derivation of non silent energy \t derivation of voiced energy \t derivation of unvoiced energy \n'); + % % firstfileOpen = 0; + % % end + % % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % % getOtherEnergyStatistics( newDirName, x, fs, statsFileID, frameLength, noOfFrames ); + % % fclose( statsFileID ); + % + % + % + % case 'getBasicDescriptors' + % % includes audio waveform and audio power + % statsFileName = '../../../../FeatureSets/singingbasicDescriptors.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % if( firstfileOpen == 1 ) + % fprintf( statsFileID, 'name \t meanAWF \t STDAWF \t minAWF \t maxAWF \t rangeAWF \t meanAP \t STDAP \t minAP \t maxAP \t rangeAP \t meanHR \t STDHR \t minHR \t maxHR \t rangeHR \n'); + % firstfileOpen = 0; + % end + % + % getBasicDescriptors( newDirName, statsFileID ); + % fclose( statsFileID ); + % case 'getBasicSpectralDescriptors' + % % includes audio spectrum envelope, audio spectrum + % % centroid, audio spectrum spread, audio spectrum flatness + % statsFileName = '../../../../FeatureSets/singingbasicSpectralDescriptors.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % % Add a header column + % if( firstfileOpen == 1 ) + % fprintf( statsFileID, ' name \t meanASE \t STDASE \t minASE \t maxASE \t rangeASE \t meanASC \t STDASC \t minASC \t maxASC \t rangeASC \t meanASS \t STDASS \t minASS \t maxASS \t rangeASS \t meanASF \t STDASF \t minASF \t maxASF \t rangeASF \n'); + % % fprintf( statsFileID, ' name \t meanAWE \t STDAWE \t minAWE \t maxAWE \t rangeAWE \t meanAP \t STDAP \t minAP \t maxAP \t rangeAP \n'); + % firstfileOpen = 0; + % end + % + % getBasicSpectralDescriptors1( newDirName, statsFileID ); + % fclose( statsFileID ); + % disp('NOT FINISHED'); + % case 'getTemporalTimbralDescriptors' + % + % % includes log attack time and temporal + % % centroid + % statsFileName = '../../../../FeatureSets/singingtemporalTimbralDescriptors.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % % Add a header column + % if( firstfileOpen == 1) + % fprintf( statsFileID, ' name \t temporal centriod \t log attack time \n'); + % firstfileOpen = 0; + % end + % + % getTemporalTimbralDescriptors( newDirName, statsFileID ); + % fclose( statsFileID ); + % case 'getSpectralTimbralDescriptors' + % % collate the Harmonic Spectral Centroid, + % % Harmonmic Spectral Deviation, Harmonic + % % Spectral Spread, Harmonic Spectral Variation + % % and Spectral Centroid + % + % statsFileName = '../../../../FeatureSets/singingspectralTimbralDescriptors.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % % Add a header column + % if( firstfileOpen == 1) + % fprintf( statsFileID, ' name \t HSC \t HSD \t HSS \t HSV \t SC \n'); + % firstfileOpen = 0; + % end + % + % getSpectralTimbralDescriptors( newDirName, statsFileID ); + % fclose( statsFileID ); + % case 'getMFCCs' + % if( firstfileOpen ) + % statsFileName = '../../../../FeatureSets/singingMFCCDescriptorsHeader.txt'; + % statsFileID = fopen( statsFileName, 'w' ); + % + % % Add a header file + % % putMFCCHeader( statsFileID ); + % fclose( statsFileID ); + % end + % + % statsFileName = '../../../../FeatureSets/singingMFCCDescriptors.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % + % getMFCCStatistics( newDirName, statsFileID ); + % fclose( statsFileID ); + % + % %lin + case 'getPRAAT' + % collate all Shimmer values + statsFileName = '../../../../../../../TestResults/singingPRAATStats.txt'; + statsFileID = fopen( statsFileName, 'w' ); + % Add a header row + if( firstfileOpen == 1 ) + fprintf( statsFileID, 'Metrics calculated using the PRAAT software. \n' ); + fprintf( statsFileID, ' name \t JITTER: ddp \t local \t ppq5 \t rap \t SHIMMER: local \t dda \t apq3 \t apq5 \t apq11 \t FORMANTS: number of formants \t positions \t bandwidths \n'); + firstfileOpen = 0; + end + + get_PRAAT( newDirName, statsFileID ); + fclose( statsFileID ); + + % case 'getJitter' + % % collate Jitter + % statsFileName = '../../../../FeatureSets/singingJitterStats.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % + % getJitter( newDirName, statsFileID ); + % fclose( statsFileID ); + % + % + % + % %%%%% get all the pitch and energy related features + + + % case 'getFeaturesOfDuration' + % % collate number, mean and ratio duration of + % % unvoiced and voiced sounds, median and + % % standard deviation number of voiced sounds + % statsFileName = '../../../../FeatureSets/singingDurationStats.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % getFeaturesOfDuration( newDirName, statsFileID ); + % fclose( statsFileID ); + % + % case 'getFeaturesOfRelativePitch' + % % collate relative pitch maximum, minimum and + % % the position of those pitches + % statsFileName = '../../../../FeatureSets/singingRelativePitchStats.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % getFeaturesOfRelativePitch( newDirName, statsFileID ); + % fclose( statsFileID ); + % + % case 'getFeaturesOfRelativeEnergyGradient' + % % collate the value and position of relative + % % maximum of energy gradient + % statsFileName = '../../../../FeatureSets/singingRelativeEnergyGradientStats.txt'; + % statsFileID = fopen( statsFileName, 'a' ); + % getFeaturesOfRelativeEnergyGradient( newDirName, statsFileID ); + % fclose( statsFileID ); + % + % + % + case 'calculatePitch' + detect_pitch( currentSampleName, OVERWRITE ); + + case 'calculateSilence' + detect_silence( currentSampleName, OVERWRITE ); + + case 'calculateVUV' + % find the voiced and unvoiced frames using + % the harmonic ratio and pitch information + detect_VoicedUnvoiced( currentSampleName, OVERWRITE ); + + % case 'calculateAudioWaveform' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % AW_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % case 'calculateAudioSpectrumEnvelope' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % ASE_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % case 'calculateAudioSpectrumCentriodAndSpread' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % ASC_ASS_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % case 'calculateAudioSpectrumBasisAndProjection' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % ASB_ASP_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + case 'calculateAudioPower' + +% [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + detect_AudioPower( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % case 'calculateLogAttackTIme' + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % LAT_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % case 'calculateTemporalCentroid' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % TC_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % case 'calculateHarmonicRatio' + % + % % can't find this function yet + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % HR_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + case 'calculateHarmonicNoiseRatio' + detect_HNR( currentSampleName, OVERWRITE ); + % + % case 'calculateSpectralTimbralDescriptors' + % + % % these are best calculated together. + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % spectralTimbralDescriptors( currentSampleName, x, fs, frameLength, noOfFrames ); + + case 'calculateAudioSpectrumFlatness' + + [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + AudioSpectralFlatness_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + + % case 'calculateMFCCs' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % MFCC_Detection( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % case 'calculatesampleDependantThresholds' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % sampleDependantThresholds( currentSampleName, x, fs, frameLength, noOfFrames ); + % + % + % %%%%% pitch related features %%%%% + % case 'calculatesampleVoicedSoundDuration' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getVoicedSoundDuration( currentSampleName, x, fs, noOfFrames, frameLength ); + % + % case 'calculatesampleRelativePitch' + % + % singingpitchStatsFileName = 'singingpitchStats.txt'; + % singingpitchStatsFileID = fopen( singingpitchStatsFileName, 'r' ); + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getRelativePitch( singingpitchStatsFileID, currentSampleName, x, fs, noOfFrames, frameLength, 7 ); + % fclose(singingpitchStatsFileID); + % + % case 'calculatesampleRelativePitchFeatures' + % + % singingpitchStatsFileName1 = '../../../../FeatureSets/singingpitchStats.txt'; + % singingpitchStatsFileID1 = fopen( singingpitchStatsFileName1, 'r' ); + % singingpitchStatsFileName2 = 'singingpitchStats.txt'; + % singingpitchStatsFileID2 = fopen( singingpitchStatsFileName2, 'r' ); + % [averageValues1 averageValues2 averageValues3] = getRelativePitchFeatures( singingpitchStatsFileID1, 7 ); + % getEachRelativePitchFeature( singingpitchStatsFileID2, currentSampleName, 7, averageValues1, averageValues2, averageValues3 ); + % fclose(singingpitchStatsFileID1); + % fclose(singingpitchStatsFileID2); + % + % case 'calculatesamplePitchPosition' + % + % singingpitchStatsFileName = 'singingpitchStats.txt'; + % singingpitchStatsFileID = fopen( singingpitchStatsFileName, 'r' ); + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getPitchPosition( singingpitchStatsFileID, currentSampleName, x, fs, noOfFrames, frameLength, 7 ); + % fclose(singingpitchStatsFileID); + % + % + % + % %%%%% energy related features %%%%% + % case 'calculatesampleRelativeEnergyByMeanAll' + % + % singingenergyStatsFileName = '../../../../FeatureSets/singingenergyStats.txt'; + % singingenergyhStatsFileID = fopen( singingenergyStatsFileName, 'r' ); + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getRelativeEnergyByMeanAll( singingenergyhStatsFileID, currentSampleName, x, fs, noOfFrames, frameLength, 3 ); + % fclose(singingenergyhStatsFileID); + % + % case 'calculatesampleRelativeEnergyByMeanVoiced' + % + % singingenergyStatsFileName = '../../../../FeatureSets/singingenergyStats.txt'; + % singingenergyhStatsFileID = fopen( singingenergyStatsFileName, 'r' ); + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getRelativeEnergy( singingenergyhStatsFileID, currentSampleName, x, fs, noOfFrames, frameLength, 3 ); + % fclose(singingenergyhStatsFileID); + % + % case 'calculatesampleRelativeEnergyGradientByMeanAll' + % + % statsFileName = '../../../../FeatureSets/singingenergyStats.txt'; + % statsFileID = fopen( statsFileName, 'r' ); + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getRelativeEnergyGradient( statsFileID, newDirName, x, fs, noOfFrames, frameLength, 3 ); + % fclose( statsFileID ); + % + % case 'calculatesampleRelativeEnergyGradientFeaturesByMeanAll' + % + % [x, fs, frameLength, noOfFrames] = openFile( currentSampleName ); + % getRelativeEnergyGradientFeatures( currentSampleName, x, fs, noOfFrames, frameLength ); + % + + end + + cd ../ + end + end + + end + cd ../ + + end + cd ../ + end + +end + +cd ../ +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Test/checkVUV.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,109 @@ +function [] = checkVUV() + +sampleName = '01_a_male01_pos'; +[x, fs, frameLength, noOfFrames] = openFile( [sampleName '.wav'] ); +y = buffer( x, frameLength ); + + %---------------- GET THE SILENT FRAME VALUES ------------------- + +% only wish to consider pitch values from voiced frames. +% silent and unvoiced frames will produce pitch values that +% are random and therefore will bias our results +segmentFrames = getSilenceDecision( sampleName ); +[ silentFrames ] = removeSilentData( segmentFrames, noOfFrames ); + + +% % open original annotation file +% fileName = [sampleName '.txt']; +% % read metrics from file +% fileID = fopen( fileName ); +% data = fscanf( fileID, '%d', inf ); +% annotatedPitch = data(3:4:end); +% fclose( fileID ); + +% open HNR file +fileName = [sampleName '_HNR.txt']; +% read metrics from file +fileID = fopen( fileName ); +data = fscanf( fileID, '%f', inf ); +harmonic2noise = data(2:2:end); +fclose( fileID ); + +largestH2NValue = max( [max(harmonic2noise) abs(min(harmonic2noise))] ); +%normalise +harmonic2noise = harmonic2noise/largestH2NValue; +%round to 1dp +harmonic2noise= (round(harmonic2noise*10))/10; + +% open audio power file +fileName = [sampleName '_AP.txt']; +% read metrics from file +fileID = fopen( fileName ); +data = fscanf( fileID, '%f', inf ); +audioPower = data(2:2:end); +fclose( fileID ); +maxPower = max(audioPower); + + +minFreq = ceil(fs/getVariables('getMinFreq')); %default minimum fundatmental frequency +maxFreq = ceil(fs/getVariables('getMaxFreq')); %default maximum fundamental frequency + +% open pitch file +fileName = [sampleName '_YIN_pitch.txt']; +% read metrics from file +fileID = fopen( fileName ); +data = fscanf( fileID, '%f', inf ); +pitch = data(2:2:end); +fclose( fileID ); + +% open band power file +% fileName = [sampleName '_ASBP_norm.txt']; +% % read metrics from file +% fileID = fopen( fileName ); +% data = fscanf( fileID, '%f', inf ); +% bandPowerRatio = data(2:2:end); +% fclose( fileID ); +% bandPowerRatioRMS = sqrt( mean(( bandPowerRatio/max(bandPowerRatio) ).^2 )); +% + + +HNRThresh = 0; +hold off; +minLength = min([length(harmonic2noise) length(audioPower) length(pitch)]); +plotUnits = 1/frameLength; +for( i=1:minLength ) + + + if( silentFrames(i) == 0 ) % harmonic2noise(i) == 0 ) + %silent +% plot( ((i-1)*frameLength)+1:(i*frameLength), y(:,i), 'y' );hold on; + plot( i-0.5 + plotUnits : plotUnits : i+0.5, y(:,i), 'y' );hold on; + vuv(i) = 0; + elseif( ((harmonic2noise(i) <= HNRThresh))... %|| (bandPowerRatio(i)/max(bandPowerRatio) < 0.001)) ... + || (pitch(i) > minFreq) ... + || (pitch(i) < maxFreq)) %... +% || (audioPower(i)/maxPower < 0.0005)) + % unvoiced +% plot( ((i-1)*frameLength)+1:(i*frameLength), y(:,i), 'r' );hold on; + plot( i-0.5 + plotUnits : plotUnits : i+0.5, y(:,i), 'r' );hold on; + vuv(i) = 2; + else + %voiced +% plot( ((i-1)*frameLength)+1:(i*frameLength), y(:,i), 'g' );hold on; + plot( i-0.5 + plotUnits : plotUnits : i+0.5, y(:,i), 'g' );hold on; + vuv(i) = 1; + end +end + + +% plot( 1 : 1 : length(harmonic2noise), harmonic2noise, 'b' );hold on; +% plot( 1 : 1 : length(pitch), pitch/max(pitch), 'w' );hold on; +% % pitch thresholds +% L=line([0 noOfFrames], [maxFreq/max(pitch) maxFreq/max(pitch)]); +% set(L,'color',[1 0 0]); +% L=line([0 noOfFrames], [minFreq/max(pitch) minFreq/max(pitch)]); +% set(L,'color',[1 0 0]); +% +% %HNR threshold +% L=line([0 noOfFrames], [HNRThresh HNRThresh]); +% set(L, 'color', [0 0 1] );
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Code/Test/plotAudioFrameByType.m Thu Jul 26 14:46:25 2012 +0100 @@ -0,0 +1,17 @@ +function plotAudioFrameByType( x, noOfValidFrames, vuv, frameLength ); + + +hold on + +for frameIdx = 1:noOfValidFrames + startPoint = ((frameIdx-1) * frameLength)+1; + endPoint = frameIdx * frameLength; +% if( vuv( frameIdx ) == 0 ) %silent +% plot( startPoint : endPoint , x( startPoint : endPoint ), 'w' ); + if( vuv( frameIdx ) == 1 ) % voiced + plot( startPoint : endPoint , x( startPoint : endPoint ), 'y' ); + elseif( vuv( frameIdx ) == 2 ) %unvoiced + plot( startPoint : endPoint , x( startPoint : endPoint ) , 'r' ); + end +end +hold off \ No newline at end of file