Untitled

from numpy import *
import scipy as scipy
import scipy.linalg
import math
import pylab

with open('signal.dat', 'r') as file:
    lines = file.readlines()
    lines = [line.rstrip() for line in lines]

points = []
for point in lines:
    points.append(point.split(', '))

X = [ float(points[i][0]) for i in range(len(points)) ]
Y = [ float(points[i][1]) for i in range(len(points)) ]


A = zeros((200,4))
for i,x in enumerate(X):
    a = array([ math.sin(x), math.cos(x), math.sin(2*x), math.cos(2*x) ])
    A[i,:] = a


ATA = A.T@A
b = A.T@Y

Q,R = scipy.linalg.qr(ATA)
print(f'R  = {R}')

def backward_sub(R, b):
    x = zeros((b.size,1))
    n = b.size - 1
    for row in R[::-1][:]:
        x[n] = b[n]
        for i in range(b.size - 1, n, -1):
            x[n] -= x[i]*R[n][i]
        x[n] = x[n]/R[n][n]
        n -= 1
    return x

x = backward_sub(R,Q.T@b)
y = lambda t : x[0,0]*math.sin(t) + x[1,0]*math.cos(t) + x[2,0]*math.sin(2*t) + x[3,0]*math.cos(2*t)

U, s, VT = scipy.linalg.svd(ATA, full_matrices = False)
S = diag(s)

VTx = backward_sub(S, U.T@b)
xSVD = VT.T@VTx

ySVD = lambda t : xSVD[0,0]*math.sin(t) + xSVD[1,0]*math.cos(t) + xSVD[2,0]*math.sin(2*t) + xSVD[3,0]*math.cos(2*t)


yp = [y(t) for t in X]
ypSVD = [ySVD(t) for t in X]
pylab.plot(X,yp)
pylab.plot(X,ypSVD, 'm')
pylab.plot(X,Y, 'co')
pylab.show()

print(f' x = {x}')
print(f' xSVD = {xSVD}')