rmeddis@38: function fach=fourierautocorrelationhistogram_direct(ANpattern,sfreq,plothandle)
rmeddis@38: 
rmeddis@38: 
rmeddis@38: time_axis = 0:1/sfreq:(size(ANpattern,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, stop20_time_index] = min(abs(0.020-time_axis));
rmeddis@38: number_of_samples20ms = stop20_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(ANpattern,2)-number_of_samples20ms;
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_samples20ms/2)),hamm_window);
rmeddis@38: norm = norm(5+1:end-5)';
rmeddis@38: win_size = number_of_samples20ms;
rmeddis@38: half_win_size = floor(win_size/2);
rmeddis@38: hop_size = number_of_samples3ms;
rmeddis@38: 
rmeddis@38: %preallocation due to speed
rmeddis@38: fach = zeros(half_win_size,size(every_3ms,2)+1);
rmeddis@38: 
rmeddis@38: for iCounter = 1:size(ANpattern,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(ANpattern(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(frames,1)
rmeddis@38:         
rmeddis@38:                
rmeddis@38:         smoothed_frame = conv(frames(frame,:),hamm_window);
rmeddis@38:         smoothed_frame = smoothed_frame(halfHamming+1:end-halfHamming);
rmeddis@38:         fsra = 20*log10(abs(fft(smoothed_frame-mean(smoothed_frame))));
rmeddis@38:         fsra = fsra(1:floor(length(fsra)/2));
rmeddis@38:         
rmeddis@38:         t = [0:1/sfreq:length(smoothed_frame)/sfreq-1/sfreq];
rmeddis@38:         frequency = [0:1/t(end):1/(2*(t(2)-t(1)))];
rmeddis@38:         %identify peaks in the fft
rmeddis@38:         df = [0 ; diff(fsra')];
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:         [sorted,sortedindex]=sort(fsra(idx),'descend');
rmeddis@38:         % just take the three highest values of the fourier-transform
rmeddis@38:          valid_peak_index = sortedindex(1:min([length(sortedindex) 1]));
rmeddis@38:          amp = sorted(1:min([length(sortedindex) 1]));
rmeddis@38:          
rmeddis@38:          %store valid peaks according to amplitude in a histogram
rmeddis@38:          if (~isempty(valid_peak_index))
rmeddis@38:             for k=1:length(valid_peak_index),
rmeddis@38:                 fach(idx(valid_peak_index(k)),frame) = fach(idx(valid_peak_index(k)),frame)+amp(k);            
rmeddis@38:             end
rmeddis@38:         end
rmeddis@38:          %transform index into frequencies
rmeddis@38:          
rmeddis@38:     end
rmeddis@38: end
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %plot the result
rmeddis@38: if ~exist('plothandle'), plothandle=figure; end
rmeddis@38: maxfrequency = 4000;
rmeddis@38: [tmp,number_of_channels_to_display] = min(abs(frequency-maxfrequency));
rmeddis@38: frequency = frequency(1:number_of_channels_to_display);
rmeddis@38: 
rmeddis@38: set(gcf,'Currentaxes',plothandle);
rmeddis@38: 
rmeddis@38: YTickIdx = 1:floor(numel(frequency)/6):numel(frequency);
rmeddis@38: XTickIdx = 1:floor(numel(every_3ms)/6):numel(every_3ms);
rmeddis@38: YTickIdxRev = numel(frequency)+1-YTickIdx;
rmeddis@38: if ~isempty(gca)
rmeddis@38:     axes(gca);  %#ok<MAXES>
rmeddis@38:     imagesc(fach(1:number_of_channels_to_display,:));
rmeddis@38:     axis xy
rmeddis@38:     set(gca, 'YTick', YTickIdx);
rmeddis@38:     set(gca, 'YTickLabel', num2str(   frequency(YTickIdx)', '%0.0f' ));
rmeddis@38:     ylabel('frequency (Hz)')
rmeddis@38:     set(gca, 'XTick', XTickIdx);
rmeddis@38:     set(gca, 'XTickLabel', XTickIdx.*3);
rmeddis@38:     xlabel('Time (ms)')
rmeddis@38: end
rmeddis@38: 
rmeddis@38: 
rmeddis@38: 
rmeddis@38: