Meta-Analysis of Face-Validity Indicators in Automatically Detected Dog Sleep Spindles.

1. %Written to search for spindles in dog EEG, special thanks to Róbert
2. %Bódizs, Gilles van Luijtelaar, Kaue Costa, and Antoine Nonclercq.
4. awakeeeg = FzCz(hypnogram == 1); %define which part of the recording is awake eeg, using the corresponding hypnogram label '1' (awake)
5. sleepeeg = FzCz(hypnogram == 4); %do the same with the non-REM part of the signal
7. sleepeeg = sleepeeg - mean(sleepeeg);
9. raweeg = sleepeeg;
11. Nyquist = SR/2; %Using SR (= Sampling Rate, provided with the file) calculate the Nyquist frequency
13. window = SR*0.5; %Define the size of the signal segments that will be used for a sliding fft, as a magnitude of sampling rates (results in 0.5 seconds window)
15. normalized_cutoff = 16/Nyquist; %passband low-pass
17. normalized_cutoff2 = 35/Nyquist; %stopband low-pass
19. nc = 30/Nyquist;
21. nc2 = 35/Nyquist;
23. highpass_cutoff2 = 5/Nyquist; %passband high-pass
25. highpass_cutoff = 3/Nyquist; %stopband low-pass
27. [order,nfreq] = buttord(normalized_cutoff,normalized_cutoff2,0.5,30); %establish the best order and cut-off frequency for a low-pass filter that attenuates with 30 dB and allows for max 0.5 dB ripple
29. [order2,nfreq2] = buttord(highpass_cutoff2,highpass_cutoff,0.5,30); %establish the best order and cut-off frequency for a high-pass filter that attenuates with 30 dB and allows for max 0.5 dB ripple
31. [order3,nfreq3] = buttord(nc,nc2,0.5,30); %establish the best order and cut-off frequency for a low-pass filter that attenuates with 30 dB and allows for max 0.5 dB ripple
33. [b,a,k1] = butter(order,nfreq,'low'); %obtaining filter coefficients with the order and cut-off specified for the low-pass filter
35. [c,d,k2] = butter(order2,nfreq2,'high'); %obtaining filter coefficients with the order and cut-off specified for the high-pass filter
37. [e,f,k3] = butter(order3,nfreq3,'low'); %obtaining filter coefficients with the order and cut-off specified for the low-pass filter
39. sos1 = zp2sos(b,a,k1);
41. sos2 = zp2sos(c,d,k2);
43. sos3 = zp2sos(e,f,k3);
45. sleepeeg = sosfilt(sos1,sleepeeg); %transform the trace with the low-pass filter
47. sleepeeg = sosfilt(sos2,sleepeeg); %transform the trace with the high-pass filter
49. raweeg = sosfilt(sos3,raweeg); %transform the trace with the low-pass filter
51. clear b a c d %clear the filter coefficients
53. [result, frequencymap, timepoints, power] = spectrogram(sleepeeg,window,round(SR*0.125),SR*10,SR); %obtain a matrix of amplitudes, the corresponding frequencies, timepoints and power for a sliding fft with 0.5 seconds window-size (using predefined window value, see above), window-overlapp of 0.125 seconds, and zero-padding for 0.1 Hz resolution (10*SR)
55. result = abs(result');
56. frequencymap = frequencymap';
57. timepoints = timepoints';
58. power = power'; %align all results to ease later indexing
61. [r,c] = size(result); %obtain size parameters of the results for subsequent loop-construction
64. for x = 1:r %loop repeats for the number of windows used in the fft
66. szelet = (frequencymap >= 0); %originally at this spot only frequencies above 2 Hz were selected for further analysis, to avoid the estimation of the max frequency power be influenced by the DC-offset, however I later learned high-pass filtering deals with it, so the value is set to zero now
67. szeletmap = frequencymap(frequencymap >= 0); %here the frequency-map is adjusted for the selection above, which is now simply above 0 (but was above 2 Hz)
68. selectresult = power(x,szelet); %the power-values of each window are collected, originally only for the selected frequency search radius (> 2 Hz)
69. [amp,pos] = max(selectresult); %amplitude and position on the frequency-map of the highest power is selected
70. maxfreq(x) = pos; %for each window the frequency with the highest power is collected in a vector
71. absmaxpow(x) = amp;
72. clear amp pos selectresult %clear variables for a new loop
74. end
78. forcrd = timepoints*SR; %reconstruct the sample positions of the time-points :D
81. for e = 1:length(timepoints); %go through all time-points
82. crd = forcrd(e); %find each corresponding sample
83. crd1 = round(crd - (0.25*SR)); %select a point that is half of the fft-window (in samples) below the time-point (which is itself the center of the window)
84. crd2 = round(crd + (0.25*SR)); %select a point that is half of the fft-window (in samples) above the time-point (which is itself the center of the window)
85. if crd1 <= 0 %if the lower boundary of the real window is shorter than expected (the center is not symmetrical)
86. crd1 = 1; %let the beginning of the whole trace define the beginning of the window
87. end
88. if crd2 > length(sleepeeg); %if the upper boundary of the window is longer than expected (signal ends less than 0.25 seconds after the center of the window)
89. crd2 = length(sleepeeg); %upper boundary of the window coincides with the end of the trace
90. end
92. ampest(e) = rms(sleepeeg(crd1:crd2)); %calculate the root mean square (rms) for the entire window as an amplitude estimate
93. clear crd crd1 crd2 %clear variables for a new loop
94. end
97. ampest = ampest'; %align amplitude-estimates for better indexing
99. ampest = zscore(ampest); %calculate standard-scores for the amplitude-estimates
101. maxfreq = maxfreq'; %align the fequencies (per window) with the max power for better indexing
103. range = (szeletmap(maxfreq) >= 9 & szeletmap(maxfreq) <= 16)'; %select those windows whose max frequency power is within the target range
105. selector = range & (ampest >= 1); %combine that with the amplitude estimate being more than a standard deviation above the average trace rms
107. happening_times = timepoints(selector == 1,:); %use the indices that correpsond to above criteria to select the windows that likely contain spindle activity
111. sampleTimes = timepoints*SR; %this and the next few lines are a repetition of the amplitude estimation, but this time only to arrive at the amplitudes of the above selected windows
115. for k = 1:length(range);
116. coord = sampleTimes(k);
117. coord1 = round(coord - (0.25*SR));
118. coord2 = round(coord + (0.25*SR));
119. if coord1 <= 0
120. coord1 = 1;
121. end
122. if coord2 > length(sleepeeg);
123. coord2 = length(sleepeeg);
124. end
126. valuesOFtimes(k) = rms(sleepeeg(coord1:coord2));
127. clear coord coord2
128. end
131. valuesOFtimes = valuesOFtimes';
133. valuesOFtimes = zscore(valuesOFtimes);
135. samp = valuesOFtimes(selector == 1); %this repetition of steps, and the selection of specifically amplitudes of spindle-containing suspected windows, ends here
137. [estimatesamp,pci_amp] = mle(samp); %estimatesamp contains estimates of the true mean(1) and standard deviation(2) of the spindle-candidate-amplitudes
139. medspin = estimatesamp(1); %label the mean estimate of amplitudes
140. varamp = estimatesamp(2); %label the std estimate of amplitudes
142. varamp = 2*varamp; %establish a 2 std range for amplitudes
144. uplimit = medspin + varamp; %establish upper-limit for acceptable amplitude values based on the 2 std range
146. lolimit = medspin - varamp; %establish lower-limit for acceptable amplitude values based on the 2 std range
148. [estimatesfreq,pci2] = mle(szeletmap(maxfreq(selector == 1))); %estimates are now obtained for the frequencies of the spindle-suspects
150. mfv = estimatesfreq(1); %label the estimate of frequency mean
151. mfs = estimatesfreq(2); %label the estimate of frequency std
153. mfs = 2*mfs; %the frequency 2 std range
155. minif = mfv - mfs; %lowerlimit frequency
156. maxf = mfv + mfs; %upperlimit frequency
158. range2 = (szeletmap(maxfreq) >= minif & szeletmap(maxfreq) <= maxf)'; %new version of above variable 'range' but now based on the upper and lower limits of frequency (normal-modeling)
160. selector2 = range2 & (valuesOFtimes > lolimit) & (valuesOFtimes < uplimit) & (valuesOFtimes >= 1); %second selection based on the new normal modelling criteria
162. happening_times = timepoints(selector2 == 1,:); %marks the correpsonding time-points of the last selection, to use them for plotting
164. h_times = happening_times; %renaming, possibly now redundant
166. for t = 1:length(h_times) %plotting loop
168. mark1 = h_times(t); %go through all time-points corresponding to detections as defined after the second round of selection
169. dvalue = mark1*SR; %find the corresponding sample in the original signal to each time-point
170. svalue = dvalue + (SR*2.5); %extend 1.5 seconds from the time-point of detection and use the correpsonding sample for the upper boundary of the plotting window
171. evalue = dvalue - (SR*2.5); %go 1.5 seconds down from the time-point of detection and use the correpsonding sample for the lower boundary of the plotting window
172. if evalue <= 0 %if the time-point of detection is less than 1.5 seconds from the start of the signal...
173. evalue = 1; %...then start plotting from the beginning of the signal
174. elseif svalue >= length(sleepeeg) %if the signal ends less than 1.5 seconds from the time-point of detection...
175. svalue = length(sleepeeg); %...then only plot as far as the signal reaches
176. end
177. xax = 0:1/SR:5; %a scale in seconds for the x-axis
178. xax = xax - 2.5; %let the detection be time-point zero on the scale
179. plotsignal = raweeg(evalue:svalue); %determine the part of the signal you want to plot
180. if length(xax) < length(plotsignal)||length(xax) > length(plotsignal) %if the portion you want to plot is longer or shorter than the 3 second time-window around detection...
181. plot(plotsignal) %...then just plot without scaling
182. else
183. plot(xax,plotsignal) %if scale and signal align always scale
184. ydown = min(plotsignal) - 10;
185. yup = max(plotsignal) + 10;
186. axis([-2.5 2.5 ydown yup]) %limit the y-axis to diplay the signal only between 20 and -20 µV
187. end
188. pause %pause after showing every detection-image
189. clear mark1 dvalue svalue evalue %eliminate variables that could interfer with the loop
191. end
194. for w = 1:length(h_times) %for calculating negative chirping index
195. maxcoordf = find(timepoints == h_times(w)); %establish which row of the time-point vector corresponds to the center of a spindle-detection
196. prevcoord = maxcoordf - 1; %find the position preceding the center of a spindle
197. postcoord = maxcoordf + 1; %find the position following the center of a spindle
198. if prevcoord > 1
199. spstartfreq = maxfreq(prevcoord); %find the maximum frequency preceding the center of a spindle
200. else
201. spstartfreq = maxfreq(maxcoordf);
202. end
203. if postcoord < length(timepoints);
204. spstopfreq = maxfreq(postcoord); %find the maximum frequency following the center of a spindle
205. else
206. spstopfreq = maxfreq(maxcoordf);
207. end
208. negchirp_value(w) = spstartfreq/spstopfreq; %calculate how much higher is the maximum frequency preceding the center of a spindle compared to the maximum frequency that follows
209. clear maxcoordf prevcoord postcoord spstartfreq spstopfreq
210. end
213. HD = diff(h_times); %calculate a vector of distances between detection times
214. UjHD = HD > 0.5; %count 1 separate spindle event for every two detections that are more than 0.5 seconds apart
215. spindle_number = sum(UjHD); %sum the previous ones to estimate the number of spindles
217. if (length(h_times)) - 1 == sum(HD > 0) %if the number of detections and events estimated as separate is equal...
218. spindle_number = spindle_number + 1 %...take the number of estimated separate events and add a one to adjust for the shorter diff-vector
219. elseif sum(HD > 0) == 1 %if there is only one 1 in UjHD, this still means 2 separate spindles
220. spindle_number = 2
221. elseif length(h_times) > 0 & sum(HD > 0) == 0 %if no event is marked as separate, but the detections are larger zero...
222. spindle_number = 1 %this is still one spindle
223. else
224. spindle_number = spindle_number %if non of above complications occur, spindle_number is a good enough estimate :)
225. end
229. average_frequency = mean(szeletmap(maxfreq(selector2 == 1))) %frequency is averaged across all detections of the second detection-round
231. SWSminutes = (length(sleepeeg)/SR)/60 %currently not activated, but can calculate the duration of the non-REM sleep in minutes, has been used to obtain data already
233. density_spindles = spindle_number/SWSminutes %the estimate of spindle number divided by the minutes of non-REM sleep
235. amplitude = mean(valuesOFtimes(selector2 == 1)) %the amplitude is averaged across all detections of the second detection-round
237. absmaxpow = absmaxpow';
239. avg_maxpow = mean(absmaxpow(selector2 == 1))
241. overlapdur = ((sum(selector2)) - spindle_number + 1)*0.3750;
243. noverlapdur = (spindle_number)*0.5;
245. total_duration = overlapdur + noverlapdur;
247. avg_duration = total_duration/spindle_number
249. avg_negative_chirp_magnitude = mean(negchirp_value)
251. transhyp = hypnogram < 2; %these last lines estimate how often the dog slipped out of sleep
252. awakenings = sum(diff(transhyp) == -1)
253. awaketime = (length(awakeeeg)/SR)/60

spindle_hound_rawplot_extras.m at commit fa87bc9, no license · at the source