-
Notifications
You must be signed in to change notification settings - Fork 1
/
draw6.m
61 lines (50 loc) · 1.87 KB
/
draw6.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
% For real-world audio example with multiple bursts,
% find all samples T exceeding the threshold
% for any two samples in T, if abs(t2-t1)<b, we determine all samples in
% between belong to noise
% corresponding to 2.4.2 in the paper
% This script plots the corrupted signal, the detection function, and the
% detected results corresponding to various values of params K and b
clear; close
x = audioread('source_Muss_l.wav'); x = x(:,1);
% 2000 samples of audio file sampled at 44100
x = x(8001:10000);
subplot(2,1,1); plot(x); xlabel('sample number'); ylabel('amplitude')
title('Corrupted signal')
% assume maximum length of a burst
Nmax = 50;
% estimate AR parameters
p = 3*Nmax + 2;
[A, e] = aryule(x, p); % 1,a1,a2...
% compute detection function d
% d(t) = 1*x(t)+a1*x(t-1)+a2*x(t-2)...+ap*x(t-p)
d = filter(A, 1, x);
d(1:p) = d(1:p)*0; % d is only defined for t>p
d = abs(d);
% parameters to tune
b1 = 1; b2 = 30;
K1 = 1; K2 = 2; K3 = 3;
thre1 = K1*sqrt(e); % e is estimated variance of excitation
thre2 = K2*sqrt(e);
thre3 = K3*sqrt(e);
subplot(2,1,2); plot(d); xlabel('sample number'); ylabel('amplitude')
title('Detection signal with thresholds')
r = yline(thre1, 'r--');
g = yline(thre2, 'g--');
b = yline(thre3, 'b--');
r.LineWidth = 1.5;
g.LineWidth = 1.5;
b.LineWidth = 1.5;
figure
pos1 = thresholding(d, thre1, b1);
pos2 = thresholding(d, thre2, b1);
pos3 = thresholding(d, thre3, b1);
subplot(3,2,1); plot(pos1); title('K=1 b=1'); ylim([-0.1, 1.1])
subplot(3,2,3); plot(pos2); title('K=2 b=1'); ylim([-0.1, 1.1])
subplot(3,2,5); plot(pos3); title('K=3 b=1'); ylim([-0.1, 1.1])
pos4 = thresholding(d, thre1, b2);
pos5 = thresholding(d, thre2, b2);
pos6 = thresholding(d, thre3, b2);
subplot(3,2,2); plot(pos4); title('K=1 b=30'); ylim([-0.1, 1.1])
subplot(3,2,4); plot(pos5); title('K=2 b=30'); ylim([-0.1, 1.1])
subplot(3,2,6); plot(pos6); title('K=3 b=30'); ylim([-0.1, 1.1])