· 4 years ago · Sep 11, 2021, 09:16 PM
1%% 1.1.Tony INIT & build all structures
2clc
3close all
4% clear all
5clear notes hums seeds chains indicesToNanify allChains allChainsAndSeeds
6
7% select experiment version - and deduce version-specific experiment parameters
8vers = '3'; %1, 2 or 3; before 4aug21, these versions were called 1a, 1b and 2 resp.
9gender = 'm'; %m or f; v1 (formerly v1a) only exists in a f version
10switch vers
11 case {'1', '2', '9'}
12 N_notesInSeed = 5;
13 duration_seed = 5; % in seconds
14 case '3'
15 N_notesInSeed = 13;
16 duration_seed = 6.5;
17end
18
19% parameters that apply to all 3 experiment versions
20N_chains = 4;
21N_seeds = 12;
22N_gens = 8;
23
24load 'd:\Experiments\stimuli\seedStimuli_fqs.mat' % contains variables (arrays) such as seed_v1f, seeds_v2m etc
25path_hums = [ 'd:\Experiments\Analyses\piano rolls from manual export of note-data in Tony, by Jonathan\'];
26path_TIV = 'd:\Experiments\Analyses\TIV measures (diatonicity etc)\';
27
28% plot preferences: aspect etc
29colour = 'blue';
30titleX = -0.145; % coordinates for labels to be added to the left of each subplot
31titleY = 0.3;
32maxStOnY = 7;
33fontSize_suptitle = 12;
34
35% other prefs
36strSimMetric = 'edit distance'; % aka Levenshtein
37
38% define bins for histograms
39% bin centres are mid-points between edges, and vice-versa; and the mid-point (on paper/in space) between two values plotted on a log-axis is the sqrt of their product (their geometric mean)
40% the first and last bin of each histogram will therefore not have a left/right edge resp. defined for them in this way, and we convene to define them as being equal to the bin's centre; the effect (loss of accuracy) of this should be negligible.
41% thus, edge nr 1 = centre nr 1, and edge nr N+1 = centre nr N, where N is the last bin; also, edge nr i will be a right-edge for bin nr i-1, and a left-edge for bin nr i; and for N bins, there will be N+1 edges, of which 2 (first, last) defined as above
42% 1 2 3 4 5 6 7 8 9 10 11 12
43% 'C2' 'C#2' 'D2' 'D#2' 'E2' 'F2' 'F#2' 'G2' 'G#2' 'A2' 'A#2' 'B2'
44pianoNotes = [65.41 69.3 73.42 77.78 82.41 87.31 92.5 98 103.83 110 116.54 123.47] ./ 2; % divide by 2 to get the smallest possible pitches, i.e. the lowest octave range
45pianoNotes_rep = [pianoNotes 2*pianoNotes 4*pianoNotes 8*pianoNotes 16*pianoNotes]'; % all possible/reasonable pitches that can be hummed
46
47% use TIVlib (in Python/visual studio) to compute diatonicity and chromaticity directly from WAVs
48fileName_D = [ path_TIV 'diatonicity_' vers gender '.csv'];
49D = readtable(fileName_D, 'readvariablenames', false, 'FileType','text'); % the M variable is now in table format
50D = table2array(D); % the M variable is now in numerical format
51fileName_C = [ path_TIV 'chromaticity_' vers gender '.csv'];
52C = readtable(fileName_C, 'readvariablenames', false, 'FileType','text'); % the M variable is now in table format
53C = table2array(C); % the M variable is now in numerical format
54
55% init
56[ seeds(1:N_seeds).parsons ] = deal({'*'});
57r=0; % row into diatonicity/chr vector from TIV
58for i_seed = 1:N_seeds % this asignement for nested indexing sadly only works through looping, esp since cell arrays are also involved, see cz matlab
59 for i_chain = 1:N_chains
60 for i_gen = 1:N_gens
61 hums(i_seed, i_chain).parsons(i_gen) = {'*'};
62 end
63 end
64end
65
66% build seeds structure
67for i_seed = 1:N_seeds
68 % create the seeds structure
69 varName=genvarname(['seeds_' vers gender]);
70 M = eval(varName); % these variable (array) names, such as seeds_v1a_female seeds_v1b_male etc, are loaded earlier from seedStimuli_fqs.mat
71 pitches = M(i_seed, :); % pitches is a vector; in M, rows=seeds, columns=notes
72 seeds(i_seed).pitch = pitches;
73 ratio = seeds(i_seed).pitch ./ 261.63; % pitch difference (in fact ratio) between the notes' pitches, and middle C (261.63 Hz)
74 seeds(i_seed).pitch_stC = 12*log2(ratio); % convert this to semitones above/below middle C
75 seeds(i_seed).noteDuration = duration_seed/N_notesInSeed; % in seconds; for a given expt version, all notes in the seed have a constant duration, hence not indexed
76
77 % build PC field (pitch class) for seeds
78 for i_note = 1:N_notesInSeed
79 [~, index] = min(abs(pitches(i_note)-pianoNotes_rep(:))); % TO DO: this min actually needs to measure along a log-axis not a linear one, thus this needs a few more lines of code
80 if mod(index, 12)~=0 % for convenience/clarity, we index the last note of a chromatic scale (the note B) not at 0 (the result of 12 mod 12) but actually at 12. Any reminaing 0 values in there are then equivalent to NaNs , i.e. there was no more note there
81 seeds(i_seed).PC(i_note) = mod(index, 12);
82 else
83 seeds(i_seed).PC(i_note) = 12;
84 end
85 end
86
87 % one more "for i_note", to compute transition-related measures: intervals and Parsons codes; this time we start from 2 as we compare every note with its predecessor
88 for i_note = 2:N_notesInSeed % we start from 2 as we compare with i_note-1 below
89 if pitches(i_note) > pitches(i_note-1)
90 seeds(i_seed).parsons = strcat(seeds(i_seed).parsons, 'u');
91 elseif pitches(i_note) < pitches(i_note-1)
92 seeds(i_seed).parsons = strcat(seeds(i_seed).parsons, 'd');
93 else % equality.. this should never really happen..
94 seeds(i_seed).parsons = strcat(seeds(i_seed).parsons, 'r');
95 end
96
97 % define intervals structure for transitions between consecutive notes (note-bigrams)
98 seeds(i_seed).interval_ratio(i_note-1) = pitches(i_note) / pitches(i_note-1); % 8 notes gives 7 intervals, so i_note goes 2:8 while i_interval would go 1:7
99 seeds(i_seed).interval_st (i_note-1) = 12*log2(seeds(i_seed).interval_ratio(i_note-1)); % convert to semitones;
100 end
101end
102
103% first do a for loop to build the essential fields of the notes structure, directly pulled from the CSV
104for i_seed=1:N_seeds
105 for i_chain = 1:N_chains
106 for i_gen = 1:N_gens
107 fileName = [ path_hums vers gender '_' num2str(i_seed, '%02.0f') '_' num2str(i_chain) '_' num2str(i_gen) '.csv']; % Jonathan's files from Tony
108 M = readtable(fileName, 'readvariablenames', false, 'FileType','text'); % the M variable is now in table format
109 M = table2array(M); % the M variable is now in numerical format
110 hums(i_seed, i_chain).N_notes (i_gen) = size(M,1); % typically this wont stay equal to N_notesInSeed, as even gen.1 will probably hum either fewer or more notes then in the seed, and further down the line this is likely to change even more
111 hums(i_seed, i_chain).N_intervals (i_gen) = hums(i_seed, i_chain).N_notes(i_gen) - 1; % this will always be the nr of notes minus 1. NB : this and the above vae are not affected by the artefactual (matlab-padding-induced) 0s later to be replacced with NaNs, as they are extracted straight from the csv; thus, theres no need to later nanify any indices for these scalar fields as we do for the oter (2D) fields
112 hums(i_seed, i_chain).pairs.intervals.N(i_gen) = hums(i_seed, i_chain).N_intervals(i_gen) - 1;
113 hums(i_seed, i_chain).pairs.notes.N (i_gen) = hums(i_seed, i_chain).N_notes (i_gen) - 1;
114
115 for i_note = 1:hums(i_seed, i_chain).N_notes(i_gen) % rows in M correspond to notes
116 notes(i_seed, i_chain).t1 (i_gen, i_note) = M(i_note, 1); % 1st column in those CSVs is note onset
117 notes(i_seed, i_chain).pitch (i_gen, i_note) = M(i_note, 2); % 2nd is note pitch
118 notes(i_seed, i_chain).duration (i_gen, i_note) = M(i_note, 3); % 3rd is note duration
119 end
120 end
121
122 % now that all notes from all generations have been defined, for the current seed/chain..
123 % for generations with notes fewer than the max nr of notes in the chain, their rows get padded with 0 up until that max nr
124 % these 0s can create confusion later on, so we replace those with NaNs
125 % we go by the pitch field, and later correct all others; where there was 0 for that field, therell be one also for other fields computed based on it
126 A = notes(i_seed, i_chain).pitch(:,:);
127 [rows,columns] = find(A==0);
128 indicesToNanify(i_seed, i_chain).rows = rows;
129 indicesToNanify(i_seed, i_chain).columns = columns;
130 end
131end
132
133% do another 3xnested for loop, to compute the other fields for 'notes' and 'hums'; we take care of the indices to nanify in the next loop
134% cred ca de fapt poate fi pus si asta DUPA acel pas de nanificare..
135for i_seed=1:N_seeds
136 for i_chain = 1:N_chains
137 for i_gen = 1:N_gens
138 pcAirtime = []; % one for each gen
139 hums(i_seed, i_chain).N_inflections(i_gen) = 0;
140 for i_note = 1:hums(i_seed, i_chain).N_notes(i_gen) % rows in M correspond to notes
141 notes(i_seed, i_chain).t2 (i_gen, i_note) = notes(i_seed, i_chain).t1(i_gen, i_note) + notes(i_seed, i_chain).duration(i_gen, i_note);
142
143 ratio = notes(i_seed, i_chain).pitch(i_gen, i_note) / 261.63; % pitch difference (in fact ratio) between this note and middle C (261.63 Hz)
144 notes(i_seed, i_chain).pitch_stC (i_gen, i_note) = 12*log2(ratio); % convert this to semitones above/below middle C
145
146 % build PC field (pitch class) for hums
147 pitch = notes(i_seed, i_chain).pitch(i_gen, i_note);
148 [~, index] = min(abs(pitch-pianoNotes_rep(:)));
149 if mod(index, 12)~=0 % for convenience/clarity, we index the last note of a chromatic scale (the note B) not at 0 (the result of 12 mod 12) but actually at 12. Any reminaing 0 values in there are then equivalent to NaNs , i.e. there was no more note there
150 notes(i_seed, i_chain).PC (i_gen, i_note) = mod(index, 12);
151 else
152 notes(i_seed, i_chain).PC (i_gen, i_note) = 12;
153 end
154
155 % build the airtime-weighted version of the pitches_highpass vector; NOTE THAT all the _highpass variables are assumed to've been computed in a previous cell (the one with the main plots), ca sa nu mai repet aici codul ala
156 ms = round ( 1000*notes(i_seed, i_chain).duration(i_gen, i_note) );
157 pcAirtime = [ pcAirtime repmat(notes(i_seed, i_chain).PC (i_gen, i_note), [1 ms]) ]; % repeat the current note's pitch for as many times as there are ms in its duration, to make a very long vector
158 end
159
160 hums(i_seed, i_chain).pcAirtime(i_gen, 1:length(pcAirtime)) = pcAirtime;
161
162 % one more "for i_note", to compute transition-related measures: intervals, Parsons codes, and (Martins idea) inflection points; this time we start from 2 as we compare every note with its predecessor
163 for i_note = 2:hums(i_seed, i_chain).N_notes(i_gen)
164 % parsons contour
165 if notes(i_seed, i_chain).pitch(i_gen, i_note) > notes(i_seed, i_chain).pitch(i_gen, i_note-1)
166 hums(i_seed, i_chain).parsons(i_gen) = strcat(hums(i_seed, i_chain).parsons(i_gen), 'u');
167 string = cell2mat(hums(i_seed, i_chain).parsons(i_gen));
168 if string(end-1)=='d' % if the previous direction was opposite, means we've just had an inflection point
169 hums(i_seed, i_chain).N_inflections(i_gen) = hums(i_seed, i_chain).N_inflections(i_gen) + 1;
170 end
171 elseif notes(i_seed, i_chain).pitch(i_gen, i_note) < notes(i_seed, i_chain).pitch(i_gen, i_note-1)
172 hums(i_seed, i_chain).parsons(i_gen) = strcat(hums(i_seed, i_chain).parsons(i_gen), 'd');
173 string = cell2mat(hums(i_seed, i_chain).parsons(i_gen));
174 if string(end-1)=='u' % if the previous direction was opposite, means we've just had an inflection point
175 hums(i_seed, i_chain).N_inflections(i_gen) = hums(i_seed, i_chain).N_inflections(i_gen) + 1;
176 end
177 else % equality.. given the fine resolution of pitchj values, it's very unlikely we'll ever go down this arm of the if
178 hums(i_seed, i_chain).parsons(i_gen) = strcat(hums(i_seed, i_chain).parsons(i_gen), 'r');
179 end
180
181 % define intervals structure for transitions between consecutive notes (note-bigrams)
182 hums(i_seed, i_chain).interval_ratio(i_gen, i_note-1) = notes(i_seed, i_chain).pitch(i_gen, i_note) / notes(i_seed, i_chain).pitch(i_gen, i_note-1); % 8 notes gives 7 intervals, so i_note goes 2:8 while i_interval would go 1:7
183 hums(i_seed, i_chain).interval_st (i_gen, i_note-1) = 12*log2(hums(i_seed, i_chain).interval_ratio(i_gen, i_note-1)); % convert to semitones;
184 end
185
186 % hum duration, ascertained as the offset of the last note; the actual duration of the WAV might be a bit different, if there are unpitched portions of it following the offset of the last note
187 N_notes = hums(i_seed, i_chain).N_notes(i_gen); % how many notes in that gen; creating this var just for ease of reading the line below
188 hums(i_seed, i_chain).duration(i_gen) = notes(i_seed, i_chain).t2(i_gen, N_notes);
189
190 % (this doesnt require an i_note loop) compute string similarities between the Parsons code of the current hum and that of the previous one (or of the seed, if we are at gen.1)
191 if i_gen==1
192 switch strSimMetric
193 case 'edit distance'
194 hums(i_seed, i_chain).distFromPrev(i_gen) = EditDistance ( cell2mat(hums(i_seed, i_chain).parsons(i_gen)), cell2mat(seeds(i_seed).parsons) ); % the strings to be compared are cell arrays, but EditDistance needs normal strings as inouts
195 end
196 else
197 switch strSimMetric
198 case 'edit distance'
199 hums(i_seed, i_chain).distFromPrev(i_gen) = EditDistance ( cell2mat(hums(i_seed, i_chain).parsons(i_gen)), cell2mat(hums(i_seed, i_chain).parsons(i_gen-1)) );
200 end
201 end
202
203 % derive first and last pitches
204 lastNonZero = max(find(notes(i_seed, i_chain).pitch(i_gen, :)~=0)); % when you start to get 0s in this vector, means youve run out of valid notes
205 hums(i_seed, i_chain).lastPitch (i_gen) = notes(i_seed, i_chain).pitch(i_gen, lastNonZero);
206 hums(i_seed, i_chain).firstPitch(i_gen) = notes(i_seed, i_chain).pitch(i_gen, 1 );
207
208 % compute the meldevFromPrev metric of invariance of pitch intervals (to what degree do melodies stay rigid and merely get transposed, vs deviating?). This metric can then be plotted with P1.1
209 % this analysis really only makes sense for v1, as we need to assume constnat nr of notes across hums, namely 5 notes
210 if vers==1
211 if i_gen==1
212 ratios = notes(i_seed, i_chain).pitch(i_gen, 1:5) ./ seeds(i_seed).pitch(1:5);
213 else
214 ratios = notes(i_seed, i_chain).pitch(i_gen, 1:5) ./ notes(i_seed, i_chain).pitch(i_gen-1, 1:5);
215 end
216 ratiosDemeaned = ratios - mean(ratios); % demean
217 hums(i_seed, i_chain).meldevFromPrev(i_gen) = sum (abs(ratiosDemeaned)); % a scalar value; I tried with and without abs, but truth is, this should really be with signed values, i.e. no abs
218 end
219
220 % use TIVlib to compute diatonicity and chromaticity directly from WAVs
221 r=r+1; % starts at 0 up in the init commands, for each expt version
222 hums(i_seed, i_chain).diatonicity(i_gen) = D(r);
223 hums(i_seed, i_chain).chromaticity(i_gen) = C(r);
224 end
225 end
226end
227
228% turn all the flagged 0s into NaNs (see above)
229% first the notes structure, then hums
230% check that there arent any fields that get computed above but not included in the lists below; better yet, just compute any new fields in the loops below, to make sure the base fields that are used have been nanified
231for i_seed=1:N_seeds
232 for i_chain = 1:N_chains
233 for field={'pitch', 't1', 't2', 'duration', 'pitch_stC', 'PC'}
234 field = cell2mat(field); % the dynamic structure reference below must have field names referred to as strings, not cells; it's ok to change this inside its very own loop, as we simply change the format
235 rows = indicesToNanify(i_seed, i_chain).rows;
236 columns = indicesToNanify(i_seed, i_chain).columns;
237 for i=1:length(rows) % same as length(columns) i hope; note that if you were to use the following non-loop command, it doesnt take pairs of valies from rows(1)&columns(1) etc, instead it makes an unexpected mess: notes(i_seed, i_chain).(field)(rows, columns) = NaN;
238 notes(i_seed, i_chain).(field)(rows(i), columns(i)) = NaN;
239 end
240 end
241 end
242end
243for i_seed=1:N_seeds
244 for i_chain = 1:N_chains
245 for field={'interval_ratio', 'interval_st'}
246 field = cell2mat(field); % dynamic structure reference must have field names referred to as strings, not cells
247 rows = indicesToNanify(i_seed, i_chain).rows;
248 columns = indicesToNanify(i_seed, i_chain).columns - 1; %-1 bcs there is one less intervals than there are notes
249 % if you use this command, it doesnt take pairs of valies from rows(1)&columns(1) etc, instead it makes an unexpected mess: hums(i_seed, i_chain).(field)(rows, columns) = NaN;
250 for i=1:length(rows) % same as length(columns) i hope
251 hums(i_seed, i_chain).(field)(rows(i), columns(i)) = NaN;
252 end
253 end
254 end
255end
256
257% with 0s in notes and hums properly Nanified, we can now build from these the remaining structures needed
258% first, "allChainsAndSeeds" to later allow histogramming data from across all gens, all chains and all seeds
259allChainsAndSeeds.lastPitch = []; % note this structure does not have the usual (seed,chain) indices as it pools across those
260allChainsAndSeeds.firstPitch = [];
261allChainsAndSeeds.pcAirtime = [];
262allChainsAndSeeds.interval_st = [];
263for i_seed=1:N_seeds
264 allChains(i_seed).lastPitch = []; % allChains contains data from across all gens and all chains
265 allChains(i_seed).firstPitch = [];
266 allChains(i_seed).pcAirtime = [];
267 allChains(i_seed).interval_st = [];
268 for i_chain = 1:N_chains
269 allChains(i_seed).lastPitch = [ allChains(i_seed).lastPitch hums(i_seed, i_chain).lastPitch(:)' ];
270 allChains(i_seed).firstPitch = [ allChains(i_seed).firstPitch hums(i_seed, i_chain).firstPitch(:)' ];
271 allChains(i_seed).pcAirtime = [ allChains(i_seed).pcAirtime reshape(hums(i_seed, i_chain).pcAirtime(:, :), 1, [])]; % we need reshape here (and in the above 2 lines not) bcs airtime is generated by repeating values within a gen, whereas final/first notes are just single scalars
272 allChains(i_seed).interval_st = [ allChains(i_seed).interval_st reshape(hums(i_seed, i_chain).interval_st(:, :), 1, [])]; % as above, since intervals are not scalar for a given gen, but vary across the hum, thus they are vectors. NB: this interval_st field gets created by parsing the "hums(i_seed, i_chain).interval_st" matrix first across rows (generations) and then across columns (notes)
273 end
274
275 allChainsAndSeeds.lastPitch = [ allChainsAndSeeds.lastPitch allChains(i_seed).lastPitch];
276 allChainsAndSeeds.firstPitch = [ allChainsAndSeeds.firstPitch allChains(i_seed).firstPitch];
277 allChainsAndSeeds.pcAirtime = [ allChainsAndSeeds.pcAirtime allChains(i_seed).pcAirtime];
278 allChainsAndSeeds.interval_st = [ allChainsAndSeeds.interval_st allChains(i_seed).interval_st];
279end
280
281% compute structures that allow plotting (in cell 7.3) of 2D transition matrices (bigrams of notes and intervals), and also entropy etc
282for i_seed=1:N_seeds
283 for i_chain = 1:N_chains
284 for i_gen = 1:N_gens
285 % interval bigrams (interval pairs)
286 N_intervals = hums(i_seed, i_chain).N_intervals(i_gen); % shorthand variable
287 N_intervalPairs = N_intervals - 1; % for N elements there will be N-1 bigrams thereof; referring to these as interval pairs rather than bigrams to avoid confusion
288 for i_interval = 1:N_intervals - 1
289 % compute 2D transition matrices (bigrams)
290 % where the interval_st value is (has been) nanified, the corresponding intervalPair will consist of [] both for the 1st and the 2nd interval, thus these are good to plot now
291 i_pair = i_interval;
292 hums(i_seed, i_chain).pairs.intervals.element1(i_gen, i_pair) = hums(i_seed, i_chain).interval_st(i_gen, i_interval); % measured in semitones (not semitones above middle C, like notes, as we have a difference here)!
293 hums(i_seed, i_chain).pairs.intervals.element2(i_gen, i_pair) = hums(i_seed, i_chain).interval_st(i_gen, i_interval+1);
294 end
295
296 % note bigrams (note pairs)
297 N_notes = hums(i_seed, i_chain).N_notes(i_gen); % shorthand variable
298 N_notePairs = N_notes - 1; % for N elements there will be N-1 bigrams thereof; referring to these as note pairs rather than bigrams to avoid confusion
299 for i_note = 1:N_notes - 1
300 % compute 2D transition matrices (bigrams)
301 % where the pitch_stC value is (has been) nanified, the corresponding notePair will consist of [] both for the 1st and the 2nd note, thus these are good to plot now
302 i_pair = i_note;
303 hums(i_seed, i_chain).pairs.notes.element1(i_gen, i_pair) = notes(i_seed, i_chain).pitch_stC(i_gen, i_note); % measured in semitones above middle C in this case
304 hums(i_seed, i_chain).pairs.notes.element2(i_gen, i_pair) = notes(i_seed, i_chain).pitch_stC(i_gen, i_note+1);
305 end
306
307 % compute entropy for intervals
308 signal = hums(i_seed, i_chain).interval_st(i_gen, 1:end);
309 h1 = histogram(signal, 'Normalization', 'Probability'); % either histogram(signal, 'Normalization', 'Probability') or hist(signal); cf https://www.mathworks.com/matlabcentral/answers/1449894-cannot-use-histogram-to-compute-entropy
310 probabilities = h1.Values; % h1.Values if above you used 'histogram', just h1 if you used 'hist'. Either way, this returns the values that go in each of the bins (whatever those bins are determined to be)
311 probabilities(find(probabilities(:)==0)) = []; % when stored in the field matrix, there will be 0s at the end for padding, but right now in probabilities there are normally no 0s, except for some odd cases where this happens for some reason; here we remove any non-end (mid-vector) 0s, as they lead to a NaN entropy, given the logarithm. This should still lead to 4 elements for each probabilities vector (though not sure why). These 0s probably occur due to the low number of samples in the signal; cf https://www.mathworks.com/matlabcentral/answers/1449894-cannot-use-histogram-to-compute-entropy
312 hums(i_seed, i_chain).probabilities (i_gen, 1:length(probabilities)) = probabilities; % dont really need to store this as a field, doing it to debug the NaNs that I get below in the entropy field
313 hums(i_seed, i_chain).entropy_interval(i_gen) = -sum( probabilities .* log2(probabilities) ); % entropy measure based on intervals
314 end
315 end
316end
317close all % we dont need the plot that opens after the 'histogram' call above
318disp('gata...........................................................................................................................................................................')
319
320%% 2.1.Tony plot seed+hums of each chain, x axis = time (uses 'notes' struc-array) [[this cell to be made into a P cell]]
321clc
322close all
323
324% which seed and chain (of the expt version specified in the INIT cell) to test the plot on
325testSeed = 12;
326testChain = 4;
327DV = 'pitch_stC'; % measure to plot on the y axis, namely pitch expressed either as Hz as before aug21 ('pitch' field), or, as I now think it's better, in semitones above/below middle C ('pitch_stC' field)
328
329noteFqs = [65.41 69.3 73.42 77.78 82.41 87.31 92.5 98 103.83 110 116.54 123.47 130.81 138.59 146.83 155.56 164.81 174.61 185 196 207.65 220 233.08 246.94 261.63 277.18 293.66 311.13 329.63 349.23 369.99 392 415.3 440 466.16 493.88 523.25];
330noteFqs_string = [' 65.41 Hz' ' 69.3 Hz' ' 73.42 Hz' ' 77.78 Hz' ' 82.41 Hz' ' 87.31 Hz' ' 92.5 Hz' ' 98 Hz' '103.83 Hz' '110 Hz' '116.54 Hz' '123.47 Hz' '130.81 Hz' '138.59 Hz' '146.83 Hz' '155.56 Hz' '164.81 Hz' '174.61 Hz' '185 Hz' '196 Hz' '207.65 Hz' '220 Hz' '233.08 Hz' '246.94 Hz' '261.63 Hz' '277.18 Hz' '293.66 Hz' '311.13 Hz' '329.63 Hz' '349.23 Hz' '369.99 Hz' '392 Hz' '415.3 Hz' '440 Hz' '466.16 Hz' '493.88 Hz' '523.25'];
331noteNames = {'C2' 'C#2' 'D2' 'D#2' 'E2' 'F2' 'F#2' 'G2' 'G#2' 'A2' 'A#2' 'B2' 'C3' 'C#3' 'D3' 'D#3' 'E3' 'F3' 'F#3' 'G3' 'G#3' 'A3' 'A#3' 'B3' 'C4' 'C#4' 'D4' 'D#4' 'E4' 'F4' 'F#4' 'G4' 'G#4' 'A4' 'A#4' 'B4' 'C5'};
332
333switch DV
334 case 'pitch'
335 yLabel =[];
336 case 'pitch_stC'
337 yLabel = 'Semitones rel. to C_4';
338end
339
340for i_seed = 1:N_seeds; % 1:N_seeds
341 for i_chain = 1:N_chains % 1:N_chains ; if interested only in a specific seed/chain/gen number, set this as a constant here, otherwise let it loop 1:N_chains
342 switch i_chain
343 case {1,2}
344 chainType = 'mxAcc';
345 case {3,4}
346 chainType = 'hiAcc';
347 end
348 plotName = ['v' vers gender ', seed #' num2str(i_seed) ', chain #' num2str(i_chain) ' (' chainType ')'];
349 figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 600 954], 'Color',[1 1 1])
350
351 % first, plot the seed at the top of the multiplot (subplot) window
352 subplot(N_gens+1, 1, 1)
353 for i_note = 1:N_notesInSeed % plot note by note
354 sec_L = (i_note-1)*seeds(i_seed).noteDuration; % for comparability with chain plots below, we want to label the x axis in terms of seconds even though the seed is plotted in terms of notes; this variable is the x coordinate for the left part of the segment, as expressed in seconds
355 sec_R = i_note*seeds(i_seed).noteDuration;
356 y = seeds(i_seed).(DV)(i_note); % pitch, expressed in Hz ('pitch' field) or in semitones above/below middle C ('pitch_stC' field)
357 plot([sec_L sec_R], [y y] , 'color', colour, 'LineWidth', 2)
358 hold on
359 end
360 set(gca,'xaxisLocation','top') % or just not display labels at all: set(gca, 'XTickLabel', {'' '' '' '' ''});
361 xlabel('Time (s)')
362 ax = gca;
363 switch DV
364 case 'pitch'
365 set(gca, 'YScale', 'log'); % only needed if you are plotting pitches (Hz) and want the ticks to be semitone-apart; if the values are already log2'ed, the scale is to stay linear
366 set(gca, 'YTick', noteFqs);
367 set(gca, 'YTickLabel', noteNames);
368 % ax.YAxis.MinorTick = 'off';
369 ax.YAxis.MinorTickValues = noteFqs;
370 case 'pitch_stC'
371 ax.YAxis.MinorTick = 'on';
372 ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):1:ax.YAxis.Limits(2); % minor y ticks (horizontal lines) at semitone level
373 ax.YAxis.MinorTickValuesMode = 'manual';
374 end
375 set(gca, 'XTick', [0:30]); % high enough value so that the x-axis values don't stop but keep counting to the end
376 box off
377
378 % adjust grid
379 ax.XAxis.MinorTickValues = 0:1:30; % minor x ticks (vertical lines) at every second; 30s = long enough to make sure entire axis is so formatted
380 grid on
381 grid minor
382 set(gca, 'GridLineStyle','-') % solid major grid lines
383 set(gca, 'MinorGridLineStyle','-') % solid minor grid lines
384
385 title('SEED', 'FontSize', 9, 'Units', 'normalized', 'Position', [titleX titleY], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
386% % COD PT A AVEA NOTE LABELS SI IN NOTE NAMES SI IN FQS, NEEDS WORK BEFORE THEN APPLYING TO THE SUBPLOTS BELOW
387% yyaxis left % note names on the left y axis
388% set(gca, 'YScale', 'log');
389% set(gca, 'YTick', noteFqs);
390% set(gca, 'YTickLabel', noteNames);
391% grid on
392% % ax = gca;
393% % ax.YAxis.MinorTickValues = noteFqs;
394% % ax.YAxis.MinorTick = 'off';
395% % limsY=get(gca,'YLim'); set(gca,'Ylim',[0.95*limsY(1) 1.1*limsY(2)]);
396% yyaxis right % note fqs on the right y axis
397% set(gca, 'YScale', 'log');
398% set(gca, 'YTick', noteFqs);
399% set(gca, 'YTickLabel', noteFqs_string);
400% grid on
401
402 % now plot each generation on its own subplot row
403 for i_gen = 1:N_gens % if interested only in a specific seed/chain/gen number, set this as a constant here, otherwise let it loop 1:N_...
404 subplot(N_gens+1, 1, i_gen+1)
405 N = nnz(notes(i_seed, i_chain).t1(i_gen,:));
406 for i_note = 1:N % these are in fact the horizontal lines in the Praat plot, which of course don't *always* in fact correspond to actual notes
407 t1 = notes(i_seed, i_chain).t1(i_gen, i_note); % these vars created for ease of reading the plot command
408 t2 = notes(i_seed, i_chain).t2(i_gen, i_note);
409 y = notes(i_seed, i_chain).(DV)(i_gen, i_note);
410 plot([t1 t2], [y y] , 'color', colour, 'LineWidth', 2) % to only plot those lines that exceed the min duration, we rely on _highpassing the time and pitch variables
411 box off % avoids repeating y tickmarks on the right hand side of the axes
412 hold on
413 end
414
415 limsX=get(gca,'XLim');
416 set(gca,'Xlim',[ 0 limsX(2)]); % 2nd value in vector is max value on x axis; so that the time axis always starts at 0, and not just when the first note starts
417 limsY=get(gca,'YLim');
418 set(gca,'Ylim',[0.95*limsY(1) 1.1*limsY(2)]); % y padding; otherwise the top-&bottom-most notes are right on the top&bottom border of plot, and look as if they dont belong to the plot
419 title(['Gen.' num2str(i_gen)], 'FontSize', 9, 'Units', 'normalized', 'Position', [titleX titleY], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
420 set(gca, 'XTick', [0:30]); % high enough value so that the x-axis values don't stop but keep counting to the end
421 if i_gen == N_gens % apply labels only to the x axis of the final subplot line
422 xlabel('Time (s)')
423 ylabel(yLabel)
424 else
425 set(gca, 'XTickLabel', {'' '' '' '' ''});
426% set(gca,'xtick',[])
427 end
428
429 % only if we plot pitch (Hz) values, not semitones above middle C
430 switch DV
431 case 'pitch'
432 % label y axis in terms of semitones
433 set(gca, 'YTick', noteFqs);
434 set(gca, 'YTickLabel', {'C2' 'C#2' 'D2' 'D#2' 'E2' 'F2' 'F#2' 'G2' 'G#2' 'A2' 'A#2' 'B2' 'C3' 'C#3' 'D3' 'D#3' 'E3' 'F3' 'F#3' 'G3' 'G#3' 'A3' 'A#3' 'B3' 'C4' 'C#4' 'D4' 'D#4' 'E4' 'F4' 'F#4' 'G4' 'G#4' 'A4' 'A#4' 'B4' 'C5'})
435 set(gca, 'YScale', 'log');
436
437 ax = gca;
438 ax.YAxis.MinorTickValues = noteFqs;
439 case 'pitch_stC'
440 ax = gca;
441 ax.YAxis.MinorTick = 'on';
442 ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):1:ax.YAxis.Limits(2); % minor y ticks (horizontal lines) at semitone level
443 end
444
445 % add grid
446 grid on
447 grid minor
448 set(gca, 'GridLineStyle','-') % solid major grid lines
449 set(gca, 'MinorGridLineStyle','-') % solid minor grid lines
450 ax.XAxis.MinorTickValues = 0:1:30; % minor x ticks (vertical lines) at every second; 30s = long enough to make sure entire axis is so formatted
451 end
452
453 % use same x ranges across all subplots
454 ax_handles=get(gcf,'children');
455 if max(hums(i_seed, i_chain).duration(1:N_gens)) > duration_seed % duration_seed is set at the very beginning, as a fct of expt version; i_chain set by for loop, whereas i_seed is fixed at the top of the script
456 xMax = max(hums(i_seed, i_chain).duration(1:N_gens));
457 else
458 xMax = duration_seed;
459 end
460 set(ax_handles, 'XLim', [0, xMax]);
461
462 if strcmp(DV, 'pitch')==1
463 % if the range is too high, then make y labels less dense (plot only (C,G) or (C,E,G) from every octave), so that labels dont overlap
464 r = max(allYLim) / min(allYLim); % frequency range, to be transformed into semitone difference
465 range_St = 12 * log2(r);
466 if range_St > maxStOnY
467 for i_subplot = 1:N_gens+1
468 subplot(N_gens+1, 1, i_subplot)
469 set(gca, 'YTick', [65.41 82.41 98 130.81 164.81 196 261.63 329.63 392 523.25]);
470 set(gca, 'YTickLabel', {'C2' 'E2' 'G2' 'C3' 'E3' 'G3' 'C4' 'E4' 'G4' 'C5'})
471% set(gca, 'YTick', [65.41 98 130.81 196 261.63 392 523.25]);
472% set(gca, 'YTickLabel', {'C2' 'G2' 'C3' 'G3' 'C4' 'G4' 'C5'})
473 end
474 end
475
476 % use same y ranges across all subplots, as a function of the widest existing range (alltime min, alltime max)
477 ax_handles=get(gcf,'children');
478 allYLim = get(ax_handles, {'YLim'});
479 allYLim = cat(2, allYLim{:});
480 set(ax_handles, 'YLim', [min(allYLim), max(allYLim)]);
481 end
482
483 h_suptitle = suptitle(plotName, 0.98); % this command, if called earlier, ruins the title's horizontal alignment
484 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
485
486 % save figure as raster image and as FIG
487 print('-dpng', ['d:\_matlabPlots\' plotName ', DV=' DV '.png' ] );
488% savefig(['d:\' plotName '.fig' ])
489 end
490end
491disp('gata...........................................................................................................................................................................')
492
493%% 5.1.Tony final-note analyses [[partea de plot fa-o cu un cell P doar dc reiese ca are de fapt sens acest tip de analiza]]
494clc
495close all
496
497testSeed = 3;
498% there's no testChain, since below we pool across chains at some point
499
500for i_seed = 1:N_seeds % 1:N_seeds, or testSeed
501 plotName = ['Final notes from chains of seed #' num2str(i_seed) ', expt ' vers gender ', pianorolls from Tony'];
502 figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 885 380], 'Color',[1 1 1])
503
504 for i_chain = 1:N_chains
505 switch i_chain
506 case {1,2}
507 chainType = 'mxAcc';
508 case {3,4}
509 chainType = 'hiAcc';
510 end
511
512 % plot current chain's distribution (we pool across all 8 gens)
513 subplot(N_chains+1, 1, i_chain)
514 plot(squeeze(hums(i_seed, i_chain).lastPitch(:))', 1, 'ko') % the last notes of all 8 gens
515 hold on
516 plot( seeds(i_seed).pitch(end), 1, 'o', 'MarkerFaceColor', 'b') % the last pitch of the seed, highlighted
517 plot(2*seeds(i_seed).pitch(end), 1, 'o', 'MarkerFaceColor', [151 151 255]./255) % ..and its version one octave higher
518 box off % avoids repeating y tickmarks on the right hand side of the axes
519 ax1 = gca; % gca = get current axis
520 ax1.YAxis.Visible = 'off'; % remove y-axis
521 pbaspect([1 0.001 1]) %squash the plot vertically
522 title(['Chain#' num2str(i_chain) ' (' chainType ')'], 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
523 end
524
525 % last row: pool across all chains
526 subplot(N_chains+1, 1, N_chains+1)
527 plot(allChains(i_seed).lastPitch, 1, 'ko') % the last notes of the 8 generations
528 hold on
529 plot( seeds(i_seed).pitch(end), 1, 'o', 'MarkerFaceColor', 'b') % just like above
530 plot(2*seeds(i_seed).pitch(end), 1, 'o', 'MarkerFaceColor', [151 151 255]./255)
531 box off % avoids repeating y tickmarks on the right hand side of the axes
532 ax1 = gca; % gca = get current axis
533 ax1.YAxis.Visible = 'off'; % remove y-axis
534 pbaspect([1 0.001 1]) %squash the plot vertically
535 title('ALL', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
536 xlabel('Hz')
537
538 % uniform xlims for all subplots
539 AX_handles=get(gcf,'children');
540 allXLim = get(AX_handles, {'XLim'});
541 allXLim = cat(2, allXLim{:});
542 set(AX_handles, 'XLim', [min(allXLim), max(allXLim)]);
543
544 h_suptitle = suptitle(plotName, 0.95); % as the plot's title (suptitle), you may want to use a shortened version of the plotName; this command, if called earlier, ruins the title's horizontal alignment
545 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
546
547 % save figure as raster image and as FIG
548% set(gcf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
549% print(gcf, '-painters', '-dpng', ['d:\' plotName '.png','-r0' ] );
550% savefig(['d:\' plotName '.fig' ])
551 print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] ); % save figure
552end
553disp('gata...........................................................................................................................................................................')
554
555%% 5.1.Prosogram final-note analyses [[partea de plot fa-o cu un cell P doar dc reiese ca are de fapt sens acest tip de analiza]]
556clc
557close all
558
559testSeed = 3; % momentan tb sa fol acelasi seed pt care ai calculat lastPitch_seed in cell 2.1 :(
560% there's no testChain, since below we pool across chains at some point
561
562tMin_lastNote = 1; % how long should the final note be (in s) for us to consider it real (not an artefact)? for Prosogram piano rolls, 1s is good, but with Tony ones it doesn't make sense to use this variable, so just set it to 0 so the condition is always met
563lastPitch = NaN(N_seeds, N_chains, N_gens);
564
565for i_seed = testSeed % 1:N_seeds
566
567 plotName = ['Final notes from chains of seed #' num2str(i_seed) ', expt ' vers gender ', pianorolls from Prosogram'];
568 figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 885 380], 'Color',[1 1 1])
569 lastPitch_allChains = []; % for current seed
570
571 for i_chain = 1:N_chains
572 switch i_chain
573 case {1,2}
574 chainType = 'mxAcc';
575 case {3,4}
576 chainType = 'hiAcc';
577 end
578
579 % derive last pitches
580 for i_gen=1:N_gens
581 lastNonNan = find(isnan(pitches_highpass(i_seed, i_chain, i_gen, :)));
582 lastNonNan = min(lastNonNan)-1; % index of the last note in the vector; next one is already a NaN
583 if noteDurations_highpass(i_seed, i_chain, i_gen, lastNonNan) > tMin_lastNote % if condition not met, then for this i_gen, lastPitch skips an element
584 lastPitch(i_seed, i_chain, i_gen) = pitches_highpass(i_seed, i_chain, i_gen, lastNonNan);
585 end
586 end
587 lastPitch_allChains = [ lastPitch_allChains squeeze(lastPitch(i_seed, i_chain, :))' ]; % add values from current chain
588
589 % plot current chain's distribution (we pool across all gens)
590 subplot(N_chains+1, 1, i_chain)
591 plot(squeeze(lastPitch(i_seed, i_chain, :))', 1, 'ko') % the last notes of the 8 generations
592 hold on
593 plot(lastPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', 'b') % the last of the seed, highlighted. ASTA E CALCULAT IN CELL 2.1.x, SUPER PROST INSA CODAT, TB ADAPTAT MASIV
594 plot(2*lastPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', [151 151 255]./255) % ..and its version one octave higher
595 box off % avoids repeating y tickmarks on the right hand side of the axes
596 ax1 = gca; % gca = get current axis
597 ax1.YAxis.Visible = 'off'; % remove y-axis
598 pbaspect([1 0.001 1]) %squash the plot vertically
599 title(['Chain#' num2str(i_chain) ' (' chainType ')'], 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
600 end
601
602 % last row: pool across all chains
603 subplot(N_chains+1, 1, N_chains+1)
604 plot(lastPitch_allChains, 1, 'ko') % the last notes of the 8 generations
605 hold on
606 plot(lastPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', 'b') % the last of the seed, highlighted
607 plot(2*lastPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', [151 151 255]./255) % ..and its version one octave higher
608 box off % avoids repeating y tickmarks on the right hand side of the axes
609 ax1 = gca; % gca = get current axis
610 ax1.YAxis.Visible = 'off'; % remove y-axis
611 pbaspect([1 0.001 1]) %squash the plot vertically
612 title('ALL', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
613 xlabel('Hz')
614
615 % uniform xlims for all subplots
616 AX_handles=get(gcf,'children');
617 allXLim = get(AX_handles, {'XLim'});
618 allXLim = cat(2, allXLim{:});
619 set(AX_handles, 'XLim', [min(allXLim), max(allXLim)]);
620
621 h_suptitle = suptitle(plotName, 0.95); % as the plot's title (suptitle), you may want to use a shortened version of the plotName; this command, if called earlier, ruins the title's horizontal alignment
622 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
623
624 % save figure as raster image and as FIG
625% set(gcf,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
626% print(gcf, '-painters', '-dpng', ['d:\' plotName '.png','-r0' ] );
627% savefig(['d:\' plotName '.fig' ])
628 print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] ); % save figure
629end
630disp('gata...........................................................................................................................................................................')
631
632%% 5.2.Tony first-note analyses [[partea de plot fa-o cu un cell P doar dc reiese ca are de fapt sens acest tip de analiza]]
633clc
634close all
635
636testSeed = 3; % there's no testChain, since below we pool across chains at some point
637
638for i_seed = 1:N_seeds % 1:N_seeds, or testSeed
639 plotName = ['First notes from chains of seed #' num2str(i_seed) ', expt ' vers gender ', pianorolls from Tony'];
640 figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 885 380], 'Color',[1 1 1])
641
642 for i_chain = 1:N_chains
643 switch i_chain
644 case {1,2}
645 chainType = 'mxAcc';
646 case {3,4}
647 chainType = 'hiAcc';
648 end
649
650 % plot current chain's distribution (we pool across all 8 gens)
651 subplot(N_chains+1, 1, i_chain)
652 plot(squeeze(hums(i_seed, i_chain).firstPitch(:))', 1, 'ko') % the first notes of all 8 gens
653 hold on
654 plot( seeds(i_seed).pitch(1), 1, 'o', 'MarkerFaceColor', 'b') % the first pitch of the seed, highlighted
655 plot(2*seeds(i_seed).pitch(1), 1, 'o', 'MarkerFaceColor', [151 151 255]./255) % ..and its version one octave higher
656 box off % avoids repeating y tickmarks on the right hand side of the axes
657 ax1 = gca; % gca = get current axis
658 ax1.YAxis.Visible = 'off'; % remove y-axis
659 pbaspect([1 0.001 1]) %squash the plot vertically
660 title(['Chain#' num2str(i_chain) ' (' chainType ')'], 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
661 end
662
663 % last row: pool across all chains
664 subplot(N_chains+1, 1, N_chains+1)
665 plot(allChains(i_seed).firstPitch, 1, 'ko')
666 hold on
667 plot( seeds(i_seed).pitch(1), 1, 'o', 'MarkerFaceColor', 'b') % just like above
668 plot(2*seeds(i_seed).pitch(1), 1, 'o', 'MarkerFaceColor', [151 151 255]./255)
669 box off % avoids repeating y tickmarks on the right hand side of the axes
670 ax1 = gca; % gca = get current axis
671 ax1.YAxis.Visible = 'off'; % remove y-axis
672 pbaspect([1 0.001 1]) %squash the plot vertically
673 title('ALL', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
674 xlabel('Hz')
675
676 % uniform xlims for all subplots
677 AX_handles=get(gcf,'children');
678 allXLim = get(AX_handles, {'XLim'});
679 allXLim = cat(2, allXLim{:});
680 set(AX_handles, 'XLim', [min(allXLim), max(allXLim)]);
681
682 h_suptitle = suptitle(plotName, 0.95); % as the plot's title (suptitle), you may want to use a shortened version of the plotName; this command, if called earlier, ruins the title's horizontal alignment
683 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
684
685 % save figure as raster image and as FIG
686 print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] ); % save figure
687end
688disp('gata...........................................................................................................................................................................')
689
690%% 5.2.Prosogram first-note analyses
691clc
692close all
693% load 'd:\Experiments\analyses\saved variables\pitches_highpass.mat'
694
695tMin_firstNote = 0; % how long should the note be (in s) for us to consider it real (not an artefact)?for Prosogram piano rolls, used to have this at 1s, but with Tony ones it doesn't make sense to use this variable, so just set it to 0 so the condition is always met
696firstPitch = NaN(N_seeds, N_chains, N_gens);
697
698for i_seed = testSeed % 1:N_seeds
699
700 plotName = ['First notes, chains of seed #' num2str(i_seed) ', expt v' vers ' ' gender];
701 figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 660 380], 'Color',[1 1 1])
702 firstPitch_allChains = []; % for current seed
703
704 for i_chain = 1:N_chains
705 switch i_chain
706 case {1,2}
707 chainType = 'mxAcc';
708 case {3,4}
709 chainType = 'hiAcc';
710 end
711
712 % derive first pitches
713 for i_gen=1:N_gens
714 if noteDurations_highpass(i_seed, i_chain, i_gen, 1) > tMin_firstNote % if condition not met, then for this i_gen, firstPitch skips an element
715 firstPitch(i_seed, i_chain, i_gen) = pitches_highpass(i_seed, i_chain, i_gen, 1);
716 end
717 end
718 firstPitch_allChains = [ firstPitch_allChains squeeze(firstPitch(i_seed, i_chain, :))' ]; % add values from current chain
719
720 % plot current chain's distribution
721 subplot(N_chains+1, 1, i_chain)
722 plot(squeeze(firstPitch(i_seed, i_chain, :))', 1, 'ko') % the first notes of the 8 generations
723 hold on
724 plot(firstPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', 'b') % the first of the seed, highlighted
725 plot(2*firstPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', [151 151 255]./255) % ..and its version one octave higher
726 box off % avoids repeating y tickmarks on the right hand side of the axes
727 ax1 = gca; % gca = get current axis
728 ax1.YAxis.Visible = 'off'; % remove y-axis
729 pbaspect([1 0.001 1]) %squash the plot vertically
730 title(['Chain#' num2str(i_chain) ' (' chainType ')'], 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
731 end
732
733 % last row: pool across all chains
734 subplot(N_chains+1, 1, N_chains+1)
735 plot(firstPitch_allChains, 1, 'ko') % the first notes of the 8 generations
736 hold on
737 plot(firstPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', 'b') % the first of the seed, highlighted
738 plot(2*firstPitch_seed(i_seed), 1, 'o', 'MarkerFaceColor', [151 151 255]./255) % ..and its version one octave higher
739 box off % avoids repeating y tickmarks on the right hand side of the axes
740 ax1 = gca; % gca = get current axis
741 ax1.YAxis.Visible = 'off'; % remove y-axis
742 pbaspect([1 0.001 1]) %squash the plot vertically
743 title('ALL', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 -5], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5], 'FontWeight','bold')
744 xlabel('Hz')
745
746
747 % uniform xlims for all subplots
748 AX_handles=get(gcf,'children');
749 allXLim = get(AX_handles, {'XLim'});
750 allXLim = cat(2, allXLim{:});
751 set(AX_handles, 'XLim', [min(allXLim), max(allXLim)]);
752
753 h_suptitle = suptitle(plotName, 0.95); % as the plot's title (suptitle), you may want to use a shortened version of the plotName; this command, if called earlier, ruins the title's horizontal alignment
754 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
755 % save figure as raster image and as FIG
756 print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] );
757 savefig(['d:\' plotName '.fig' ])
758end
759disp('gata...........................................................................................................................................................................')
760
761%% 6.1 histograms for pitch-class -> airtime estimates
762% has yet to be generalised to work with DVs other than PC/airtime
763clc
764close all
765clear binRanges binSteps
766
767binSteps.PC = 1; % semitones
768binRanges.PC = binSteps.PC:binSteps.PC:12; % bcs there are 12 pitch classes; this is for the 2nd input arg of histc, called binranges. SHould really start from 1, dar ca sa fie acelasi nr de elem lasa-l sa inceapa cu steps
769
770testSeed = 4;
771testChain = 1;
772
773N_plots = N_gens+3;
774PCs = {'C' 'C#' 'D' 'D#' 'E' 'F' 'F#' 'G' 'G#' 'A' 'A#' 'B'};
775reorder = [1 8 3 10 5 12 7 2 9 4 11 6]; % from circle of fifths, or equivalently from Fabian's thesis Fig 10.2. Gets adjusted to binSteps below
776
777for i_seed = 1:N_seeds; % 1:N_seeds
778 for i_chain = 1:N_chains % testChain or 1:N_chains
779 switch i_chain
780 case {1,2}
781 chainType = 'mxAcc';
782 case {3,4}
783 chainType = 'hiAcc';
784 end
785 plotName = ['Airtime (pitch classes), binSteps=' num2str(binSteps.PC) ' semitones, seed #' num2str(i_seed) ', chain #' num2str(i_chain) ' (' chainType '), expt v' vers gender];
786 figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 900 954], 'Color',[1 1 1])
787
788 % first, plot the histo for the seed at the top of the multiplot (subplot) window
789 [binCounts] = histc(seeds(i_seed).PC(:), binRanges.PC);
790 subplot(N_plots, 1, 1)
791 h = bar(binRanges.PC, binCounts, 'histc'); % 'histc' makes the bar plot in the style of a histogram
792 h.FaceColor = [.5 .5 .5];
793 set(gca,'xaxisLocation','top') % or just not display labels at all: set(gca, 'XTickLabel', {'' '' '' '' ''});
794 set(gca, 'XTick', [1.5:12.5]);
795 set(gca, 'XTickLabel', PCs)
796 title('SEED', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5])
797 xlim([0.9,13.5])
798
799 % now plot each generation on its own subplot row
800 for i_gen = 1:N_gens % if interested only in a specific seed/chain/gen number, set this as a constant here, otherwise let it loop 1:N_...
801 [binCounts] = histc(hums(i_seed, i_chain).pcAirtime(i_gen, :), binRanges.PC); % the 2nd input arg is called binranges in the help, although i find that term confusing!
802 binCounts = binCounts / 1000; % make airtime values from ms into seconds
803 subplot(N_plots, 1, i_gen+1)
804 bar(binRanges.PC, binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
805
806% h = findobj(gca,'Type','line'); set(h,'Marker','none'); % remove the asterisks/stars that the bar function for some reason places on the x axis
807 box off % avoids repeating y tickmarks on the right hand side of the axes
808 title(['Gen.' num2str(i_gen)], 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight', 'normal')
809 set(gca, 'XTickLabel', {'' '' '' '' ''});
810
811 h_suptitle = suptitle(plotName, 0.95); % as the plot's title (suptitle), you may want to use a shortened version of the plotName; this command, if called earlier, ruins the title's horizontal alignment
812 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
813 xlim([0.9,13.5])
814 end
815
816 % now plot pooling across all generations
817 [binCounts] = histc(reshape(hums(i_seed, i_chain).pcAirtime(:, :), 1, []), binRanges.PC); % pool from all across the hum, across all gens
818 binCounts = binCounts / 1000; % make airtime values from ms into seconds
819 subplot(N_plots, 1, N_gens+2)
820 bar(binRanges.PC, binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
821 xlabel('Pitch class')
822 ylabel('Seconds')
823 xlim([0.9,13.5])
824 set(gca, 'XTick', [1.5:12.5]);
825 set(gca, 'XTickLabel', PCs)
826 box off
827 title({'ALL GENS' '(monotonous)'}, 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight','bold')
828
829 % same plot but with pitch class on the x axis going up in fifths (7-semitone-steps) rather than in 1-semitone-steps
830 % first, adjust the reorder vector to fit whatever subdivision of the semitones is given by the chosen binStep
831 % if this step is 1 (default), the vector should stay the same
832 % ACTUALLY, THIS DOESNT GIVE THE EXPECTED RESULTS FOR STEPS DIFFERENT FROM 1, SO THIS CODE NEEDS MORE WORK!
833 reorderNew = [];
834 for i_PC=12:-1:1
835 repetitions(i_PC, 1:1/binSteps.PC) = repmat(reorder(i_PC), [1 1/binSteps.PC]);
836 reorderNew = [reorderNew repetitions(i_PC, 1:1/binSteps.PC)];
837 end
838 reorder = flip(reorderNew);
839 % now continue to make the plot
840 binCounts = binCounts(reorder);
841 subplot(N_plots, 1, N_plots)
842 bar(binRanges.PC, binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
843 xlabel('Pitch class')
844 ylabel('Seconds')
845 xlim([0.9,13.5])
846 set(gca, 'XTick', [1.5:12.5]);
847 set(gca, 'XTickLabel', PCs(reorder)) % x-axis labels reordered too
848 box off
849 title({'ALL GENS' '(in fifths)'}, 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight','bold')
850
851 % save figure as raster image and as FIG
852 print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] );
853% savefig(['d:\' plotName '.fig' ])
854 end
855end
856disp('gata...........................................................................................................................................................................')
857
858%% 6.2 - as 6.1, but pool across all of the testSeed's chains; you don't need to've run 6.1 before
859% modularised to do the work that 7.2 would have done
860clc
861close all
862
863testSeed = 2;
864DV = 'interval_st'; %'PC' or 'interval_st'
865
866switch DV
867 case 'PC'
868 reorder = [1 8 3 10 5 12 7 2 9 4 11 6]; % from circle of fifths, or equivalently from Fabian's thesis Fig 10.2
869 PCs = {'C' 'C#' 'D' 'D#' 'E' 'F' 'F#' 'G' 'G#' 'A' 'A#' 'B'};
870 N_plots = 3; % seed, allChains_monotonous, allChains_inFifths
871 plotName = ['Airtime (pitch classes), all chains of seed #' num2str(testSeed) ', expt v' vers gender];
872 case 'interval_st'
873 N_plots = 2; % seed, allChains
874 plotName = ['Interval distribution, all chains of seed #' num2str(testSeed) ', expt v' vers gender];
875end
876
877figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 972 954], 'Color',[1 1 1])
878h_suptitle = suptitle(plotName, 0.97);
879set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
880
881% first, plot the histo for the seed at the top of the multiplot (subplot) window
882[binCounts] = histc(seeds(testSeed).(DV)(:), binRanges.(DV));
883subplot(N_plots, 1, 1)
884h = bar(binRanges.(DV), binCounts, 'histc'); % 'histc' makes the bar plot in the style of a histogram
885h.FaceColor = [.5 .5 .5];
886% set(gca,'xaxisLocation','top') % or just not display labels at all: set(gca, 'XTickLabel', {'' '' '' '' ''});
887title('SEED', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5])
888switch DV
889 case 'PC'
890 xlim([0.9,13.5])
891 set(gca, 'XTick', [1.5:12.5]);
892 set(gca, 'XTickLabel', PCs)
893 case 'interval_st'
894 xlim([-12, 12])
895end
896
897% plot pooling across all generations AND all chains
898switch DV % can't use dynamic field names as I have both 'pcAirtime' for when it's repeated and 'PC' when it's just pitch classes
899 case 'PC'
900 [binCounts] = histc(allChains(testSeed).pcAirtime, binRanges.(DV));
901 binCounts = binCounts / 1000; % make airtime values from ms into seconds
902 case 'interval_st'
903 [binCounts] = histc(allChains(testSeed).(DV), binRanges.(DV));
904end
905subplot(N_plots, 1, 2)
906bar(binRanges.(DV), binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
907switch DV % need to have another switch, post-plot-command
908 case 'PC'
909 xlabel('Pitch class')
910 ylabel('Seconds')
911 xlim([0.9,13.5])
912 set(gca, 'XTick', [1.5:12.5]);
913 set(gca, 'XTickLabel', PCs)
914 case 'interval_st'
915 xlabel('Interval (semitones)')
916 ylabel('Freq of occurrence')
917 xlim([-12 12])
918end
919box off
920title({'ALL CHAINS' 'ALL GENS'}, 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight','bold')
921
922% for PCs, same plot but with pitch class on the x axis going up in fifths (7-semitone-steps) rather than in 1-semitone-steps
923if strcmp(DV, 'PC')==1
924 binCounts = binCounts(reorder);
925 subplot(N_plots, 1, N_plots)
926 bar(binRanges.PC, binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
927 xlabel('Pitch class')
928 ylabel('Seconds')
929 xlim([0.9,13.5])
930 set(gca, 'XTick', [1.5:12.5]);
931 set(gca, 'XTickLabel', PCs(reorder)) % x-axis labels reordered too
932 box off
933 title({'ALL CHAINS' 'ALL GENS' '(in fifths)'}, 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight','bold')
934end
935
936% save figure as raster image and as FIG
937print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] );
938disp('gata...........................................................................................................................................................................')
939
940%% 6.3 - as 6.1, but pool across both chains and seeds
941clc
942close all
943clear binRanges binSteps
944
945DV = 'PC'; %'PC' or 'interval_st'
946
947binSteps.interval_st = 0.1; % step determines the 1D histogram's resolution: 1 = semitone, 0.01 = cent
948binSteps.PC = 1;
949binRanges.PC = 1:binSteps.PC:12; % bcs there are 12 pitch classes; this is for the 2nd input arg of histc, called binranges in the help, although i find that term confusing!
950binRanges.interval_st = -12:binSteps.interval_st:12;
951
952switch DV
953 case 'PC'
954 reorder = [1 8 3 10 5 12 7 2 9 4 11 6]; % from circle of fifths, or equivalently from Fabian's thesis Fig 10.2
955 PCs = {'C' 'C#' 'D' 'D#' 'E' 'F' 'F#' 'G' 'G#' 'A' 'A#' 'B'};
956 N_plots = 2; % allChains_monotonous, allChains_inFifths
957 plotName = ['Airtime (pitch classes), all chains and all seeds, expt v' vers gender];
958 [binCounts] = histc(allChainsAndSeeds.pcAirtime, binRanges.(DV)); % can't use dynamic field names here, as I have both 'pcAirtime' for when it's repeated and 'PC' when it's just pitch classes
959 binCounts = binCounts / 1000; % make airtime values from ms into seconds
960 case 'interval_st'
961 N_plots = 1; % allChains
962 plotName = ['Interval distribution, all chains and all seeds, expt v' vers gender];
963 [binCounts] = histc(allChainsAndSeeds.(DV), binRanges.(DV));
964
965 % Assess symmetry of rising and falling intervals
966 [N_risingIntervals] = length( find(allChainsAndSeeds.(DV) > 0 ) );
967 [N_fallingIntervals] = length( find(allChainsAndSeeds.(DV) < 0 ) );
968 [N_intervals] = sum ( ~isnan(allChainsAndSeeds.(DV) ) );
969 disp(['We have ' num2str(N_risingIntervals) ' rising intervals (' num2str(100*N_risingIntervals/N_intervals, '%2.1f') '%), and ' num2str(N_fallingIntervals) ' falling intervals (' num2str(100*N_fallingIntervals/N_intervals, '%2.1f') '%)'])
970end
971
972% plot pooling across all generations AND all chains
973figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 972 954], 'Color',[1 1 1])
974h_suptitle = suptitle(plotName, 0.97);
975set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
976subplot(N_plots, 1, 1)
977bar(binRanges.(DV), binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
978switch DV % need to have another switch, post-plot-command
979 case 'PC'
980 xlabel('Pitch class')
981 ylabel('Seconds')
982 xlim([0.9,13.5])
983 set(gca, 'XTick', [1.5:12.5]);
984 set(gca, 'XTickLabel', PCs)
985 case 'interval_st'
986 xlabel(['Interval (semitones), bin step = ' num2str(binSteps.interval_st) ' semitones'])
987 ylabel('Frequency of occurrence')
988 xlim([-13 13])
989 set(gca, 'Xtick', -12:12)
990end
991box off
992title({'ALL SEEDS' 'ALL CHAINS' 'ALL GENS'}, 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight','bold')
993
994% for PCs, same plot but with pitch class on the x axis going up in fifths (7-semitone-steps) rather than in 1-semitone-steps
995if strcmp(DV, 'PC')==1
996 binCounts = binCounts(reorder);
997 subplot(N_plots, 1, N_plots)
998 bar(binRanges.PC, binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
999 xlabel('Pitch class')
1000 ylabel('Seconds')
1001 xlim([0.9,13.5])
1002 set(gca, 'XTick', [1.5:12.5]);
1003 set(gca, 'XTickLabel', PCs(reorder)) % x-axis labels reordered too
1004 box off
1005 title({'ALL SEEDS' 'ALL CHAINS' 'ALL GENS' '(in fifths)'}, 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight','bold')
1006end
1007
1008% save figure as raster image and as FIG
1009print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] );
1010
1011disp('gata...........................................................................................................................................................................')
1012
1013%% 7.1 histograms for intervals, seed-wise
1014% doesnt make sense to modularise this cell together with 6.1, to plot histograms of both airtime and intervals, since too much of the code is different, more trbl than it s wrth
1015% actually, maybe I should modularise it, as I did with 6.2!
1016clc
1017close all
1018
1019testSeed = 3;
1020testChain = 1;
1021
1022N_plots = N_gens+2;
1023
1024for i_seed = testSeed; % 1:N_seeds
1025 for i_chain = 1:N_chains % testChain or 1:N_chains
1026 switch i_chain
1027 case {1,2}
1028 chainType = 'mxAcc';
1029 case {3,4}
1030 chainType = 'hiAcc';
1031 end
1032 plotName = ['Intervals, seed #' num2str(i_seed) ', chain #' num2str(i_chain) ' (' chainType '), expt v' vers gender];
1033 figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 900 954], 'Color',[1 1 1])
1034
1035 % first, plot the histo for the seed at the top of the multiplot (subplot) window
1036 [binCounts] = histc(seeds(i_seed).interval_st(:), binRanges.interval_st);
1037 subplot(N_plots, 1, 1)
1038 h = bar(binRanges.interval_st, binCounts, 'histc'); % 'histc' makes the bar plot in the style of a histogram
1039 h.FaceColor = [.5 .5 .5];
1040 set(gca,'xaxisLocation','top') % or just not display labels at all: set(gca, 'XTickLabel', {'' '' '' '' ''});
1041 set(gca, 'XTick', [-12:2:12]);
1042% set(gca, 'XTickLabel', PCs)
1043 title('SEED', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'Color', [.5 .5 .5])
1044 xlim([-12 12])
1045
1046 % now plot each generation on its own subplot row
1047 for i_gen = 1:N_gens % if interested only in a specific seed/chain/gen number, set this as a constant here, otherwise let it loop 1:N_...
1048 [binCounts] = histc(hums(i_seed, i_chain).interval_st(i_gen, :), binRanges.interval_st); % the 2nd input arg is called binranges in the help, although i find that term confusing!
1049 subplot(N_plots, 1, i_gen+1)
1050 bar(binRanges.interval_st, binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
1051
1052 box off % avoids repeating y tickmarks on the right hand side of the axes
1053 title(['Gen.' num2str(i_gen)], 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight', 'normal')
1054 set(gca, 'XTickLabel', {'' '' '' '' ''});
1055
1056 h_suptitle = suptitle(plotName, 0.95); % as the plot's title (suptitle), you may want to use a shortened version of the plotName; this command, if called earlier, ruins the title's horizontal alignment
1057 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
1058 xlim([-12 12])
1059 end
1060
1061 % now plot pooling across all generations
1062 [binCounts] = histc(reshape(hums(i_seed, i_chain).interval_st(:, :), 1, []), binRanges.interval_st); % pool from all across the hum, across all gens
1063 subplot(N_plots, 1, N_gens+2)
1064 bar(binRanges.interval_st, binCounts, 'histc') % 'histc' makes the bar plot in the style of a histogram
1065 xlabel('Interval (semitones)')
1066 ylabel('Freq of occurrence')
1067 xlim([-12 12])
1068 set(gca, 'XTick', [-12:2:12]);
1069% set(gca, 'XTickLabel', PCs)
1070 box off
1071 title('ALL GENS', 'FontSize', 9, 'Units', 'normalized', 'Position', [-.16 .37], 'HorizontalAlignment', 'left', 'FontWeight','bold')
1072
1073 % save figure as raster image and as FIG
1074 print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] );
1075 end
1076end
1077disp('gata...........................................................................................................................................................................')
1078
1079%% 7.2 here I actually did modularise 6.2 ,so use that :)
1080
1081%% 7.3 2D histograms for pairs (bigrams) of notes or intervals
1082clc
1083close all
1084
1085poolAcross = 'seeds'; % 'chains', 'seeds' or 'chains and seeds'
1086DV = 'notes'; % bigrams of 'notes' or of 'intervals'?
1087
1088fontSize_labels = 12;
1089fontSize_values = 10;
1090fontSize_legendLabel = 11;
1091
1092switch poolAcross % if we decided to pool across chains, the 'for' below loops across chains, and vice-versa
1093 case 'chains'
1094 N = N_seeds;
1095 N_plots_x = 3;
1096 N_plots_y = 4;
1097 figX = 1697;
1098 figY = 963;
1099 case 'seeds'
1100 N = N_chains;
1101 N_plots_x = 2;
1102 N_plots_y = 2;
1103 figX = 1165;
1104 figY = 954;
1105end
1106
1107switch DV
1108 case 'intervals'
1109 aboveMiddleC = [];
1110 case 'notes'
1111 aboveMiddleC = ' rel. to C_4';
1112end
1113
1114plotName = ['2D histograms of bigrams (pairs) of ' DV ', pooling across ' poolAcross ', expt v' vers gender];
1115figure('Name', plotName, 'NumberTitle','off', 'Position', [0 44 figX figY], 'Color',[1 1 1])
1116h_suptitle = suptitle(plotName, 0.97);
1117set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
1118
1119% % find out the max number of the colorbar
1120% counts_max = 1; % initialise
1121% for i_genre=1:N_genre
1122% pls = ratings(:, i_genre, :, :, 1); % pleasantness ratings, values from across all ss, exemplars and repetitions, i.e. x matrix has dim1=ss, dim2=exemplars, dim3=reps; size(x), just like size(y), is 30x1x20x2
1123% rgh = ratings(:, i_genre, :, :, 2); % same but for roughness (these are flipped values, i.e. 8-..); positions across the x and y matrices correspond, e.g. x(3,4,2) and y(3,4,2) represent the pls and rgh ratings respectively for subj3, exemplar4, repetition2, given the current i_genre; it makes no sense to average across reps at this stage
1124% pls = reshape(pls,1,[]); % since at this stage we don't care about the subject, exemplar and repetition identity of each rating (data point) but instead want them to all look/contribute the same, we therefore rearrange the values in the 3D matrices (x and y) to a single row vector (1x1200), for easier scatter-plotting; the correspondence between elements of x and y of equal index still remains - a sanity check for that is to look at the two vectors side by side and observe that the NaNs always occur for both ratings simultaneously, given a certain index into the vectors
1125% rgh = reshape(rgh,1,[]);
1126% ratingPairs = [pls(:) rgh(:)];
1127% ratingPairs(isnan(ratingPairs(:, 1)), :) = []; % delete rows containing NaNs, as they cannot be processed by accumarray
1128% counts = accumarray(ratingPairs, 1); % matrix that counts how many times each rating pair occurs; e.g., at position (3,5) it shows how many exemplars(including unique reps) have been rated pls=3,rgh=5
1129% if max(max(counts)) > counts_max
1130% counts_max = max(max(counts));
1131% end
1132% end
1133
1134% main for
1135for i = 1:N % if we decided to pool across chains, this 'for' loops across seeds, and vice-versa
1136 pairs = []; % one such 2-columned matrix for every subplot, from which the the 2D histo is computed
1137 switch poolAcross
1138 case 'chains'
1139 subplotTitle = {['Seed ' num2str(i)]}; % this wouldnt necessarily have to be inside of an i_gen loop, but it's fine like that also, to reduce the number of switch statements
1140 for i_chain=1:N_chains % asta s-ar putea face si fara for loop, cu arrayfun, cf cz matlab
1141 for i_gen=1:N_gens
1142 for i_pair = 1:hums(i, i_chain).pairs.(DV).N(i_gen)
1143 pairs = [pairs ; ...
1144 hums(i, i_chain).pairs.(DV).element1(i_gen, i_pair), ...
1145 hums(i, i_chain).pairs.(DV).element2(i_gen, i_pair)];
1146 end
1147 end
1148 end
1149 case 'seeds'
1150 switch i
1151 case {1,2}
1152 chainType = 'mixed accuracy';
1153 case {3,4}
1154 chainType = 'high accuracy';
1155 end
1156 subplotTitle = {['Chain ' num2str(i) ' (' chainType ')']};
1157 for i_seed=1:N_seeds % asta s-ar putea face si fara for loop, cu arrayfun, cf cz matlab
1158 for i_gen=1:N_gens
1159 for i_pair = 1:hums(i_seed, i).pairs.(DV).N(i_gen)
1160 pairs = [pairs ; ...
1161 hums(i_seed, i).pairs.(DV).element1(i_gen, i_pair), ...
1162 hums(i_seed, i).pairs.(DV).element2(i_gen, i_pair)];
1163
1164 end
1165 end
1166 end
1167 end
1168 subplot(N_plots_x, N_plots_y, i);
1169
1170 % the 'histogram2' way of plotting the histograms
1171 h(i) = histogram2(pairs(:,1), pairs(:,2), -12:12, -12:12, 'DisplayStyle','tile', 'ShowEmptyBins','on'); % , 'XBinLimits', [-12 12], 'YBinLimits', [-12 12]
1172
1173% % the 'imagesc' way of plotting the histograms
1174% h(i) = imagesc(counts, [0 counts_max]) % sep21: in cazul meu i_genre ar fi pss i!
1175% gridN = 12; %refinement of map
1176% minvals = min(intervalPairs);
1177% maxvals = max(intervalPairs);
1178% rangevals = maxvals - minvals;
1179% xidx = 1 + round((intervalPairs(:,1) - minvals(1)) ./ rangevals(1) * (gridN-1));
1180% yidx = 1 + round((intervalPairs(:,2) - minvals(2)) ./ rangevals(2) * (gridN-1));
1181% density = accumarray([yidx, xidx], 1, [gridN,gridN]); % NB: y is rows, x is cols
1182% imagesc(density, 'xdata', [minvals(1), maxvals(1)], 'ydata', [minvals(2), maxvals(2)]);
1183
1184 % add heatmap legend, one for each subplot (either comment the row below or the further colormap ones below)
1185% h_colorbar(i) = colorbar;
1186% ylabel(h_colorbar(i), 'Bigram count', 'FontSize', fontSize_legendLabel)
1187
1188 cx(i,:) = caxis; % save away the intensity limits
1189 axis equal
1190 xlabel([DV(1:end-1) ' 1 (semitones' aboveMiddleC ')'], 'FontSize', fontSize_labels, 'FontName', fontName) % we do DV(1:end-1) to remove the plural, eg DV=notes becomes 'note' as the axis label
1191 set(gca, 'XTick', -12:4:12, 'FontSize', fontSize_values);
1192 set(gca, 'XTickLabel', [-12 -8 -4 0 4 8 12])
1193 ylabel([DV(1:end-1) ' 2 (semitones' aboveMiddleC ')'], 'FontSize', fontSize_labels, 'FontName', fontName)
1194 set(gca, 'YTick', -12:4:12, 'FontSize', fontSize_values); % inverseaza asta cumva, ca sa fie the same way up as in Zivic et al 2013; daca faci histograma cu imagesc, ala makes y axis increase top-bottom by default
1195 set(gca, 'YTickLabel', [-12 -8 -4 0 4 8 12])
1196% colormap gray
1197 title(subplotTitle, 'FontSize', 12, 'Units', 'normalized', 'Position', [0.55, 1.02], 'HorizontalAlignment', 'center', 'Color', [.5 .5 .5], 'FontWeight','bold')
1198end
1199
1200% normalise cell colour in each 2D histogram (subplot), that matches the colours shown in the heatmap (colormap)
1201% solutie de la https://www.mathworks.com/matlabcentral/answers/1449169-many-2d-histograms-single-colorbar
1202for i = 1:N
1203 subplot(N_plots_x, N_plots_y, i)
1204 caxis([min(cx(:,1)),max(cx(:,2))]) % From the smallest low to the largest high
1205end
1206
1207% add heatmap legend, namely a single one for all subplots, its range equal to the the widest range available among the various subplots
1208h_colorbar = colorbar;
1209h_colorbar.Position = [ 0.9224 0.1093 0.0133 0.8159 ];
1210ylabel(h_colorbar, 'Bigram count', 'FontSize', fontSize_legendLabel)
1211% colorbar('Ticks',[0 max(max(counts))], 'TickLabels',{'0',[num2str(max(max(counts)))]}) % this leads to different max values for the colorbar of each subplot
1212% colorbar('Ticks',[0 counts_max], 'TickLabels',{'0',[num2str(counts_max)]}); % this leads to a single range for all subplots (easier to compare visually)
1213
1214% save figure as raster image and as FIG
1215print('-dpng', ['d:\_matlabPlots\' plotName '.png' ] );
1216disp('gata...........................................................................................................................................................................')
1217
1218%% P1.1 PLOT various DVs that characterise each generation (or each inter-generational transmission) down the chain. By default opens 4 figures (one per chain), each with its 12 seeds
1219close all
1220clc
1221
1222% which measure do we want to pool across seeds/chains and plot? change both DV and DV_type!
1223% make sure that each of these DVs is only indexed by generation, and not also by note/interval number, as below we only use one index, namely for generation
1224DV = 'N_inflections'; % can be: 'meldevFromPrev' (intergen); 'distFromPrev' (intergen); 'diatonicity' or 'chromaticity' (gen); 'entropy_interval' (gen); 'N_inflections' (gen)
1225DV_type = 'gen'; % DVs that contain FromPrev refer to inter-generational transmission/transition ('intergen'), otherwise it's a feature of that generation itself ('gen')
1226
1227switch DV_type
1228 case 'gen'
1229 xLabel = 'Generation';
1230 xTicks = {1:8};
1231 case 'intergen'
1232 xLabel = 'Inter-generational transmission';
1233 xTicks = {'Seed/1' '1/2' '2/3' '3/4' '4/5' '5/6' '6/7' '7/8'};
1234end
1235
1236for i_chain = 1:N_chains
1237 % figure params that are DV-specific
1238 switch i_chain
1239 case {1,2}
1240 chainType = 'mxAcc';
1241 case {3,4}
1242 chainType = 'hiAcc';
1243 end
1244 switch DV
1245 case 'distFromPrev'
1246 plotName = {'Edit distances btw consecutive hums'' contours'...
1247 ['Expt. v' vers gender ', chain #' num2str(i_chain) ' (' chainType ')'] };
1248 yLabel = 'Edit actions';
1249 case 'meldevFromPrev'
1250 plotName = {'Melody deviation from invariance; low values = accurate transposition (high rigidity)'...
1251 ['Expt. v' vers gender ', chain #' num2str(i_chain) ' (' chainType ')'] };
1252 yLabel = 'Deviation';
1253 case {'diatonicity', 'chromaticity'}
1254 plotName = {[DV ', Expt. v' vers gender ', chain #' num2str(i_chain) ' (' chainType ')'] };
1255 yLabel = DV;
1256 case {'entropy_interval'}
1257 plotName = {[DV ', Expt. v' vers gender ', chain #' num2str(i_chain) ' (' chainType ')'] };
1258 yLabel = 'Entropy (bits)';
1259 case 'N_inflections'
1260 plotName = {[DV ', Expt. v' vers gender ', chain #' num2str(i_chain) ' (' chainType ')'] };
1261 yLabel = DV;
1262 end
1263 plotName_fig = cell2mat(plotName); % figure() doesnt accept cell names
1264 figure('Name', plotName_fig, 'NumberTitle','off', 'Color',[1 1 1], 'Position', [0 44 553 954])
1265
1266 for i_seed = 1:N_seeds
1267 subplot(N_seeds, 1, i_seed)
1268 varToPlot = hums(i_seed, i_chain).(DV)(:); % all these DVs are only indexed by generation, and not also by note/interval number, so it's ok that we only use one index, which gets us the values from all generations. We use var name constructed dyanmically - https://www.mathworks.com/help/matlab/matlab_prog/generate-field-names-from-variables.html
1269 plot(squeeze(varToPlot), 'o-')
1270 axis tight
1271 title({'Seed' num2str(i_seed)}, 'FontSize', 8, 'Units', 'normalized', 'Position', [-0.12 .25], 'HorizontalAlignment', 'center', 'Color', [.5 .5 .5], 'FontWeight','bold')
1272 box off
1273
1274 % apply labels only to the x axis of the final subplot line
1275 if i_seed == N_seeds
1276 xlabel(xLabel)
1277 ylabel(yLabel)
1278 set(gca, 'XTickLabel', xTicks);
1279 else
1280 set(gca, 'XTickLabel', {'' '' '' '' ''});
1281 end
1282 end
1283
1284 h_suptitle = suptitle(plotName, 0.95); % as the plot's title (suptitle), you may want to use a shortened version of the plotName; this command, if called earlier, ruins the title's horizontal alignment
1285 set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
1286
1287 print('-dpng', ['d:\_matlabPlots\' plotName_fig '.png' ] ); % save figure
1288end
1289disp('gata...........................................................................................................................................................................')
1290
1291%% P1.2 PLOT the above pooling across seeds OR chains
1292clc
1293close all
1294
1295clear n y err
1296
1297% which measure do we want to pool across seeds AND chains and plot? (see P1.1 for explanations)
1298DV = 'meldevFromPrev';
1299DV_type = 'intergen';
1300
1301poolAcross = 'seeds'; % 'chains', 'seeds' or 'chains and seeds'
1302
1303switch DV_type
1304 case 'gen'
1305 xLabel = 'Generation';
1306 xTicks = {1:8};
1307 case 'intergen'
1308 xLabel = 'Inter-generational transmission';
1309 xTicks = {'Seed/1' '1/2' '2/3' '3/4' '4/5' '5/6' '6/7' '7/8'};
1310end
1311
1312switch poolAcross % if we decided to pool across chains, the 'for' below loops across chains, and vice-versa
1313 case 'chains'
1314 N = N_seeds;
1315 case 'seeds'
1316 N = N_chains;
1317end
1318
1319switch DV
1320 case 'distFromPrev'
1321 plotName = {'Edit distances btw consecutive hums'' contours'...
1322 ['Expt. v' vers gender '; errorbars show variation across ' poolAcross] };
1323 yLabel = 'Edit actions';
1324 case 'meldevFromPrev'
1325 plotName = {'Melody deviation from invariance; low values = accurate transposition (high rigidity)'...
1326 ['Expt. v' vers gender '; errorbars show variation across ' poolAcross] };
1327 yLabel = 'Deviation';
1328end
1329
1330plotName_fig = cell2mat(plotName); % figure() doesnt accept cell names
1331figure('Name', plotName_fig, 'NumberTitle','off', 'Color',[1 1 1], 'Position', [0 44 766 740])
1332h_suptitle = suptitle(plotName, 0.95); % this command, if called earlier, ruins the title's horizontal alignment
1333set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
1334
1335for i = 1:N % if we decided to pool across chains, this 'for' loops across seeds, and vice-versa
1336 for i_gen = 1:N_gens
1337 values = []; % one such vector for every gen
1338 switch poolAcross
1339 case 'chains'
1340 subplotTitle = {['Seed ' num2str(i)]}; % this wouldnt necessarily have to be inside of an i_gen loop, but it's fine like that also, to reduce the number of switch statements
1341 for i_chain=1:N_chains
1342 values = [values hums(i, i_chain).(DV)(i_gen)]; % all these DVs are only indexed by generation, and not also by note/interval number, so it's ok that we only use one index, which gets us the values from all generations
1343 end
1344 case 'seeds'
1345 switch i
1346 case {1,2}
1347 chainType = 'mxAcc';
1348 case {3,4}
1349 chainType = 'hiAcc';
1350 end
1351 subplotTitle = {['Chain ' num2str(i)] ['(' chainType ')']};
1352 for i_seed=1:N_seeds
1353 values = [values hums(i_seed, i).(DV)(i_gen)]; % asta s-ar putea face si fara for loop, cu arrayfun, cf cz matlab
1354 end
1355 end
1356
1357 n(i_gen) = length(values);
1358 y(i_gen) = mean(values);
1359 err(i_gen) = std(values)/sqrt(length(values));
1360 end
1361 subplot(N, 1, i); errorbar(y,err)
1362 axis tight; box off
1363 xlim([0.75 N_gens])
1364 title(subplotTitle, 'FontSize', 8, 'Units', 'normalized', 'Position', [-0.1 .45], 'HorizontalAlignment', 'center', 'Color', [.5 .5 .5], 'FontWeight','bold')
1365
1366 % apply labels only to the x axis of the final subplot line
1367 if i == N
1368 xlabel(xLabel)
1369 ylabel(yLabel)
1370 set(gca, 'XTickLabel', xTicks);
1371 else
1372 set(gca, 'XTick', []);
1373 end
1374
1375 % add sample size above each data point, using the variable defined earlier, n(i_gen) = length(values); this can be used as a sanity-check but otherwise not really needed since this will always be n=4 when pooling across chains and n=12 when pooling across across seeds
1376 for i_gen=1:N_gens
1377 txt = ['N=' num2str(n(i_gen))];
1378 text(i_gen, y(i_gen) ,txt, 'FontSize', 6, 'fontweight','bold')
1379 end
1380end
1381print('-dpng', ['d:\_matlabPlots\' plotName_fig '.png' ] ); % save figure
1382disp('gata...........................................................................................................................................................................')
1383
1384%% P1.3 PLOT the above pooling across seeds AND chains. Doesn't require having previously run P1.2 or P1.1, only the base cell, 1.1
1385clc
1386close all
1387
1388clear n y err
1389
1390% which measure do we want to pool across seeds AND chains and plot? (see P1.1 for explanations)
1391DV = 'N_inflections';
1392DV_type = 'gen';
1393
1394switch DV_type
1395 case 'gen'
1396 xLabel = 'Generation';
1397 xTicks = {1:8};
1398 case 'intergen'
1399 xLabel = 'Inter-generational transmission';
1400 xTicks = {'Seed/1' '1/2' '2/3' '3/4' '4/5' '5/6' '6/7' '7/8'};
1401end
1402
1403switch DV
1404 case 'distFromPrev'
1405 plotName = {'Edit distances btw consecutive hums'' contours'...
1406 ['Expt. v' vers gender '; errorbars show variation across chains and seeds'] };
1407 yLabel = 'Edit actions';
1408 case 'meldevFromPrev'
1409 plotName = {'Melody deviation from invariance; low values = accurate transposition (high rigidity)'...
1410 ['Expt. v' vers gender '; errorbars show variation across chains and seeds'] };
1411 yLabel = 'Deviation';
1412 case {'diatonicity', 'chromaticity'}
1413 plotName = {[DV ', Expt. v' vers gender '; errorbars show variation across chains and seeds'] };
1414 yLabel = DV;
1415 case {'entropy_interval'}
1416 plotName = {['Entropy for intervals, Expt. v' vers gender '; errorbars show variation across chains and seeds'] };
1417 yLabel = 'Entropy (bits)';
1418 case 'N_inflections'
1419 plotName = {['Number of melodic inflections, Expt. v' vers gender '; errorbars show variation across chains and seeds'] };
1420 yLabel = 'Number of melodic inflections';
1421end
1422
1423% obtain the means & error bars for each gen, pooling across seeds and chains
1424for i_gen = 1:N_gens % out-most for since we want to get a grand average for each gen/transition
1425 values = [];
1426 for i_chain=1:N_chains
1427 for i_seed = 1:N_seeds
1428 values = [values hums(i_seed, i_chain).(DV)(i_gen)]; % all these DVs are only indexed by generation, and not also by note/interval number, so it's ok that we only use one index, which gets us the values from all generations
1429 end
1430 end
1431 n(i_gen) = length(values); % one such scalar value for every gen
1432 y(i_gen) = mean(values);
1433 err(i_gen) = std(values)/sqrt(length(values));
1434end
1435
1436% now plot those
1437plotName_fig = cell2mat(plotName); % figure() doesnt accept cell names
1438figure('Name', plotName_fig, 'NumberTitle','off', 'Color',[1 1 1], 'Position', [0 44 704 460])
1439h_suptitle = suptitle(plotName, 0.95); % this command, if called earlier, ruins the title's horizontal alignment
1440set(h_suptitle, 'FontSize',fontSize_suptitle, 'FontWeight','normal')
1441errorbar(y,err)
1442axis tight; box off
1443xlim([0.75 N_gens])
1444xlabel(xLabel)
1445ylabel(yLabel)
1446set(gca, 'XTickLabel', xTicks);
1447
1448% add sample size above each data point, using the variable defined earlier, n(i_gen) = length(values); this can be used as a sanity-check but otherwise not really needed since this will always be n=4 when pooling across chains and n=12 for seeds
1449for i_gen=1:N_gens
1450 txt = ['N=' num2str(n(i_gen))];
1451 text(i_gen, y(i_gen) ,txt, 'FontSize', 6, 'fontweight','bold')
1452end
1453
1454print('-dpng', ['d:\_matlabPlots\' plotName_fig '.png' ] ); % save figure
1455disp('gata...........................................................................................................................................................................')
1456
1457%% plot through-chain evolution of tune length (recording length in seconds)
1458clc
1459clear recordingDuration
1460
1461% plot
1462for i_seed = testSeed
1463 for i_chain = testChain
1464 clear y
1465 plotName = ['seed #' num2str(i_seed) ', chain #' num2str(i_chain) ' (' chainType '), in expt v' vers ' ' gender];
1466 figure('Name', plotName, 'NumberTitle','off', 'Color',[1 1 1])
1467 y(1:N_gens) = duration_hum(i_seed, i_chain, :);
1468 y = [ duration_seed y]; % first element of vector to be plotted is seed duration, thereafter hum duration at each gen
1469 plot(1:9, [NaN y(2:end)], '-o')
1470 hold on
1471 plot(1, y(1), 'sr', 'MarkerSize', 12)
1472 plot(1:2, y(1:2), '--r')
1473 xlabel('Generation #')
1474 set(gca, 'XTickLabel', {'Seed' '1' '2' '3' '4' '5' '6' '7' '8'})
1475 ylabel('recording duration (s)')
1476 end
1477end
1478disp('gata')
1479
1480%% plot through-chain evolution of tune length (number of "notes")
1481
1482