Blind Source Separation

In [1]:
import Mads
In [2]:
import NMF
In [5]:
srand(2015)
nk = 3
s1 = (sin(0.05:0.05:5)+1)/2
s2 = (sin(0.3:0.3:30)+1)/2
s3 = rand(100);

Source matrix (assumed unknown)

In [6]:
S = [s1 s2 s3]
Out[6]:
100×3 Array{Float64,2}:
 0.52499      0.64776     0.518466 
 0.549917     0.782321    0.249446 
 0.574719     0.891663    0.244758 
 0.599335     0.96602     0.304535 
 0.623702     0.998747    0.0407869
 0.64776      0.986924    0.258012 
 0.671449     0.931605    0.0868265
 0.694709     0.837732    0.882364 
 0.717483     0.71369     0.850562 
 0.739713     0.57056     0.216851 
 0.761344     0.421127    0.333606 
 0.782321     0.27874     0.285782 
 0.802593     0.156117    0.889261 
 ⋮                                 
 0.0171135    0.999997    0.577248 
 0.0112349    0.978188    0.956294 
 0.00657807   0.913664    0.62366  
 0.0031545    0.812189    0.914383 
 0.000972781  0.682826    0.578899 
 3.83712e-5   0.537133    0.515968 
 0.000353606  0.388122    0.415338 
 0.0019177    0.249105    0.961519 
 0.00472673   0.1325      0.795982 
 0.00877369   0.0487227   0.73677  
 0.0140485    0.00525646  0.653462 
 0.0205379    0.00598419  0.607898
In [9]:
Mads.plotseries(S, title="Original sources", name="Source", combined=true)
Out[9]:
X -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 Source 1 Source 2 Source 3 Original sources -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 -1 0 1 2 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 Y

Mixing matrix (assumed unknown)

In [7]:
H = [[1,1,1] [0,2,1] [1,0,2] [1,2,0]]
Out[7]:
3×4 Array{Int64,2}:
 1  0  1  1
 1  2  0  2
 1  1  2  0

Data matix (known)

In [10]:
X = S * H
Out[10]:
100×4 Array{Float64,2}:
 1.69122   1.81399   1.56192   1.82051  
 1.58168   1.81409   1.04881   2.11456  
 1.71114   2.02808   1.06423   2.35805  
 1.86989   2.23657   1.2084    2.53137  
 1.66324   2.03828   0.705276  2.6212   
 1.8927    2.23186   1.16378   2.62161  
 1.68988   1.95004   0.845102  2.53466  
 2.4148    2.55783   2.45944   2.37017  
 2.28173   2.27794   2.41861   2.14486  
 1.52712   1.35797   1.17341   1.88083  
 1.51608   1.17586   1.42856   1.6036   
 1.34684   0.843262  1.35389   1.3398   
 1.84797   1.20149   2.58111   1.11483  
 ⋮                                      
 1.59436   2.57724   1.17161   2.01711  
 1.94572   2.91267   1.92382   1.96761  
 1.5439    2.45099   1.2539    1.83391  
 1.72973   2.53876   1.83192   1.62753  
 1.2627    1.94455   1.15877   1.36663  
 1.05314   1.59023   1.03198   1.0743   
 0.803814  1.19158   0.83103   0.776598 
 1.21254   1.45973   1.92496   0.500128 
 0.933209  1.06098   1.59669   0.269727 
 0.794267  0.834216  1.48231   0.106219 
 0.672767  0.663975  1.32097   0.0245614
 0.63442   0.619866  1.23633   0.0325062
In [13]:
Mads.plotseries(X, title="Mixed signals", name="Signal", combined=true)
Out[13]:
X -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 Signal 1 Signal 2 Signal 3 Signal 4 Mixed signals -4 -3 -2 -1 0 1 2 3 4 5 6 7 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 -3 0 3 6 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 Y

Blind signal reconstruction

In [15]:
Wipopt, Hipopt, pipopt = Mads.NMFipopt(X, nk, retries=1);

OF = 4.5864774801052004e-14

Reconstructed sources

In [16]:
Mads.plotseries(Wipopt, title="Reconstructed sources", name="Source", combined=true)
Out[16]:
X -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 Source 1 Source 2 Source 3 Reconstructed sources -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 -1 0 1 2 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 Y

Reproduced signals

In [17]:
Mads.plotseries(Wipopt * Hipopt, title="Reproduced signals", name="Signal", combined=true)
Out[17]:
X -150 -100 -50 0 50 100 150 200 250 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 -100 0 100 200 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 Signal 1 Signal 2 Signal 3 Signal 4 Reproduced signals -4 -3 -2 -1 0 1 2 3 4 5 6 7 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 -3 0 3 6 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 Y
In [ ]: