-
Notifications
You must be signed in to change notification settings - Fork 1
/
PulseEnvelope.m
88 lines (70 loc) · 2.64 KB
/
PulseEnvelope.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
function [Start,Center,Stop] = PulseEnvelope(Data, Pulses, MaxExtend, step, NumStd)
%Pulses = most inclusive Pulses structure that includes start and stop
%times
%MaxExtend = maximum extension to start and stop of original pulse to try
%(efault = 150)
%step = step size to try from original pulse to MaxExtend (default = 10)
%NumStd= # standard deviations of noise to use as cutoff (default = 2)
addpath(genpath('./chronux'))
if nargin < 3
MaxExtend = 150;
end
if nargin < 4
step = 10;
end
if nargin < 5
NumStd = 2;
end
%fprintf(1,' Finding Noise Level\n');
if numel(Data.d) > 1e6
shortdata = Data.d(1:1e6);
else
shortdata = Data.d;
end
[ssf] = sinesongfinder(shortdata,Data.fs,12,20,.1,.01,.05,[0 1000]); %returns ssf, which is structure containing the following fields
noise = findnoise(ssf,3,80,1000);
StDev = noise.sigma * NumStd;
pulseData = Pulses;
%check here for pulses that start or stop too close to end (within ±
%Extend) -- delete these pulses
KeepIdx = find(pulseData.w0 > MaxExtend & pulseData.w1 < (numel(Data.d) - MaxExtend)); %idx of pulses that are < length of song - MaxExtend
pulseData.w0 = pulseData.w0(KeepIdx);
pulseData.wc = pulseData.wc(KeepIdx);
pulseData.w1 = pulseData.w1(KeepIdx);
AllStarts = pulseData.w0;
NumClips = numel(AllStarts);
for j = 0:step:MaxExtend
Extend = j;
clips = GetClips(pulseData.w0-Extend,pulseData.w1+Extend,Data.d);
k = cell2mat(clips);
SmoothedClips = smooth(k);
SmoothedClips= reshape(SmoothedClips,size(k,1),size(k,2));
Env = smooth(abs(hilbert(SmoothedClips)));
Env = reshape(Env,size(k,1),size(k,2));
[~,MaxIdx]=max(Env);
%[~,center] = max(abs(SmoothedClips));
Start = NaN(NumClips,1);
Stop = Start;
Center = Start;
for i = 1:NumClips;
start = find(Env(1:MaxIdx,i) < (StDev),1,'last');
stop = find(Env(MaxIdx:end,i) < (StDev),1,'first');
% if strcmp(StartOrCenter,'Start')
if ~isempty(stop) && ~isempty(start)
Start(i) = AllStarts(i) - Extend + start; %start at estimated pulse start
Stop(i) = AllStarts(i) - Extend + stop + MaxIdx(i);
% end
% elseif strcmp(StartOrCenter,'Center')
% if ~isempty(stop)
Center(i) = pulseData.wc(i); %AllStarts(i) -Extend + center(i); %start at pulse center
% Stop(i) = AllStarts(i) - Extend + stop + MaxIdx(i);
% end
end
end
Start = Start(~isnan(Start));
Stop = Stop(~isnan(Stop));
Center = Center(~isnan(Center));
if numel(Start) > NumClips * 0.99
break
end
end