rmeddis@38: function [iih,IPIhisttime,IPIhistweight]=track_formants_from_IPI_guy(IFRAN_pattern, sfreq) rmeddis@38: % rmeddis@38: % tracks the formants according to an analysis proposed in Secker-Walker rmeddis@38: % JASA 1990, section V.A rmeddis@38: % Tim Jürgens, February 2011, code from Guy Brown included rmeddis@38: % rmeddis@38: % input: IFRAN_pattern: pattern of the auditory model (dependend on the number of modules used) rmeddis@38: % first dimension: frequency channel, rmeddis@38: % second dimension: time (samples) rmeddis@38: % sfreq: sampling frequency rmeddis@38: % output: iih: interpeak-interval histogram, matrix very similar rmeddis@38: % the plot 5 in the Secker-Walker paper rmeddis@38: % rmeddis@38: % rmeddis@38: % rmeddis@38: rmeddis@38: rmeddis@38: time_axis = 0:1/sfreq:(size(IFRAN_pattern,2)-1)/sfreq; rmeddis@38: rmeddis@38: %find how many samples of AN_pattern are 10ms and 3ms rmeddis@38: %one_sample_is_a_time_of = time_axis(2); rmeddis@38: [tmp, start_time_index] = min(abs(0-time_axis)); rmeddis@38: [tmp, stop10_time_index] = min(abs(0.01-time_axis)); rmeddis@38: number_of_samples10ms = stop10_time_index - start_time_index; rmeddis@38: rmeddis@38: [tmp, stop3_time_index] = min(abs(0.003-time_axis)); rmeddis@38: number_of_samples3ms = stop3_time_index - start_time_index; rmeddis@38: every_3ms = 1:number_of_samples3ms:size(IFRAN_pattern,2)-number_of_samples10ms; rmeddis@38: rmeddis@38: hamm_window = hamming(11); rmeddis@38: halfHamming = (length(hamm_window)-1)/2; rmeddis@38: rmeddis@38: % window normalization rmeddis@38: rmeddis@38: norm = conv(ones(1,floor(number_of_samples10ms/2)),hamm_window); rmeddis@38: norm = norm(5+1:end-5)'; rmeddis@38: win_size = number_of_samples10ms; rmeddis@38: half_win_size = floor(win_size/2); rmeddis@38: hop_size = number_of_samples3ms; rmeddis@38: rmeddis@38: rmeddis@38: %pre-allocation due to speed rmeddis@38: %Acorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1); rmeddis@38: %RAcorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1); rmeddis@38: %SRAcorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1-10); rmeddis@38: IPIhisttime = zeros(size(IFRAN_pattern,1),size(every_3ms,2),3); rmeddis@38: IPIhistweight = zeros(size(IFRAN_pattern,1),size(every_3ms,2),3); %maximum 3 peaks from the SRA rmeddis@38: iih = zeros(half_win_size,size(every_3ms,2)+1); rmeddis@38: rmeddis@38: rmeddis@38: rmeddis@38: rmeddis@38: for iCounter = 1:size(IFRAN_pattern,1) %each channel rmeddis@38: fprintf('Channel No. %i\n',iCounter); rmeddis@38: %time_counter = 1; rmeddis@38: %for jCounter = every_3ms %every 3ms time segment rmeddis@38: rmeddis@38: rmeddis@38: rmeddis@38: %% Guy's code rmeddis@38: % enframe this signal rmeddis@38: rmeddis@38: frames = enframe(IFRAN_pattern(iCounter,:),win_size,hop_size); rmeddis@38: rmeddis@38: % compute the autocorrelation rmeddis@38: rmeddis@38: acf = real(ifft(abs(fft(frames,[],2)).^2,[],2)); rmeddis@38: acf(acf<0)=0; rmeddis@38: acf = sqrt(acf(:,1:half_win_size)); rmeddis@38: rmeddis@38: % smooth with hamming window and take the root rmeddis@38: rmeddis@38: for frame=1:size(acf,1) rmeddis@38: rmeddis@38: %%debug rmeddis@38: %if iCounter == 130 rmeddis@38: % disp('here'); rmeddis@38: %end rmeddis@38: rmeddis@38: rmeddis@38: sra = conv(acf(frame,:),hamm_window); rmeddis@38: sra = sra(halfHamming+1:end-halfHamming)./norm'; rmeddis@38: df = [0 ; diff(sra')]; rmeddis@38: idx = find((df(1:end-1)>=0)&(df(2:end)<0)); rmeddis@38: % interpolate rmeddis@38: a=df(idx); rmeddis@38: b=df(idx+1); rmeddis@38: idx = (idx-1+a./(a-b)); rmeddis@38: % get rid of a zero peak, if it exists rmeddis@38: idx = idx(idx>1); rmeddis@38: % peak values corresponding to these intervals rmeddis@38: amp = interp1(1:length(sra),sra,idx,'linear'); rmeddis@38: % if required, remove peaks that lie below the mean sra rmeddis@38: % note that we disregard the value at zero delay rmeddis@38: %if (params.removePeaksBelowMean) rmeddis@38: valid = find(amp>mean(sra(2:end))); rmeddis@38: idx = idx(valid); rmeddis@38: amp = amp(valid); rmeddis@38: %end rmeddis@38: % only use the first four peaks (three intervals) rmeddis@38: idx = idx(1:min(4,length(idx))); rmeddis@38: % find the intervals rmeddis@38: interval = diff(idx); rmeddis@38: % now histogram the intervals rmeddis@38: if (~isempty(interval)) rmeddis@38: for k=1:length(interval), rmeddis@38: iih(round(interval(k)),frame) = iih(round(interval(k)),frame)+amp(k); rmeddis@38: IPIhisttime(iCounter,frame,k) = interval(k)/sfreq; rmeddis@38: IPIhistweight(iCounter,frame,k) = amp(k); rmeddis@38: end rmeddis@38: end rmeddis@38: rmeddis@38: end rmeddis@38: rmeddis@38: rmeddis@38: rmeddis@38: rmeddis@38: %% end Guy's code rmeddis@38: rmeddis@38: rmeddis@38: % %take the autocorrelation (ACF) of a 10ms-segment of each channel rmeddis@38: % Acorr(iCounter,time_counter,:) = xcorr(IFRAN_pattern(iCounter,jCounter:number_of_samples10ms+jCounter),'biased'); %biased scales the ACF by the reciprocal of the length of the segment rmeddis@38: % %root calculation rmeddis@38: % RAcorr(iCounter,time_counter,:) = sqrt(abs(Acorr(iCounter,time_counter,:))); rmeddis@38: % rmeddis@38: % %smoothing using the 11-point hamming window rmeddis@38: % for kCounter = 6:size(RAcorr(iCounter,time_counter,:),3)-5 %start with 6 and end with 5 samples rmeddis@38: % %less the length of time_axis not to get in conflict with the length of rmeddis@38: % %the hamm_window rmeddis@38: % SRAcorr(iCounter,time_counter,kCounter-5) = ... rmeddis@38: % squeeze(RAcorr(iCounter,time_counter,(kCounter-5):(kCounter+5)))'*hamm_window./sum(hamm_window); rmeddis@38: % end rmeddis@38: % rmeddis@38: % %mean value of actual SRA rmeddis@38: % SRA_mean = mean(SRAcorr(iCounter,time_counter,:)); rmeddis@38: % rmeddis@38: % %find signed zero-crossings of the first derivative (=difference) rmeddis@38: % z_crossings_indices = find(diff(sign(diff(squeeze(SRAcorr(iCounter,time_counter,:))))) < 0)+1; %+1 is necessary, because diff shortens vector by 1 rmeddis@38: % middle_index = ceil(size(SRAcorr(iCounter,time_counter,:),3)/2); rmeddis@38: % rmeddis@38: % validCounter = 1; rmeddis@38: % valid_z_crossings_indices = []; rmeddis@38: % %find valid zero-crossings (peak higher than meanvalue and within first 5 ms of SRA) rmeddis@38: % for lCounter = 1:length(z_crossings_indices) rmeddis@38: % if (SRAcorr(iCounter,time_counter,z_crossings_indices(lCounter)) > SRA_mean) && ... rmeddis@38: % (abs(z_crossings_indices(lCounter)-middle_index) < round(number_of_samples10ms/2)); rmeddis@38: % valid_z_crossings_indices(validCounter) = z_crossings_indices(lCounter); rmeddis@38: % validCounter = validCounter+1; rmeddis@38: % end rmeddis@38: % end rmeddis@38: % rmeddis@38: % %find main peak in the ACF rmeddis@38: % [tmp,index_of_z_crossings_main_index] = min(abs(middle_index-valid_z_crossings_indices)); rmeddis@38: % if ~tmp == 0 rmeddis@38: % disp('middle peak not appropriately found'); rmeddis@38: % end rmeddis@38: % rmeddis@38: % %%% for debugging rmeddis@38: % % if iCounter == 130 rmeddis@38: % % disp('here'); rmeddis@38: % % figure, plot(squeeze(SRAcorr(iCounter,time_counter,:))); rmeddis@38: % % hold on, plot([1 length(squeeze(SRAcorr(iCounter,time_counter,:)))],[SRA_mean SRA_mean],'r-'); rmeddis@38: % % end rmeddis@38: % %%% rmeddis@38: % rmeddis@38: % %generate IPI-histogram: take the first 3 intervals of SRAcorr rmeddis@38: % %(positive delay) in the first 5 ms rmeddis@38: % histcounter = 1; rmeddis@38: % for lCounter = index_of_z_crossings_main_index+1:min([length(valid_z_crossings_indices(index_of_z_crossings_main_index+1:end)) 3])+index_of_z_crossings_main_index rmeddis@38: % sampledifference = abs(valid_z_crossings_indices(lCounter)-valid_z_crossings_indices(lCounter-1)); rmeddis@38: % %the difference between two adjacent peaks in the SRA is taken rmeddis@38: % %as IPI estimate rmeddis@38: % IPIhisttime(iCounter,time_counter,histcounter) = sampledifference*one_sample_is_a_time_of; rmeddis@38: % %the amplitude of the SRA at the start of the SRA interval is rmeddis@38: % %taken as the IPIweight rmeddis@38: % IPIhistweight(iCounter,time_counter,histcounter) = SRAcorr(iCounter,time_counter,valid_z_crossings_indices(lCounter-1)); rmeddis@38: % histcounter = histcounter + 1; rmeddis@38: % end rmeddis@38: rmeddis@38: %time_counter = time_counter+1; rmeddis@38: end rmeddis@38: