-
Notifications
You must be signed in to change notification settings - Fork 0
/
PfNmf.m
127 lines (105 loc) · 3.22 KB
/
PfNmf.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
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
%% Partially Fixed Nonnegative Matrix Factorization
% [WD, HD, WH, HH, err] = PfNmf(X, WD, HD, WH, HH, rh, sparsity)
% input:
% X = float, numFreqX*numFrames matrix, input magnitude spectrogram
% WD = float, numFreqD*rd matrix, drum dictionary
% HD = float, rd*numFrames matrix, drum activation matrix
% WH = float, numFreqH*rh matrix, harmonic dictionary
% HH = float, rh*numFrames matrix, harmonic activation matrix
% rh = int, rank of harmonic matrix
% sparsity = float, sparsity coefficient
% output:
% WD = float, numFreqD*rd matrix, updated drum dictionary
% HD = float, rd*numFrames matrix, updated drum activation matrix
% WH = float, numFreqH*rh matrix, updated harmonic dictionary
% HH = float, rh*numFrames matrix, updated harmonic activation matrix
% err = error vector (numIter * 1)
% usage:
% To randomly initialized different matrix, please give [] as input.
% For example, [WD,HD,WH,HH,err] = PfNmf(X, WD, [], [], [], 0, 0)
% is the basic NMF approach given only the drum template.
% [WD,HD,WH,HH,err] = PfNmf(X, WD, [], [], [], 50, 0) is the
% partially fixed NMF with 50 random intialized entries
%
% CW @ GTCMT 2015
function [WD, HD, WH, HH, err] = PfNmf(X, WD, HD, WH, HH, rh, sparsity)
X = X + realmin; %make sure there's no zero frame
[numFreqX, numFrames] = size(X);
[numFreqD, rd] = size(WD);
%initilization
WD_update = 0;
HD_update = 0;
WH_update = 0;
HH_update = 0;
if WH
[numFreqH, rh] = size(WH);
else
WH = rand(numFreqD, rh);
[numFreqH, ~] = size(WH);
WH_update = 1;
end
if (numFreqD ~= numFreqX)
error('Dimensionality of the WD does not match X');
elseif (numFreqH ~= numFreqX)
error('Dimensionality of the WH does not match X');
end
if HD
WD_update = 1;
else
HD = rand(rd, numFrames);
HD_update = 1;
end
if HH
else
HH = rand(rh, numFrames);
HH_update = 1;
end
alpha = (rh + rd)/rd;
beta = rh/(rh + rd);
%normalize W / H matrix
for i = 1:rd
WD(:,i) = WD(:,i)./(norm(WD(:,i),1));
end
for i = 1:rh
WH(:,i) = WH(:,i)./(norm(WH(:,i),1));
end
count = 0;
rep = ones(numFreqX, numFrames);
%start iteration
while (count < 300)
approx = alpha*WD*HD + beta*WH*HH;
%update
if WD_update
WD = WD .* ((X./approx)*(alpha * HD)')./(rep*(alpha * HD)');
else
end
if HD_update
HD = HD .* ((alpha * WD)'* (X./approx))./((alpha * WD)'*rep + sparsity);
else
end
if WH_update
WH = WH .* ((X./approx)*(beta * HH)')./(rep*(beta * HH)');
else
end
if HH_update
HH = HH .* ((beta * WH)'* (X./approx))./((beta * WH)'*rep);
else
end
%normalize W matrix
for i = 1:(rh)
WH(:,i) = WH(:,i)./(norm(WH(:,i),1));
end
for i = 1:(rd)
WD(:,i) = WD(:,i)./(norm(WD(:,i),1));
end
%normalize H matrix
%calculate variation between iterations
count = count + 1;
err(count) = KlDivergence(X, (alpha * WD*HD + beta * WH*HH)) + sparsity * norm(HD, 1);
if (count >=2)
if (abs(err(count) - err(count -1 )) / (err(1) - err(count) + realmin)) < 0.001
break;
end
end
end
end