-
Notifications
You must be signed in to change notification settings - Fork 329
/
sparsity_9b_music.py
112 lines (101 loc) · 2.66 KB
/
sparsity_9b_music.py
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
def exo1():
"""
Compute and display $d_y$
VD
_y
isp
"""
[U, S, V] = svd(MusicHankel(y), 0); S = diag(S)
Ubot = U(: , N + 1: end)
d = Ubot'*exp(-2i*pi*(0: L-1)'*z(: )')
d = sum(abs(d).^2) / L
clf; hold on
plot(z, d, 'b')
stem(x0, 1 + x0*0, 'r.')
axis([0 1 0 1]); box on
def exo2():
"""
Display the roots of $P_y$ that are inside the unit disk
oefficients
oots
eep those inside
isplay
"""
B = []
for j in 1: size(Ubot, 2):
u = Ubot(: , j)
v = flipud(conj(u))
B(: , j) = conv(u, v)
C = sum(B, 2)
R = roots(C(end: -1: 1))
R = R(abs(R) <= 1)
clf; hold on
plot(exp(2i*pi*z), 'k')
plot(R, 'b.', 'MarkerSize', ms)
axis equal; box on
axis([-1 1 -1 1]*1.1)
axis off
def exo3():
"""
Keep only the best $N$ ones, i.e. those that are the closest from the
unit circle.
We denote those as $\tilde x \in \mathbb{T}^N$.
Recover an approximation $\tilde a$ of the amplitudes $a_0$
By solving in least squares (using backslash operator) the system
$ \Phi_{\tilde x} \tilde a = y. $
Display the recovered measure $ m_{\tilde a, \tilde x} $.
osition
eep only the best N ones.
ompute amplitude by solving a least square.
isplay the recovered measure.
"""
[x1, I] = sort(mod(angle(R), 2*pi)/ (2*pi))
R = R(I)
[~, I] = sort(abs(abs(R)-1))
x1 = x1(I(1: N))
R = R(I(1: N))
a1 = real(Phi(x1)\y)
clf; hold on
plot(z, f)
stem(x0, a0, 'r.')
stem(x1, a1, 'k--')
def exo4():
"""
Display the evolution of roots as the noise level $\sigma$ increases.
"""
slist = linspace(1e-9, 1, 1000)
RL = []
for is in 1: length(slist):
% observations
sigma = slist(is)
y = y0 + sigma*norm(y0)*w
% SVD
[U, S, V] = svd(MusicHankel(y), 0); S = diag(S)
Ubot = U(: , N + 1: end)
% Coefficients
B = []
for j in 1: size(Ubot, 2):
u = Ubot(: , j)
v = flipud(conj(u))
B(: , j) = conv(u, v)
C = sum(B, 2)
% Roots
R = roots(C(end: -1: 1))
% keep those inside
R = R(abs(R) <= 1)
% position
[x1, I] = sort(mod(angle(R), 2*pi)/ (2*pi))
R = R(I)
% Keep only the best N ones.
if 0
[~, I] = sort(abs(abs(R)-1))
x1 = x1(I(1: N))
RL(: , end + 1) = R
clf; hold on
plot(exp(2i*pi*z), 'k')
cm = jet(length(slist))
for is in 1: length(slist):
plot(RL(: , is), '.', 'MarkerSize', ms, 'color', cm(is, : ))
axis equal; box on
axis([-1 1 -1 1]*1.1)
axis off