Как оценить симпатичные кусочно-подобные выражения, чтобы было легко получить только числовой вывод? - PullRequest
Новая
       25

Как оценить симпатичные кусочно-подобные выражения, чтобы было легко получить только числовой вывод?

0 голосов
/ 12 марта 2020

У меня есть следующий код, который является частью практики, над которой я работаю: вот функция, которая получает 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

Проблема заключается в том, что при использовании функция с кодом ниже рисунков У меня все еще есть странная кусочная функция в качестве выхода (особенно с пилообразной волной)

The other functions apparently doesn't have the same problem Другие функции не имеют того же задача

But the sawtooth shows directly the piecewise output Но пилообразная схема показывает кусочно-прямой вывод

это выходные выражения для k = 0: description

#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

Заранее спасибо!

...