In [6]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import ifft, idct, dct

%matplotlib inline

# np.random.seed(0)
In [58]:
x = np.arange(0.0, 3.0, 0.01)

# compute the value (amplitude) of the sin wave at the for each sample
t_697  = np.sin(2*np.pi*697 * (x)) 
t_1209 = np.sin(2*np.pi*1209 * (x)) 
sum_t = t_697 + t_1209

sum_t.shape


plt.figure(1)
plt.subplot(311)
plt.plot(x,t_697)

plt.subplot(312)
plt.plot(x,t_1209)

plt.subplot(313)
plt.plot(x,sum_t)
Out[58]:
[<matplotlib.lines.Line2D at 0x7f60b5de88d0>]
In [3]:
from scipy.fftpack import ifft, idct, dct

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.dct.html
#                      N-1
# y[k] = [scaling]* 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
#                      n=0


f_i  = dct(sum_t, type=2, norm='ortho') #this does the same transform as neil's

plt.plot(x,f_i)
Out[3]:
[<matplotlib.lines.Line2D at 0x7f0b5fee1518>]
In [4]:
f_i_inverse  = idct(f_i, type=2, norm='ortho') #this does the same transform as neil's
diff_t = sum_t - f_i_inverse
plt.plot(x,diff_t)
Out[4]:
[<matplotlib.lines.Line2D at 0x7f0b5f52a2e8>]
In [16]:
np.random.seed(0)
M = 100

sum_t_sample = np.zeros(len(sum_t))
sum_t_sample[sum_t_sample==0]=np.nan

sum_t_sample_index = np.random.choice(range(len(sum_t)), M, replace = False)
sum_t_sample[sum_t_sample_index] = sum_t[sum_t_sample_index]

plt.plot(x,sum_t_sample,linestyle="",marker="o", color = 'blue')
plt.plot(x,sum_t, color = 'blue')

sum_t_sample.shape
Out[16]:
(300,)
In [17]:
sum_t_sample_new = sum_t_sample[np.isfinite(sum_t_sample)]
plt.plot(range(M), sum_t_sample_new, color = 'orange')
plt.plot(range(M), sum_t_sample_new, color = 'orange', marker="o")

f_i_dash  = dct(sum_t_sample_new, type=2, norm='ortho') #this does the same transform as neil's
plt.plot(range(M),sum_t_sample_new, color = 'orange')
plt.plot(range(M),idct(f_i_dash, type=2, norm='ortho'), marker="o", linestyle="", color = 'blue')
Out[17]:
[<matplotlib.lines.Line2D at 0x7f0b5f299828>]

minimizing MSE

In [22]:
np.random.seed(0)
alpha = 0.0005          #learning rate
iterations = 10000
f_i_optim =  np.ones(M) #random starting guess


cost_list = []
for i in range(0, iterations):
    hypothesis = idct(f_i_optim, type=2, norm='ortho')
#     loss = 2*(hypothesis - sum_t_sample_new)
    cost = np.sum((hypothesis - sum_t_sample_new) ** 2) 
    
    if i % 100 == 0: 
#         print("Iteration %d | Cost: %f" % (i, cost))
        cost_list.append(cost)
    
    gradient = idct(2*(hypothesis - sum_t_sample_new), type=2, norm='ortho')  # from derivation
    
    f_i_optim = f_i_optim - alpha * gradient
#     break

plt.plot(range(M),f_i_dash)
plt.plot(range(M),f_i_optim)
# print(f_i_optim, f_i_dash, (f_i_dash - f_i_optim)/f_i_dash*100 )

plt.figure(1)
plt.subplot(121)
plt.plot(range(M),f_i_dash, color = 'orange')
plt.plot(range(M),f_i_optim, linestyle="", marker="o", color = 'blue')

plt.subplot(122)
plt.plot(cost_list)

print('min cost:', min(cost_list))

min cost: 0.005628660173090521

minimizing MSE with L2 norm

In [21]:
np.random.seed(0)
alpha = 0.0005          #learning rate
iterations = 10000
f_i_optim =  np.ones(M) #random starting guess


cost_list = []
for i in range(0, iterations):
    hypothesis = idct(f_i_optim, type=2, norm='ortho')
#     loss = 2*(hypothesis - sum_t_sample_new)
    cost = np.sum((hypothesis - sum_t_sample_new) ** 2) 
    
    if i % 100 == 0: 
#         print("Iteration %d | Cost: %f" % (i, cost))
        cost_list.append(cost)
    
    gradient = idct(2*(hypothesis - sum_t_sample_new), type=2, norm='ortho') + 2*f_i_optim
    
    f_i_optim = f_i_optim - alpha * gradient
#     break

plt.plot(range(M),f_i_dash)
plt.plot(range(M),f_i_optim)
# print(f_i_optim, f_i_dash, (f_i_dash - f_i_optim)/f_i_dash*100 )

