Signal Processing Asked on February 8, 2021
Are there any established techniques for estimating the envelope of a swept filter? I tried a first order piece-wise polynomial approximation using the SHGO global minimizer algorithm (Python code below), but it fails. Setting the points by hand approximately succeeds.
The purpose for the estimation is to obtain a synthesizer’s filter decay envelope, after the filter is known. The input is thereby one of the typical waveforms a synthesizer can generate, in the code a sawtooth waveform was chosen.
import numpy
import scipy.optimize
import scipy.signal
#2nd order low pass filter
def filt(X,F,X0=0,X1=0,Y0=0,Y1=0,resonance=1):
X=numpy.append(X,[X1,X0])
Y=numpy.zeros(len(F))
Y=numpy.append(Y,[Y1,Y0])
for i in range(len(F)):
wd=F[i]*2*numpy.pi
T=1/88200
d=resonance
C=numpy.tan(wd*T/2)
D=(1+d*C+C**2)
a1=2*(C**2-1)/D
a2=(1-d*C+C**2)/D
b0=C**2/D
b1=2*b0
b2=b0
Y[i]=-a1*Y[i-1]-a2*Y[i-2]+b0*X[i]+b1*X[i-1]+b2*X[i-2]
return Y[0:len(Y)-2]
#X = filter input, Y = Output.
freq=600
t=numpy.linspace(0, 1, 88200)
X=scipy.signal.sawtooth(2 * numpy.pi * freq * t)
F=freq+5*freq*numpy.e**(-t*40) #Envelope
Y=filt(X,F)
Y+=numpy.random.normal(0,0.075,len(Y))
#Estimation (Result in Fest):
def PolMin(X,P,I,X0,X1,Y0,Y1):
F=numpy.zeros(I[1]-I[0])
F=numpy.polyval(P,range(I[1]-I[0]))
Y=filt(X[I[0]:I[1]],F,X0,X1,Y0,Y1)
return Y
def Target(P,I,X0,X1,Y0,Y1):
return numpy.sum((PolMin(X,P,I,X0,X1,Y0,Y1)-Y[I[0]:I[1]])**2)
chunkl=200
length=25000
chunks=int(25000/200)
Results = numpy.zeros([chunks,2])
B=[(-2,0),(0,10000)]
Y0=0
Y1=0
X0=0
X1=0
Z=numpy.zeros(length+1)
Fest=numpy.zeros(length+1)
for i in range(chunks):
chunkst=i*chunkl
chunke=chunkl*(i+1)+1
f=lambda x: Target(x,(chunkst,chunke),X0,X1,Y0,Y1)
Results[i,:]=scipy.optimize.shgo(f,B,iters=4).x
endpoint=numpy.polyval(Results[i,:],chunkl*(i+1))
B=[(-1,1),(endpoint,endpoint)]
Fest[chunkst:chunke]=numpy.polyval(Results[i,:],range(chunkl+1))
Z[chunkst:chunke]=filt(X[chunkst:chunke],Fest[chunkst:chunke],X0,X1,Y0,Y1)
Y0=Z[chunke-2]
Y1=Z[chunke-3]
X0=X[chunke-2]
X1=X[chunke-3]
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP