Mercurial > hg > map

comparison userProgramsTim/IPIHextract.m @ 38:c2204b18f4a2 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents
children
comparison
equal deleted inserted replaced
37:771a643d5c29 38:c2204b18f4a2
1 function [iih,IPIhisttime,IPIhistweight]=IPIHextract(IFRAN_pattern, sfreq)
2 %
3 % tracks the formants according to an analysis proposed in Secker-Walker
4 % JASA 1990, section V.A
5 % Tim Jürgens, February 2011, code from Guy Brown included
6 %
7 % input: IFRAN_pattern: pattern of the auditory model (dependend on the number of modules used)
8 % first dimension: frequency channel,
9 % second dimension: time (samples)
10 % sfreq: sampling frequency
11 % output: iih: interpeak-interval histogram, matrix very similar
12 % the plot 5 in the Secker-Walker paper
13 %
14 %
15 %
16
17
18 time_axis = 0:1/sfreq:(size(IFRAN_pattern,2)-1)/sfreq;
19
20 %find how many samples of AN_pattern are 10ms and 3ms
21 %one_sample_is_a_time_of = time_axis(2);
22 [tmp, start_time_index] = min(abs(0-time_axis));
23 [tmp, stop20_time_index] = min(abs(0.020-time_axis));
24 number_of_samples20ms = stop20_time_index - start_time_index;
25
26 [tmp, stop10_time_index] = min(abs(0.010-time_axis));
27 number_of_samples10ms = stop10_time_index - start_time_index;
28 every_10ms = 1:number_of_samples10ms:size(IFRAN_pattern,2)-number_of_samples20ms;
29
30 hamm_window = hamming(11);
31 halfHamming = (length(hamm_window)-1)/2;
32
33 % window normalization
34
35 norm = conv(ones(1,number_of_samples20ms),hamm_window);
36 norm = norm(5+1:end-5)';
37 win_size = number_of_samples20ms;
38 half_win_size = floor(win_size/2);
39 hop_size = number_of_samples10ms;
40
41 %parameters of the autocorrelation
42 params.acfTau=0.1;
43 params.lags=[0:1/sfreq:0.02-1/sfreq];
44 sampledacf = runningACF(IFRAN_pattern,sfreq,params);
45 sampledacf(sampledacf<0)=0;
46 sampledacf = sqrt(sampledacf);
47
48
49 %pre-allocation due to speed
50 %Acorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1);
51 %RAcorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1);
52 %SRAcorr = zeros(size(IFRAN_pattern,1),size(every_3ms,2),number_of_samples10ms*2+1-10);
53 IPIhisttime = zeros(size(IFRAN_pattern,1),size(every_10ms,2),3);
54 IPIhistweight = zeros(size(IFRAN_pattern,1),size(every_10ms,2),3); %maximum 3 peaks from the SRA
55 iih = zeros(half_win_size,size(sampledacf,1));
56
57
58
59
60
61 for iCounter = 1:size(sampledacf,2) %each channel
62 fprintf('Channel No. %i\n',iCounter);
63 %time_counter = 1;
64 %for jCounter = every_3ms %every 3ms time segment
65
66 for frame=1:size(sampledacf,1)
67
68 %%debug
69 %if iCounter == 130
70 % disp('here');
71 %end
72
73
74 sra = conv(squeeze(sampledacf(frame,iCounter,:)),hamm_window);
75 sra = sra(halfHamming+1:end-halfHamming)./norm;
76 df = [0 ; diff(sra)];
77 idx = find((df(1:end-1)>=0)&(df(2:end)<0));
78 % interpolate
79 a=df(idx);
80 b=df(idx+1);
81 idx = (idx-1+a./(a-b));
82 % get rid of a zero peak, if it exists
83 idx = idx(idx>1);
84 %include the zeroth peak
85 idx = [1 idx']';
86 % peak values corresponding to these intervals
87 amp = interp1(1:length(sra),sra,idx,'linear');
88 % if required, remove peaks that lie below the mean sra
89 % note that we disregard the value at zero delay
90 %if (params.removePeaksBelowMean)
91 valid = find(amp>1.2*mean(sra(1:floor(length(sra)/2))));
92 %valid = find(amp>mean(sra));
93 %just take the mean of the first half of the sra as a comparison
94 idx = idx(valid);
95 amp = amp(valid);
96 %end
97 % only use the first four peaks (three intervals)
98 idx = idx(1:min(4,length(idx)));
99 % find the intervals
100 interval = diff(idx);
101 % now histogram the intervals
102 if (~isempty(interval))
103 for k=1:length(interval)
104 if interval(k)<=half_win_size
105 iih(round(interval(k)),frame) = iih(round(interval(k)),frame)+amp(k);
106 IPIhisttime(iCounter,frame,k) = interval(k)/sfreq;
107 IPIhistweight(iCounter,frame,k) = amp(k);
108 end
109 end
110 end
111
112 end
113
114
115
116
117 %% end Guy's code
118
119
120 % %take the autocorrelation (ACF) of a 10ms-segment of each channel
121 % 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
122 % %root calculation
123 % RAcorr(iCounter,time_counter,:) = sqrt(abs(Acorr(iCounter,time_counter,:)));
124 %
125 % %smoothing using the 11-point hamming window
126 % for kCounter = 6:size(RAcorr(iCounter,time_counter,:),3)-5 %start with 6 and end with 5 samples
127 % %less the length of time_axis not to get in conflict with the length of
128 % %the hamm_window
129 % SRAcorr(iCounter,time_counter,kCounter-5) = ...
130 % squeeze(RAcorr(iCounter,time_counter,(kCounter-5):(kCounter+5)))'*hamm_window./sum(hamm_window);
131 % end
132 %
133 % %mean value of actual SRA
134 % SRA_mean = mean(SRAcorr(iCounter,time_counter,:));
135 %
136 % %find signed zero-crossings of the first derivative (=difference)
137 % z_crossings_indices = find(diff(sign(diff(squeeze(SRAcorr(iCounter,time_counter,:))))) < 0)+1; %+1 is necessary, because diff shortens vector by 1
138 % middle_index = ceil(size(SRAcorr(iCounter,time_counter,:),3)/2);
139 %
140 % validCounter = 1;
141 % valid_z_crossings_indices = [];
142 % %find valid zero-crossings (peak higher than meanvalue and within first 5 ms of SRA)
143 % for lCounter = 1:length(z_crossings_indices)
144 % if (SRAcorr(iCounter,time_counter,z_crossings_indices(lCounter)) > SRA_mean) && ...
145 % (abs(z_crossings_indices(lCounter)-middle_index) < round(number_of_samples10ms/2));
146 % valid_z_crossings_indices(validCounter) = z_crossings_indices(lCounter);
147 % validCounter = validCounter+1;
148 % end
149 % end
150 %
151 % %find main peak in the ACF
152 % [tmp,index_of_z_crossings_main_index] = min(abs(middle_index-valid_z_crossings_indices));
153 % if ~tmp == 0
154 % disp('middle peak not appropriately found');
155 % end
156 %
157 % %%% for debugging
158 % % if iCounter == 130
159 % % disp('here');
160 % % figure, plot(squeeze(SRAcorr(iCounter,time_counter,:)));
161 % % hold on, plot([1 length(squeeze(SRAcorr(iCounter,time_counter,:)))],[SRA_mean SRA_mean],'r-');
162 % % end
163 % %%%
164 %
165 % %generate IPI-histogram: take the first 3 intervals of SRAcorr
166 % %(positive delay) in the first 5 ms
167 % histcounter = 1;
168 % 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
169 % sampledifference = abs(valid_z_crossings_indices(lCounter)-valid_z_crossings_indices(lCounter-1));
170 % %the difference between two adjacent peaks in the SRA is taken
171 % %as IPI estimate
172 % IPIhisttime(iCounter,time_counter,histcounter) = sampledifference*one_sample_is_a_time_of;
173 % %the amplitude of the SRA at the start of the SRA interval is
174 % %taken as the IPIweight
175 % IPIhistweight(iCounter,time_counter,histcounter) = SRAcorr(iCounter,time_counter,valid_z_crossings_indices(lCounter-1));
176 % histcounter = histcounter + 1;
177 % end
178
179 %time_counter = time_counter+1;
180 end
181