plt.figure(1)
plt.subplot(121)
plt.plot(range(M),f_i_dash, color = 'orange')
plt.plot(range(M),f_i_optim, linestyle="", marker="o", color = 'blue')

plt.subplot(122)
plt.plot(cost_list)

print('min cost:', min(cost_list))

min cost: 28.202160800110388

minimizing MSE with L1 norm

In [23]:
np.random.seed(0)
alpha = 0.0005          #learning rate
iterations = 10000
f_i_optim =  np.ones(M) #random starting guess


cost_list = []
for i in range(0, iterations):
    hypothesis = idct(f_i_optim, type=2, norm='ortho')
#     loss = 2*(hypothesis - sum_t_sample_new)
    cost = np.sum((hypothesis - sum_t_sample_new) ** 2) 
    
    if i % 100 == 0: 
#         print("Iteration %d | Cost: %f" % (i, cost))
        cost_list.append(cost)
    
    gradient = idct(2*(hypothesis - sum_t_sample_new), type=2, norm='ortho') + abs(f_i_optim)/f_i_optim
    
    f_i_optim = f_i_optim - alpha * gradient
#     break

plt.plot(range(M),f_i_dash)
plt.plot(range(M),f_i_optim)
# print(f_i_optim, f_i_dash, (f_i_dash - f_i_optim)/f_i_dash*100 )

plt.figure(1)
plt.subplot(121)
plt.plot(range(M),f_i_dash, color = 'orange')
plt.plot(range(M),f_i_optim, linestyle="", marker="o", color = 'blue')

plt.subplot(122)
plt.plot(cost_list)

print('min cost:', min(cost_list))

min cost: 17.71874674153577

MSE sample freq vs error

In [55]:
from scipy.fftpack import ifft, idct, dct

np.random.seed(0)
alpha = 0.0005          #learning rate
iterations = 10000
x = np.arange(0.0, 3.0, 0.0001)
# compute the value (amplitude) of the sin wave at the for each sample
t_697  = np.sin(2*np.pi*697 * (x)) 
t_1209 = np.sin(2*np.pi*1209 * (x)) 
sum_t = t_697 + t_1209

cost_list_orig_min = []
list_M = np.arange(0.1, 5, 0.1)*1000

for M in list_M:
    M = int(M)
#     print(M)
    cost_list = []
    
    sum_t_sample = np.zeros(len(sum_t))
    sum_t_sample[sum_t_sample==0]=np.nan
    sum_t_sample_index = np.random.choice(range(len(sum_t)), M, replace = False)
    sum_t_sample[sum_t_sample_index] = sum_t[sum_t_sample_index]
    sum_t_sample_new = sum_t_sample[np.isfinite(sum_t_sample)]
    f_i_dash  = dct(sum_t_sample_new, type=2, norm='ortho') #this does the same transform as neil's

#     sum_t_sample.shape
    
    f_i_optim =  np.ones(M)
#     print(f_i_optim, M)

    for i in range(0, iterations):
        
        hypothesis = idct(f_i_optim, type=2, norm='ortho')
#         print(hypothesis.shape, sum_t_sample_new.shape)
        cost_orig = np.sum((hypothesis - sum_t_sample_new) ** 2) 
        

        if i % 100 == 0: 
#             print("Iteration %d | Cost: %f" % (i, cost_orig))
            cost_list.append(cost_orig)

        gradient = idct(2*(hypothesis - sum_t_sample_new), type=2, norm='ortho')

        f_i_optim = f_i_optim - alpha * gradient
    #     break

    if M%500 == 0: print('M:', M, 'min cost:', np.mean(cost_list[-5:]) )
    cost_list_orig_min.append(np.mean(cost_list[-5:]))
    
#     break

plt.plot(list_M, cost_list_orig_min)

M: 500 min cost: 0.40094599777774553
M: 1000 min cost: 1.0719564597431832
M: 1500 min cost: 1.6809698724977014
M: 2000 min cost: 3.77033667904295
M: 2500 min cost: 4.137695670056404
M: 3000 min cost: 5.054189335544822
M: 3500 min cost: 5.945624779959034
M: 4000 min cost: 5.937425432053777
M: 4500 min cost: 7.3576733525986855
Out[55]:
[<matplotlib.lines.Line2D at 0x7f60b5fe6160>]

L1 sample freq vs error

In [48]:
from scipy.fftpack import ifft, idct, dct

np.random.seed(0)
alpha = 0.0005          #learning rate
iterations = 10000
x = np.arange(0.0, 3.0, 0.0001)
# compute the value (amplitude) of the sin wave at the for each sample
t_697  = np.sin(2*np.pi*697 * (x)) 
t_1209 = np.sin(2*np.pi*1209 * (x)) 
sum_t = t_697 + t_1209

cost_list_L1_min = []
list_M = np.arange(0.1, 5, 0.1)*1000

for M in list_M:
    M = int(M)
