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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
// Date of creation: 19 Dec, 2015
function [outputData,msg] = musicBase(inputData)
// Implements the core of the MUSIC algorithm
// Used by pmusic and rootmusic algorithm
// TODO: complete docs
msg = "";
outputData = struct();
[eigenvects,eigenvals, msg] = computeEig(inputData.x, ...
inputData.isCorrFlag, inputData.windowLength, ...
inputData.noverlap, inputData.windowVector, ...
inputData.isWindowSpecified);
// disp("eigenvects in musicBase");
// disp(eigenvects);
// disp("eigenvals in musicBase");
// disp(eigenvals);
if length(msg)~=0 then
return
end
pEffective = determineSignalSpace(inputData.p, eigenvals);
if length(msg)~=0 then
return
end
// Separating the eigenvects into signal and noise subspace
signalEigenvects = eigenvects(:,1:pEffective);
noiseEigenvects = eigenvects(:,pEffective+1:$);
outputData.signalEigenvects = signalEigenvects;
outputData.noiseEigenvects = noiseEigenvects;
outputData.eigenvals = eigenvals;
outputData.pEffective = pEffective;
endfunction
function [eigenvects,eigenvals,msg] = computeEig(x,isCorrFlag, windowLength, noverlap, window,isWindowSpecified)
// Computes the eigenvalues for the correlation matrix
// If x is a correlation matrix, which is specified using the isCorrFlag,
// spec() is used for eigenvalue decomposition.
// Otherwise, SVD is used for a proper restructure of x
// (i.e windowed version)
eigenvects = 0;
eigenvals = 0;
msg = "";
// determine if input is a matrix
xIsMatrix = ~or(size(x)==1);
if xIsMatrix & isCorrFlag then
// TODO: check the order of eigenvalues
[eigenvects,d] = spec((R+R')/2); // ensuring hermitian property
eigenvals = diag(d);
// sorting in decreasing order
[eigenvals,order] = gsort(eigenvals);
// TODO: nonnegative eigenvals check
// rearragning in decreasing order of eigenvalues
eigenvects = eigenvects(:,order);
else
if xIsMatrix then
// TODO: check for dimenion constraints
else
// x is vector
Lx = length(x);
if Lx<=windowLength then
msg = "Incorrrect value for window length; must be smaller than the signal length";
return
end
if ~isWindowSpecified then
// disp("window not specified");
// TODO: understand
[x,msg] = createBufferMatrix(x,Lx-windowLength+1,Lx-windowLength);
if length(msg)~=0 then
return
end
// reversing the column order and scaling
x = x(:,$:-1:1)./sqrt(Lx-windowLength+1);
else
// disp("window specified");
[x,msg] = createBufferMatrix(x, windowLength, noverlap);
if length(msg)~=0 then
return
end
// scaling so as to get the correct value of R
x = x'./sqrt(Lx-windowLength);
end
end
// **applying window to each row of the data matrix**
// replicating window along the rows
if ~isempty(window) then
window = repmat(window(:)',size(x,1),1);
x = x.*window;
end
// computing eignevals and eigenvectors of R using SVD of x
// disp("X = (before SVD)");
// disp(x);
[temp,S,eigenvects] = svd(x,0);
// squaering the eigenvalues
eigenvals = diag(S).^2;
// disp("eigenvals in computeEig");
// disp(eigenvals);
end
endfunction
function [xMat,msg] = createBufferMatrix(x,windowLength,noverlap)
// creates a matrix where each row represents a section which has to be
// windowed
//
// Input Arguments
// x - input signal as a column vector
// windowLength
// noverlap
//
// will perform task similar to that performed by MATLAB's
// buffer(x,windowLength,noverlap) with nodelay option
msg="";
xMat = [];
// check input to be a vector
xIsVector = or(size(x)==1) & ndims(x)==2;
if ~xIsVector then
msg = "createBufferMatrix: x should be a vector";
return
end
if size(x,2)~=1 then
// convert to column vector
x = x';
end
L = length(x);
temp = windowLength - noverlap;
numOfSections = ceil((L-noverlap)/temp);
// performing zero padding of x
zeroPadLength = numOfSections*temp + noverlap - L;
zeroPad = zeros(zeroPadLength,1);
x = [x;zeroPad];
xMat = zeros(windowLength, numOfSections);
// disp(size(xMat));
// disp(size(x));
for i=1:numOfSections
xMat(1:temp,i) = x(1+(i-1)*temp:i*temp,1);
xMat(temp+1:windowLength,i) = x(1+i*temp:i*temp+noverlap,1);
end
endfunction
function pEffective = determineSignalSpace(p, eigenvals)
// Determines the effective dimension of the signal subspace
// Inputs:
// p: p(1) - signal subspace dimension
// p(2) (optional) - desired threshold
// eigenvals: vector conatining eigenvalues of the correlation
// matrix in descreasing order
// Output:
// pEffective - the effective dimension of the signal subspace. If
// a threshold is given as p(2), the signal subspace will
// be equal to the number of eigenvalues greater than the
// threshold times the smallest eigenvalue. Also,
// pEffective<=p(1). So, minimum of the two values is
// considered. If the threshold criteria results in an
// empty signal subspace, we take pEffective = p(1).
//
if length(p)==2 then
threshold = p(2)*eigenvals($);
signalIndices = find(eigenvals>threshold);
if ~isempty(signalIndices) then
pEffective = min(p(1), length(signalIndices));
else
// dont change p
pEffective = p(1);
end
else
pEffective = p;
end
endfunction
|
