python code

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

T = 12; #[s] durata
N = 32; # numarul de armonici
Ne = 1000; # numar de esantioana de timp
Te = T/Ne; #[s] perioada de esantionare

# Axa timp
t = np.arange(0,T,Te);
t = np.asmatrix(t);


# Punctele desanate vor fi de coordonate (x,y,z), respectiv (t,n,u)
# unde:
# t - timp (pe x)
# n - armonica corespunzatoare pulsatie n*omega (a n-a pulsatie fundamentala)
# u - semnalul u(t) corespunzator

# Coordonata x corespunzatoare timpului
# (sunt necesare N+1 astfel de "axe" - N armonici + 1 suma lor)

t = t.transpose();
ones_matrix = np.ones((1, N+1), dtype=np.float);
X = np.dot(t,ones_matrix);

# Coordonata y corespunzatoare indicelui armonicii
# (sunt necesare N+1 astfel de "axe" - N armonici + 1 suma lor)

n_axa = np.arange(0,N+1,1);
n_axa = np.asmatrix(n_axa);

ones_matrix = np.ones((Ne, 1), dtype=np.float);
Y = np.dot(ones_matrix,n_axa);

# Coordonata z
# Amplitudinile corespunzatoare
n = np.arange(1,N+1,1);
n = np.asmatrix(n);


pi_n = np.pi * n;
an =(np.cos((4*pi_n)/5)-1+np.cos((7*pi_n)/5)-np.cos(pi_n)+np.cos(2*pi_n)-np.cos((9*pi_n)/5))/(-pi_n);
bn =(np.sin((4*pi_n)/5)+np.sin((7*pi_n)/5)-np.sin(pi_n)+np.sin(2*pi_n)-np.sin((9*pi_n)/5))/pi_n;

ones_matrix = np.ones((Ne, 1), dtype=np.float);
An = np.dot(ones_matrix,an);
Bn = np.dot(ones_matrix,bn);


# matricea de fecvente ce se aplica
w = 2*np.pi/T;
Y_part = Y[:,1:(N+1)];
W_n=  np.multiply(w,Y_part);
X_part = X[:,1:(N+1)];
fi = np.multiply(W_n,X_part);

Z = X; # initializare
Z[:,1:(N+1)] = np.multiply(An,np.sin(fi)) # + np.multiply(Bn,np.cos(fi));

# suma
ones_matrix = np.ones((N,1), dtype=np.float);
Z_part = Z[:,1:(N+1)];
Z[:,0] = np.dot(Z_part,ones_matrix);


ones_matrix =  np.ones((1,N), dtype=np.float);
one_array = [1, 1];
one_array = np.asmatrix(one_array);
one_array = one_array.transpose();
cx = (T+Te) * one_array * ones_matrix;
one_matrix = [[1] , [1]];
cy = np.dot(one_matrix,n);
zeros_matrix = np.zeros((1,N), dtype=np.float);
abs_an = np.abs(an);
cz = np.concatenate((zeros_matrix,abs_an), axis = 0);


fig = plt.figure(figsize=(9,10));
ax = fig.add_subplot(111, projection='3d');
ax.plot_wireframe(X[:,1:(N+1)],Y[:,1:(N+1)],Z[:,1:(N+1)], label = '',color="black");
ax.plot_wireframe(X[:,1],Y[:,1],Z[:,1], label = '',color="blue");
ax.legend();
ax.set_xlabel('Timp');
ax.set_ylabel('Indice armonica n ( n*2*pi*f = frecventa)');
ax.set_zlabel('Amplitudine');
ax.grid();
ax.plot_wireframe(cx,cy,cz,color="red");
plt.show();