#     print(M)
    cost_list = []
    
    sum_t_sample = np.zeros(len(sum_t))
    sum_t_sample[sum_t_sample==0]=np.nan
    sum_t_sample_index = np.random.choice(range(len(sum_t)), M, replace = False)
    sum_t_sample[sum_t_sample_index] = sum_t[sum_t_sample_index]
    sum_t_sample_new = sum_t_sample[np.isfinite(sum_t_sample)]
    f_i_dash  = dct(sum_t_sample_new, type=2, norm='ortho') #this does the same transform as neil's

#     sum_t_sample.shape
    
    f_i_optim =  np.ones(M)
#     print(f_i_optim, M)

    for i in range(0, iterations):
        
        hypothesis = idct(f_i_optim, type=2, norm='ortho')
#         print(hypothesis.shape, sum_t_sample_new.shape)
        cost = np.sum((hypothesis - sum_t_sample_new) ** 2 + abs(hypothesis) ) 
        

        if i % 100 == 0: 
#             print("Iteration %d | Cost: %f" % (i, cost_orig))
            cost_list.append(cost)

        gradient = idct(2*(hypothesis - sum_t_sample_new), type=2, norm='ortho') + abs(f_i_optim)/f_i_optim

        f_i_optim = f_i_optim - alpha * gradient
    #     break

    if M%500 == 0: print('M:', M, 'min cost:', np.mean(cost_list[-5:]) )
    cost_list_L1_min.append(np.mean(cost_list[-5:]))
    
#     break

plt.plot(list_M, cost_list_L1_min)

M: 500 min cost: 373.80968700996914
M: 1000 min cost: 762.320865241569
M: 1500 min cost: 1143.9742817008282
M: 2000 min cost: 1526.3704609225572
M: 2500 min cost: 1889.7698295366984
M: 3000 min cost: 2315.797069000605
M: 3500 min cost: 2663.654401697786
M: 4000 min cost: 3027.1428522030324
M: 4500 min cost: 3489.4347726195047
Out[48]:
[<matplotlib.lines.Line2D at 0x7f60b5f12b38>]
In [56]:
from scipy.fftpack import ifft, idct, dct

np.random.seed(0)
alpha = 0.0005          #learning rate
iterations = 10000
x = np.arange(0.0, 3.0, 0.0001)
# compute the value (amplitude) of the sin wave at the for each sample
t_697  = np.sin(2*np.pi*697 * (x)) 
t_1209 = np.sin(2*np.pi*1209 * (x)) 
sum_t = t_697 + t_1209

cost_list_L2_min = []
list_M = np.arange(0.1, 5, 0.1)*1000

for M in list_M:
    M = int(M)
#     print(M)
    cost_list = []
    
    sum_t_sample = np.zeros(len(sum_t))
    sum_t_sample[sum_t_sample==0]=np.nan
    sum_t_sample_index = np.random.choice(range(len(sum_t)), M, replace = False)
    sum_t_sample[sum_t_sample_index] = sum_t[sum_t_sample_index]
    sum_t_sample_new = sum_t_sample[np.isfinite(sum_t_sample)]
    f_i_dash  = dct(sum_t_sample_new, type=2, norm='ortho') #this does the same transform as neil's

#     sum_t_sample.shape
    
    f_i_optim =  np.ones(M)
#     print(f_i_optim, M)

    for i in range(0, iterations):
        
        hypothesis = idct(f_i_optim, type=2, norm='ortho')
#         print(hypothesis.shape, sum_t_sample_new.shape)
        cost = np.sum((hypothesis - sum_t_sample_new) ** 2 + hypothesis**2 ) 
        

        if i % 100 == 0: 
#             print("Iteration %d | Cost: %f" % (i, cost_orig))
            cost_list.append(cost)

        gradient = idct(2*(hypothesis - sum_t_sample_new), type=2, norm='ortho') + 2*f_i_optim

        f_i_optim = f_i_optim - alpha * gradient
    #     break

    if M%500 == 0: print('M:', M, 'min cost:', np.mean(cost_list[-5:]) )
    cost_list_L2_min.append(np.mean(cost_list[-5:]))
    
#     break

plt.plot(list_M, cost_list_L2_min)

M: 500 min cost: 277.1616545710751
M: 1000 min cost: 554.4919543045163
M: 1500 min cost: 854.2662888756416
M: 2000 min cost: 1112.676149371171
M: 2500 min cost: 1381.390825488809
M: 3000 min cost: 1686.0580189516154
M: 3500 min cost: 1966.656452153197
M: 4000 min cost: 2214.960680984254
M: 4500 min cost: 2583.2476648734423
Out[56]:
[<matplotlib.lines.Line2D at 0x7f60b5e5f080>]
In [57]:
plt.plot(list_M, cost_list_orig_min, color = 'blue')
plt.plot(list_M, cost_list_L1_min, color = 'orange')
plt.plot(list_M, cost_list_L2_min, color = 'red')
Out[57]:
[<matplotlib.lines.Line2D at 0x7f60b5e8c4a8>]