# Untitled

unknown
python
a year ago
1.4 kB
1
Indexable
Never
```from numpy import *
import scipy as scipy
import scipy.linalg
import math
import pylab

with open('signal.dat', 'r') as file:
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}')
```