I have the following code which is part of a practice I'm working on: Here is a function that obtains 2*N+1 coefficients of the Fourier series for the input function "ft" written in terms of sym.Heaviside(t), the goal is to get from this function only the numerical values so it is possible to perform operations with the obtained Power Spectral Density (PSD):
def ck_power(ft,f,N=5):
c_k_coeff=f*sym.integrate(ft*sym.exp(-2*j*np.pi*t*k*f),(t,-1/(2*f),1/(2*f)))
#C_k_coeff=sym.lambdify(k,c_k_coeff,modules=['numpy'])
#PSD=[abs((c_k_coeff(i)))**2 for i in range(-N,N+1)]
##PSD=[(c_k_coeff.subs(k,i))*(c_k_coeff.subs(k,-i)) for i in range(-N,N+1)]
### "worked" but keeps the same behavior. PSD=[(c_k_coeff.evalf(subs={k:i}))*(c_k_coeff.evalf(subs={k:-i})) for i in range(-N,N+1)]
PSD=[(c_k_coeff.evalf(subs={k:i}))*(c_k_coeff.evalf(subs={k:-i})) for i in range(-N,N+1)]
PSD=[sym.N(abs(i.subs({u(0):0.5})).evalf(3)) for i in PSD]
positive_PSD=PSD[len(PSD)//2:]
# because the k=0 c_k isn't twice in the k axis apporting power:
positive_PSD[1:]=[2*i for i in positive_PSD[1:]]
PSD=np.array(PSD)
positive_PSD=np.array(positive_PSD)
return PSD,positive_PSD
The problem is that when I use the function with the code below the figures I'm still having a weird piecewise function as the output (especially with the sawtooth wave)
The other functions don't have the same problem
But the sawtooth shows directly the piecewise output
these are the output expressions for k=0:
#vertical
type_hash=("Frequency","Power")
frequency_hash=("100KHz","1MHz","100MHz")
frequency=(100e+3,1e+6,100e+6)
waveform_hash=("Pulsetrain","Sawtooth","Triangle")
#horizontal
armonic=list(range(6))
signal={
"Pulsetrain":my.pulsetrain,
"Sawtooth":my.sawtooth,
"Triangle":my.triangle
}
t1={"waveform":[],"frequency":[],"type\\armonic":[]}
for i in armonic:
t1[str(i)]=[]
(A,D,Aoff,toff)=(1,0.5,0,0)
for iwav in waveform_hash:
for ifreq in range(len(frequency)):
ftemp=frequency[ifreq]
for itype in range(len(type_hash)):
f_now=signal[iwav](A,1/ftemp,D,Aoff,toff)
PSD,pPSD=my.ck_power(f_now,ftemp,N=5) #<<--right here is where the functions is called
for i in armonic:
if not itype:
t1[str(i)].append(f"{ftemp*i:.2e}")
else:
t1[str(i)].append(pPSD[i])# add here Power
t1["type\\armonic"].append(type_hash[itype])
t1["frequency"].append(frequency_hash[ifreq])
t1["waveform"].append(iwav)
d1=pd.DataFrame(t1)
d1
these are my imports:
import numpy as np
import matplotlib.pyplot as plt
from scipy.io.wavfile import read
from scipy import signal
import sympy as sym
import sigplot as my
import pandas as pd
Lastly, this is part of my sigplot.py file
import numpy as np
import matplotlib.pyplot as plt
from scipy.io.wavfile import read
from scipy import signal
import sympy as sym
from sympy import Heaviside as u
from IPython.display import display, Math
import pandas as pd
t=sym.symbols('t')
k=sym.symbols('k', integer=True)
j=sym.I
#Fourier
def ck(ft,T):
c_k=1/T*sym.integrate(ft*sym.exp(-j*k*(2*np.pi/T)*t),(t,-T/2,T/2))
c_0=c_k.subs(k,0)
#c_k_mag =[sym.simplify((abs(c_n_coeff.subs(n,i)))) for i in range (-N,N+1)]
return (c_k,c_0)
def espectral(c_k,f,N):
PSD=[abs((c_k.subs(k,i)))**2 for i in range(-N,N+1)]
positive_PSD=[2*abs((c_k.subs(k,i)))**2 if i!=0 else (c_k.subs(k,0))**2 for i in range(N+1)]
return PSD,positive_PSD
def ck_power(ft,f,N=5):
c_k_coeff=f*sym.integrate(ft*sym.exp(-2*j*np.pi*t*k*f),(t,-1/(2*f),1/(2*f)))
#C_k_coeff=sym.lambdify(k,c_k_coeff,modules=['numpy'])
#PSD=[abs((c_k_coeff(i)))**2 for i in range(-N,N+1)]
##PSD=[(c_k_coeff.subs(k,i))*(c_k_coeff.subs(k,-i)) for i in range(-N,N+1)]
### worked. PSD=[(c_k_coeff.evalf(subs={k:i}))*(c_k_coeff.evalf(subs={k:-i})) for i in range(-N,N+1)]
PSD=[(c_k_coeff.evalf(subs={k:i}))*(c_k_coeff.evalf(subs={k:-i})) for i in range(-N,N+1)]
PSD=[sym.N(abs(i.subs({u(0):0.5})).evalf(3)) for i in PSD]
positive_PSD=PSD[len(PSD)//2:]
# because the k=0 c_k isn't twice in the k axis apporting power.
positive_PSD[1:]=[2*i for i in positive_PSD[1:]]
PSD=np.array(PSD)
positive_PSD=np.array(positive_PSD)
return PSD,positive_PSD
def espectrump(c_k,f,N):
# this function should not be used for performance purposes
PSD=[abs((c_k.subs(k,i)))**2 for i in range(-N,N+1)]
positive_PSD=[2*abs((c_k.subs(k,i)))**2 if i!=0 else (c_k.subs(k,0))**2 for i in range(N+1)]
plt.stem(list(range(-N,N+1)),PSD,markerfmt='^',use_line_collection=True)
plt.xlabel("$k$")
plt.grid()
plt.show()
return positive_PSD
def psdplot(N,PSD):
plt.stem(list(range(-N,N+1)),PSD,markerfmt='^',use_line_collection=True)
plt.xlabel("$k$")
plt.ylabel("$y=10{log()}$")
plt.grid()
plt.show()
return None #excuse me, I can't hide it. Is that simple
#waveforms
def pulsetrain(A,T,D,off=0,toff=0):
f=A*(-u(t+T/2)+2*(u(t+D*T/2)-u(t-D*T/2))+u(t-T/2))
foff=f+off*(u(t+T/2)-u(t-T/2))
foff=foff.subs(t,t-toff)
return foff
def triangle(A,T,D=None,off=0,toff=0):
m = (2*A)/T
toff_temp=0.5
tri=m*t*u(t+T/2)-2*m*t*u(t)+m*t*u(t-T/2)
foff = tri+(off+toff_temp)*(u(t+T/2)-u(t-T/2))
foff = foff.subs(t,t-toff)
return foff
def sawtooth(A,T,D=None,off=0,toff=0):
m = A/T
saw =m*t*u(t+T/2)-m*t*u(t-T/2)-(m/2)*u(t-T/2)
foff = saw+off*(u(t+T/2)-u(t-T/2))
foff = foff.subs(t,t-toff)
return foff
Thanks in advance!
