<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Foruier Series Approximation for discrete time serires
Dependancy
Step1:  Final Code
Step2:  How and Why it work
We start with the foruier series
Step3:  we start with f that is a sawthoot wave
Step4:  we choose the number of armonichs that we want in this case 50 but $ n \in [1,N)$ where $N$ is len(f)
Step5:  first we get rid of the mean
Step6:  we transforme the array into a coloumn vector
Step7:  we create the C matrix
Step8:  the QR decomposition
Step9:  we truncate the matrix so that it's possible to solve the system
Step10:  we solve the system and get the coeff vector
Step11:  We separate the matrix into the sin and cos coeff list
Step12:  we convert to $c$ and $\phi$ and now we have all the coefficent ready to go
Step13:  we calculate the result function and re-add the mean to the funciton
Step14:  Result
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from scipy import linalg
Explanation: Foruier Series Approximation for discrete time serires
Dependancy:
End of explanation
#the function to calculate the coefficent
def approx(x, y, n, w):
    x = np.matrix(x).transpose()
    y = np.matrix(y).transpose()
    f, b = x.shape
    c, d = y.shape
    if c != f or b != d:
        print('The Input vector have wrong dimension')
        return -1
    j = np.matrix(np.arange(1,n+1))
    V1 = np.cos(w*x*j)
    j = np.matrix(np.arange(n+1,2*n+1))
    V2 = np.sin(w*x*(j-n))
    V = np.concatenate([V1,V2],axis=1)
    Q, R = linalg.qr(V)
    R = R[:2 * n, :2 * n]
    Q = Q[:f, :2 * n]
    # coeff = linalg.solve_triangular(R, (np.dot(Q.transpose(), y)),check_finite=False)
    coeff = linalg.solve_triangular(R, (np.dot(Q.transpose(), y)))
    n = int(len(coeff) / 2)
    mag = np.sqrt(coeff[:n]**2+coeff[n:]**2)
    angle = np.arctan2(coeff[:n],coeff[n:])
    r = []
    for i,(m,a) in enumerate(zip(mag,angle)):
        r.append([float(m),i+1,float(a)])
    return r
#the function to calculate the reconstructed function from the coefficent
def calc_fourier(X,coeff,vmed,w=0.5):
    y = np.zeros_like(X) + vmed
    for (m,i,p) in coeff:
            y += m*np.sin(w*i*X+p)
    return y
#approx a function and get both the coefficent and the reconstructed function
def fourier_approx(funzione,n=0,w=0.5):
    fmean = np.mean(funzione)
    funzione = list(funzione)
    funzione = funzione + funzione[::-1] 
    mean = np.mean(funzione)
    funzione = [z - mean for z in funzione]
    T = np.linspace( 0, 4 * np.pi, num=len(funzione), endpoint=True)
    if n == 0:
        n = int(len(T) / 2) - 1
    if n < 1:
        return -1
    coeff = approx(T, funzione, n, w)
    T =  np.array(T [:int(len(T)/2)])
    funzione =  np.array(funzione [:int(len(funzione)/2)])
    
    y = calc_fourier(T,coeff[:index+1],0)
    
    return y,coeff
