带有步响应的python卷积 [英] python- convolution with step response

问题描述

我想计算这个整数$ \ frac {1} {L} \ int _ {-\ infty} ^ {t} H(t ^ {'})\ exp(-\ frac {R} {L}(tt ^ {'}))dt ^ {'} $使用numpy.convolution，其中$ H(t)$是heavside函数.我应该得到等于$ \ exp(-\ frac {R} {L} t)H(t)$以下是我所做的我通过将变量乘以不同的H(t)来将积分的限制从-inf更改为+ inf，然后将其用作与H(t)(积分中的那个)进行卷积的函数，但是输出图绝对不是exp函数，我也找不到我的代码中的任何错误，请帮忙，任何提示或建议将不胜感激！

I want to compute this integral $\frac{1}{L}\int_{-\infty}^{t}H(t^{'})\exp(-\frac{R}{L}(t-t^{'}))dt^{'}$ using numpy.convolution, where $H(t)$ is heavside function. I am supposed to get this equals to $\exp(-\frac{R}{L}t)H(t)$ below is what I did, I changed the limitation of the integral into -inf to +inf by change of variable multiplying a different H(t) then I used this as my function to convolve with H(t)(the one inside the integral), but the output plot is definitely not a exp function, neither I could find any mistakes in my code, please help, any hint or suggestions will be appreciated!

import numpy as np
import matplotlib.pyplot as plt
R = 1e3 
L = 3. 

delta = 1
Nf = 100
Nw = 200
k = np.arange(0,Nw,delta)
dt = 0.1e-3 
tk = k*dt
Ng = Nf + Nw -2
n = np.arange(0,Nf+Nw-1,delta)
tn = n*dt

#define H
def H(n):
    H = np.ones(n)
    H[0] = 0.5
    return H

#build ftns that get convoluted
f = H(Nf)
w = np.exp((-R/L)*tk)*H(Nw)

#return the value of I
It = np.convolve(w,f)/L

#return the value of Voutput, b(t)
b = H(Ng+1) - R*It
plt.plot(tn,b,'o')
plt.show()

推荐答案

您的代码所涉及的问题与其说是概念上的，还不如说是编程.将卷积重写为Integral [HeavisideTheta [t-t'] * Exp [-R/L * t']，-Inf，t](这是Mathematica代码)，经检查，您发现H(t-t')始终为在限制内为1(t'= t处除外，这是积分限制...但这并不重要).因此，实际上，您实际上并没有执行完整的卷积...您基本上只是在进行卷积的一半(或三分之一).

The issue with your code is not so much programming as it is conceptual. Rewrite the convolution as Integral[HeavisideTheta[t-t']*Exp[-R/L * t'], -Inf, t] (that's Mathematica code) and upon inspection you find that H(t-t') is always 1 within the limits (except for at t'=t which is the integration limit... but that's not important). So in reality you're not actually performing a complete convolution... you're basically just taking half (or a third) of the convolution.

如果您认为卷积是反转一个序列，然后一次移动一个移位并将其全部累加(请参见 http://en.wikipedia.org/wiki/Convolution#Derivations -卷积的可视化解释)，那么您想要的是中间的一半……即仅当它们重叠时.您不希望引入(第4个图形): http://zh.wikipedia.org/wiki/File:Convolution3.svg ).您确实需要导出.

If you think of a convolution as inverting one sequence and then going one shift at the time and adding it all up (see http://en.wikipedia.org/wiki/Convolution#Derivations - Visual Explanation of Convolution) then what you want is the middle half... i.e. only when they're overlapping. You don't want the lead-in (4-th graph down: http://en.wikipedia.org/wiki/File:Convolution3.svg). You do want the lead-out.

现在，最简单的代码修复方法如下:

Now the easiest way to fix your code is as such:

#build ftns that get convoluted
f = H(Nf)
w = np.exp((-R/L)*tk)*H(Nw)

#return the value of I
It = np.convolve(w,f)/L
max_ind = np.argmax(It)
print max_ind
It1 = It[max_ind:]

导入是卷积积分(在本例中为技术总和)增加的唯一时间...因此，在导入完成后，卷积积分遵循Exp [-x] ...所以告诉python只在达到最大值后才取值.

The lead-in is the only time when the convolution integral (technically sum in our case) increases... thus after the lead-in is finished the convolution integral follows Exp[-x]... so you tell python to only take values after the maximum is achieved.

#返回Voutput的值，b(t)现在可以正常工作了！

#return the value of Voutput, b(t) works perfectly now!

注意:由于需要引出线，因此不能使用 np.convolve(a，b，mode ='valid').

Note: Since you need the lead-out you can't use np.convolve(a,b, mode = 'valid').

因此It1看起来像:

So It1 looks like:

b(t)如下:

您不可能获得exp(-x)作为一般形式，因为b(t)的方程是由1-R * exp(-x)给出的...它在数学上不能遵循exp(-x)形式.此时，有3件事:

There is no way you can ever get exp(-x) as the general form because the equation for b(t) is given by 1 - R*exp(-x)... It can't mathematically follow an exp(-x) form. At this point there are 3 things:

1. 单位真的没有意义...检查它们.Heaviside函数为1，R * It1约为10,000.我不确定这是一个问题，但是以防万一，归一化曲线看起来像这样:

1. The units don't really make sense... check them. The Heaviside function is 1 and R*It1 is about 10,000. I'm not sure this is an issue but just in case, the normalized curve looks as such:

如果使用b(t)= R * It1-H(t)，则可以得到exp(-x)形式...此处的代码在这里(您可能必须根据自己的情况进行规范化需求):

You can get an exp(-x) form if you use b(t) = R*It1 - H(t)... the code for that is here (You might have to normalize depending on your needs):

b = R*It1 - H(len(It1))
# print len(tn)
plt.plot(tn[:len(b)], b,'o')
plt.show()

情节如下:

1. 您的问题可能仍未解决，在这种情况下，您需要解释您到底认为错了什么.有了您给我的信息，除非b(t)的方程式混乱，否则b(t)永远不会具有Exp [-x]形式.正如您原始代码中所代表的那样，It1的形式遵循Exp [-x]，但b(t)不能.

这篇关于带有步响应的python卷积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！