如何逐段计算类似于sympy的表达式，以便轻松获得仅数值的输出？

Question

我有以下代码，这是我正在进行的实践的一部分： 这里是一个函数，该函数获得输入函数“ft”的傅里叶级数的2*N+1系数，输入函数“ft”用sym.Heaviside（t）表示，目标是从该函数中仅获得数值，因此可以使用获得的功率谱密度（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

问题是，当我使用数字下面的代码的函数时，我仍然有一个奇怪的分段函数作为输出（特别是锯齿波）

其他函数没有相同的问题

但是锯齿波直接显示了分段输出

以下是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

以下是我的作品：

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

最后，这是我的sigplot.py文件的一部分

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

提前谢谢