Explanation: Final Code
End of explanation
import matplotlib.pyplot as plt
from pprint import pprint
Explanation: How and Why it work
We start with the foruier series:
$$f(t) = \frac{a_0}{2} + \sum_{n = 1}^{N} \left [ a_n cos(nwt) + b_n sin(nwt) \right ]$$
we get if we manage to get a function with mean value equal to zero we can simplify to:
$$f(t) = \sum_{n = 1}^{N} \left [ a_n cos(nwt) + b_n sin(nwt) \right ]$$
$$ f(t) = \sum_{n = 1}^{N} < \begin{pmatrix}cos(nwt)\ sin(nwt)\end{pmatrix},\begin{pmatrix}a_n\ b_n\end{pmatrix} >$$
$$f(t) = \sum_{n = 1}^{N} \begin{pmatrix}cos(nwt)& sin(nwt)\end{pmatrix}\begin{pmatrix}a_n\ b_n\end{pmatrix}$$
$$f(t) = \begin{pmatrix}
cos(wt)
&
cos(2wt)
&
...
&
cos((N-1)wt)
& 
sin(wt)
&
sin(2wt)
&
...
&
sin((N-1)wt)
\end{pmatrix}
\begin{pmatrix}
a_1
\
a_2
\
...
\
a_{N-1}
\ 
b_1
\
b_2
\
...
\
b_{N-1}
\end{pmatrix}
$$
now if we substitute t with the succession $T_k = {t_k} $ we get $k$ row forming a matrix $A$ with $k$ row and $2n$ column
and the succession $Y_k(T_k) = {y_k}$ that$ y_k = f(t_k)$ wich form the vector $Y$ 
$$ f(T_k) = \begin{pmatrix}
cos(wt_1)
&
cos(2wt_1)
&
...
&
cos((N-1)wt_1)
& 
sin(wt_1)
&
sin(2wt_1)
&
...
&
sin((N-1)wt_1)
\
cos(wt_2)
&
cos(2wt_2)
&
...
&
cos((N-1)wt_2)
& 
sin(wt_2)
&
sin(2wt_2)
&
...
&
sin((N-1)wt_2)
\
...
&
...
&
...
&
...
&
...
&
...
&
...
&
...
\
cos(wt_K)
&
cos(2wt_K)
&
...
&
cos((N-1)wt_K)
& 
sin(wt_K)
&
sin(2wt_K)
&
...
&
sin((N-1)wt_K)
\end{pmatrix}
\begin{pmatrix}
a_1
\
a_2
\
...
\
a_{N-1}
\ 
b_1
\
b_2
\
...
\
b_{N-1}
\end{pmatrix}
=
\begin{pmatrix}
y_1
\
y_2
\
...
\
...
\ 
...
\
...
\
...
\
y_{k}
\end{pmatrix}
$$
so now we have a linear sistem to solve for+the vector of $a_n$ and $ b_n$ called $C$
$$ A C = Y$$
the system to be solvable it have to have $det(A) > 0$
so we have to have more lineary indipendent row than coloumn, so the maximum number of coloumn it can have is when $k = 2n$ if all the row are lineary independent
in the case $ non lineary independent row < 2n$ the sistem is overdeterminated so to solve it we can use the $QR$ decomposition
so that $ A = QR$ so we can compute $$ A = \frac{R}{Q^T Y} $$
and now we have the coefficent of the fourier series ready to be used
There is a possible ottimization if the approximation is computed in preprocessing:
        if we use the geometric identity
$$ acos(wt) + bsin(wt) = \sqrt{a^2 + b^2} sin(wt + arctan(\frac{b}{a}))$$
and now we can call $c_n = \sqrt{a_n^2 + b_n^2} $ and $\phi_n = arctan(\frac{b_n}{a_n})$
so we have a final series in the form of
$$ f(t) = \sum_{n = 1}^{N} c_n sin(nwt + \phi_n)$$
Implementation
we'll import matplotlib to plot the result
End of explanation
f = np.linspace(0,1,100)
f = np.concatenate((f,f,f,f))
T = np.linspace( 0, 2 * np.pi, num=len(f), endpoint=True) - np.pi
plt.plot(T,f)
plt.show()
Explanation: we start with f that is a sawthoot wave
End of explanation
n = 50
w = 1
Explanation: we choose the number of armonichs that we want in this case 50 but $ n \in [1,N)$ where $N$ is len(f)
End of explanation
f_mean = np.mean(f)
f -= f_mean
plt.plot(T,f)
plt.show()
Explanation: first we get rid of the mean
End of explanation
x = np.matrix(T).transpose()
y = np.matrix(f).transpose()
Explanation: we transforme the array into a coloumn vector
End of explanation
j = np.matrix(np.arange(1,n+1))
C1 = np.cos(w*x*j)
j = np.matrix(np.arange(n+1,2*n+1))
C2 = np.sin(w*x*(j-n))
C = np.concatenate([C1,C2],axis=1)
print(C)
Explanation: we create the C matrix
End of explanation
Q, R = linalg.qr(C)
Explanation: the QR decomposition
End of explanation
R = R[:2 * n, :2 * n]
Q = Q[:x.shape[0], :2 * n]
Explanation: we truncate the matrix so that it's possible to solve the system
End of explanation
coeff = linalg.solve_triangular(R, (np.dot(Q.transpose(), y)))
# coeff = linalg.solve_triangular(R, (np.dot(Q.transpose(), y)),check_finite=False) Alternative way
Explanation: we solve the system and get the coeff vector
End of explanation
n_ = int(len(coeff) / 2)
sin_coeff = coeff[:n_]
cos_coeff = coeff[n_:]
Explanation: We separate the matrix into the sin and cos coeff list
End of explanation
mag = np.sqrt(cos_coeff**2+sin_coeff**2)
phi = np.arctan2(sin_coeff,cos_coeff)
Explanation: we convert to $c$ and $\phi$ and now we have all the coefficent ready to go
End of explanation
y = np.zeros_like(T)
for (m,i,p) in zip(mag,range(n),phi):
        y += m*np.sin(w*(i+1)*T+p)
y += f_mean
f += f_mean
Explanation: we calculate the result function and re-add the mean to the funciton
End of explanation
plt.plot(T,f)
plt.plot(T,y)
plt.show()
Explanation: Result
End of explanation
<END_TASK>