数字信号处理8：利用Python进行数字信号处理基础

我前两天买了本MATLAB信号处理，但是很无语，感觉自己对MATLAB的语法很陌生，看了半天也觉得自己写不出来，所以就对着MATLAB自己去写用Python进行的数字信号处理基础，我写了两天左右，基本上把matlab书上的代码全部用Python实现了，所以，今天贴的代码和图有些多，

要用到的包：

1、Scipy包：其中signal库，这个库是真的绝，很多信号处理的基础函数都有的，

2、numpy包：numpy包中也有很多进行信号处理的，比如说相关、卷积，都有相关函数

3、mmatplotlib包：这就不多说了，信号处理就是用它来展示的，这里主要用到的就是stem方法。

signal库我找了一下，csdn有个博主总结的很全，这是他的博客链接，可以看一看：

(1条消息) scipy.signal信号处理的库（笔记06）_scipy.signalku_月疯的博客-CSDN博客

然后，还可以看一下scipy官方的文档，里面也有很详细的介绍：

Signal processing (scipy.signal) — SciPy v1.10.1 Manual

这里我做了个目录，可以查看相应的方法：

目录

1、离散时间信号序列的表示：

 2、采样定理：

 3、简单离散信号序列：

 4、单位阶跃函数：

5、正弦信号序列：

6、实指数序列：

 7、复指数序列：

8、序列值累加和乘积：

9、序列反转、移位：

 10、信号的尺度变换：

 11、连续信号的奇偶分解：

12、奇函数和偶函数合并：

13、信号的微积分：

 14、积分：

15、卷积运算和相关计算

16、产生信号波形：

 17、连续矩形周期信号采样变成离散信号：

 18、随机函数，这个用numpy就可以直接生成：

 19、三角波，用signal的sawtooth：

 20、sinc曲线:

21、生成非周期三角波

 22、高斯脉冲的实现：

 23、脉冲序列发生器：

24、产生非周期方波：

25、连续时间信号的时域分析

（1）零状态响应：

 （2）冲激响应和阶跃响应：

 （3）各种信号的响应：

（4）连续时间信号的卷积：

26、离散时间系统：

27、离散时间系统的冲激和阶跃响应：

 28、卷积和运算，这都不用说啥了，前面都说过了：

 29、相关序列（自相关、互相关）：

1、离散时间信号序列的表示：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
N=np.linspace(-3,11,15,dtype=int)
x=np.array([0,2,3,3,2,3,0,-1,-2,-3,-4,-5,1,2,1])
dt=0.01
n=N*dt
fig=plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.stem(N,x)
ax2=fig.add_subplot(2,1,2)
ax2.plot(n,x)
ax2.plot(n,np.zeros(len(n)))

 2、采样定理：

#用10hz的采样频率采样
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
dt=0.01
n=np.linspace(0,89,90,dtype=int)
t=n*dt
f=10
x=np.sin(3*np.pi*f*t+0.5)#原始信号
dt=0.1
n=np.linspace(0,9,10,dtype=int)
t1=n*dt
x1=np.sin(3*f*np.pi*t1+0.5)# 采样后的信号
fig1=plt.figure()
ax1=fig1.add_subplot(3,1,1)
ax1.plot(t,x)
ax1.set_title('origional signal')
ax2=fig1.add_subplot(3,1,2)
ax2.plot(t,x)
ax2.plot(t1,x1,'*')
ax2.set_title('Sampling process')
ax3=fig1.add_subplot(3,1,3)
ax3.plot(t1,x1)
ax3.set_title('Sampling signal')
plt.show()

 3、简单离散信号序列：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
n=50
x=np.zeros(n)
x[1]=1
xn=np.linspace(0,n-1,n,dtype=int)
fig1=plt.figure()
ax1=fig1.add_subplot(2,1,1)
ax1.stem(xn,x)
x[1]=0
x[10]=1
ax2=fig1.add_subplot(2,1,2)
ax2.stem(xn,x)
plt.show()

 4、单位阶跃函数：

这个我没有在signal里找到，但是我自己写了一个：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
def u(n):
    if n>=0:
        r=1
    else:
        r=0
    return r
x=np.linspace(-10,10,21,dtype=int)
y=np.array([u(i)for i in x])
plt.stem(x,y)
plt.show()

 很简单的，最简单的是递增序列，就是一个递增函数就行，比如说上面的x，他就是个递增序列：

plt.stem(x,x,'r')

 蓝色的就是刚才的阶跃序列，

5、正弦信号序列：

这个也很简单，延续上面的代码：

# x=np.linspace(-10,10,21,dtype=int)
# z=np.linspace(-10,10,10000,dtype=float)
plt.stem(x,np.sin(x))
plt.plot(z,np.sin(z),'r')
plt.show()

 可以看到，红色的连续信号和蓝色棉签棒的离散信号。

6、实指数序列：

numpy的power函数可以很简单的实现：

x(n)=a^nu(n)

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
import sympy as sy
n=np.linspace(-10,10,21,dtype=int)
y1=np.power(1.6,n)
y2=np.power(-1.6,n)
y3=np.power(0.9,n)
y4=np.power(-0.9,n)
f=plt.figure()
ax1=f.add_subplot(2,2,1)
ax1.stem(n,y1)
ax2=f.add_subplot(2,2,2)
ax2.stem(n,y2)
ax3=f.add_subplot(2,2,3)
ax3.stem(n,y3)
ax4=f.add_subplot(2,2,4)
ax4.stem(n,y4)
ax1.set_title('1.6^n')
ax2.set_title('(-1.6)^n')
ax3.set_title('0.9^n')
ax4.set_title('(-0.9)^n')
plt.show()

 7、复指数序列：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
import cmath
n=np.linspace(0,50,51,dtype=complex)
A=3
a=-1/9
b=np.pi/5
h=-1/9+np.pi/5j
x=A*np.exp(h*n)
f=plt.figure()
ax1=f.add_subplot(2,2,1)
ax1.stem(n,x.real)
ax2=f.add_subplot(2,2,2)
ax2.stem(n,x.imag)
ax3=f.add_subplot(2,2,3)
ax3.stem(n,abs(x))
ax4=f.add_subplot(2,2,4)
y=np.array([-cmath.phase(i) for i in x])#如果直接用phase的话，和matlab计算的angle是相反的，所以我这里为了和matlab一样，就用了-phase
ax4.stem(n,y)
ax1.set_title('real')
ax2.set_title('imag')
ax3.set_title('abs')
ax4.set_title('angle')
plt.show()

8、序列值累加和乘积：

这个就是常规的numpy操作，没啥意思

9、序列反转、移位：

反转的话，可以用[:：-1]或者也可以用y=np.flipud(x)，两者是一样的效果，信号移位就是，在一个信号序列前面或者后面加一个全零数组，前面加就是延迟，后面加就是提前：

# matlab中有一个函数较fliplr(x)来进行序列反转，但是，python会比这个更简单:两种方法：第一种是用列表反转的形式，第二种是用numpy的flipud函数执行
import numpy as np
import matplotlib.pyplot as plt
nx=np.linspace(-2,5,8,dtype=int)
x=np.linspace(2,9,8,dtype=int)
ny=-np.flipud(nx)
y=np.flipud(x)
print(ny)
print(y)
# ny=-nx[::-1]
# y=x[::-1]
fig=plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.stem(nx,x,'.')
ax1.axis([-6,6,-1,9])
ax1.grid(visible=True)
ax1.set_xlabel('n')
ax1.set_ylabel('x(n)')
ax1.set_title('origional sequence')
ax2=fig.add_subplot(2,1,2)
ax2.stem(ny,y,'.')
ax2.axis([-6,6,-1,9])
ax2.grid(visible=True)
ax2.set_xlabel('n')
ax2.set_ylabel('y(n)')
ax2.set_title('Inversion sequence')
plt.show()

#序列位移：
import matplotlib.pyplot as plt
import numpy as np
nx=np.linspace(-2,5,8,dtype=int)
x=np.array([9,8,7,6,5,5,5,5])
y=x
ny1=nx+3
ny2=nx-2
fig=plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.stem(nx,x,'.')
ax1.axis([-5,9,-1,10])
ax1.grid(visible=True)
ax1.set_xlabel('n')
ax1.set_ylabel('x(n)')
ax1.set_title('origional sequence')
ax2=fig.add_subplot(2,2,3)
ax2.stem(ny1,y,'.')
ax2.axis([-5,9,-1,10])
ax2.grid(visible=True)
ax2.set_xlabel('n')
ax2.set_ylabel('y1(n)')
ax2.set_title('y1=(n+3)')
ax3=fig.add_subplot(2,2,4)
ax3.stem(ny2,y,'.')
ax3.axis([-5,9,-1,10])
ax3.grid(visible=True)
ax3.set_xlabel('n')
ax3.set_ylabel('y2(n)')
ax3.set_title('y2=(n-2)')
plt.show()

 10、信号的尺度变换：

把信号的横坐标压缩或者扩展：

t = np.linspace(-4, 4, 8000, dtype=float)
T = 2
f = np.zeros((8000),dtype=float)
t1=2*t
f1=np.zeros((8000),dtype=float)
for i in range(len(t)):
    if ((-1 <= t[i]) & (t[i] <= 1)).any():
        f[i] = 1
for i in range(len(t1)):
    if ((-0.5 <= t1[i]) & (t1[i] <= 0.5)).any():
        f1[i] = 1
fig=plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.plot(t,f)
ax1.axis([-4,4,-0.5,1.5])
ax2=fig.add_subplot(2,1,2)
ax2.plot(t1,f1)
ax2.axis([-4,4,-0.5,1.5])
plt.show()

 11、连续信号的奇偶分解：

我们都知道，一个信号可以分解成一个偶分量和一个奇分量：

f(t)=\frac{1}{2}[f(t)+f(-t)]+\frac{1}{2}[f(t)-f(-t)]

#连续信号的奇偶分解：
#对于一个信号f(n)来说，奇信号：1/2[f(t)-f(-t)]偶信号：1/2[f(t)+f(-t)]
t=np.linspace(-8,8,100000)
f=np.cos(t+1)+t
f1=np.cos(-t+1)-t
g=1/2*(f+f1)
h=1/2*(f-f1)
fig=plt.figure()
ax1=fig.add_subplot(3,1,1)
ax1.plot(t,f)
ax1.set_title('origional signal')
ax2=fig.add_subplot(3,1,2)
ax2.plot(t,g)
ax2.set_title('odd signal')
ax3=fig.add_subplot(3,1,3)
ax3.plot(t,h)
ax3.set_title('even signal')

12、奇函数和偶函数合并：

还是上面的

#将上面的两个奇偶分量合并成原函数：就是反向操作一波，也可以用g-h
t=np.linspace(-8,8,100000)
f=np.cos(t+1)+t
f1=np.cos(-t+1)-t
g=1/2*(f+f1)
h=1/2*(f-f1)
z=g+h
l=h-g
fig=plt.figure()
ax1=fig.add_subplot(4,1,3)
ax1.plot(t,z)
ax1.set_title('origional signal')
ax1=fig.add_subplot(4,1,4)
ax1.plot(t,l)
ax1.set_title('origional signal')
ax2=fig.add_subplot(4,1,1)
ax2.plot(t,g)
ax2.set_title('odd signal')
ax3=fig.add_subplot(4,1,2)
ax3.plot(t,h)
ax3.set_title('even signal')

 多画了一个原函数，将就着看吧

13、信号的微积分：

这个里面有一个heaviside函数，matlab有，但是python没有，然后我就把它自己实现了，要用的可以直接用：

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from scipy import signal
import sympy.plotting as syp
def heaviside(x):
    if x==0:
          r=0.5
    elif x>0:
          r=1
    elif x<0:
          r=0
    return r
#本来要用的，但我们没有用heaviside函数
t=sp.symbols('t')
f=sp.Function('f')(t)
f=t*t+2*t-1
f1=f.diff(t)#一阶导
f2=f.diff(t,t)#二阶导
f3=f.diff(t,t,t)#三阶导
syp.plot(f,f1,f2,f3,(t,-1,2))

 14、积分：

这里用的是：

import sympy as sy
t,f=sy.symbols('t,f')
f=2*t+2
intt=sy.integrate(f,t)
print(intt)
syp.plot(intt,(t,-7,5))

 当然，微积分这里比较水，以为找个很好看的函数很麻烦

15、卷积运算和相关计算

连续信号和离散信号都一样，我们可以用numpy中的convolve函数，singal中也有这个函数，直接计算就可以，这里给例子：y(n)=x(n)*h(n)

#离散时间信号的卷积和运算
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
def uDt(n):
    if n>=0:
        y=1
    else:
        y=0
    return y
nx=np.linspace(-1,5,7,dtype=int,endpoint=True)
nh=np.linspace(-2,10,13,dtype=int,endpoint=True)
nx2=nx-4
nh2=nh-9
x=np.array([uDt(i)for i in nx])-np.array([uDt(i)for i in nx2])
h=np.power(0.9,nh)
h1=np.array([uDt(i)for i in nh])-np.array([uDt(i)for i in nh2])
y=np.convolve(h,h1)
ny1=nx[0]+nh[0]
ny=ny1+np.linspace(0,(len(y)),len(y),dtype=int)
fig=plt.figure()
ax1=fig.add_subplot(3,1,1)
ax1.stem(nx,x)
ax1.grid(visible=True)
ax1.set_title('x(n)')
ax1.axis([-4,16,0,1.5])
ax2=fig.add_subplot(3,1,2)
ax2.stem(nh,h)
ax2.grid(visible=True)
ax2.set_title('h(n))')
ax2.axis([-4,16,0,1.5])
ax3=fig.add_subplot(3,1,3)
ax3.stem(ny,y)
ax3.set_title('y(n)')
ax3.grid(visible=True)
ax3.axis([-4,16,0,9])
plt.show()

 然后相关序列计算，切记，千万不要和person系数混淆了，我先用person那个函数计算来着，后来发现是不对的，signal.correlate和np.correlate都可以完成相关计算：

#计算自相关和互相关
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from scipy.stats import pearsonr
x=np.array([1,3,5,7,9,11,13,15,17,19])
y=np.array([1,1,1,1,2,2,2,2,2,2])
# 计算自相关函数
auto_corr = signal.correlate(x, x, mode='same')
# 计算互相关函数
cross_corr = signal.correlate(x, y, mode='same')
# np.correlate
print("自相关函数：", auto_corr)
print("互相关函数：", cross_corr)

这里就不放图了，没啥意思，就两个类似山峰的曲线，

16、产生信号波形：

from scipy import signal as signal
t=np.linspace(0,1,200,dtype=float)

# scipy.signal.chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True)

h=signal.chirp(t,0,1,120,method='linear',phi=np.pi/3)#linear线性\quadratic二次扫描、logarithmic对数扫描（这时候f0、f1均不能为零）
plt.plot(t,h)
plt.grid(visible=True)
plt.show()

我们这里用singal的chirp函数来生成波形，这里，产生的波形是时间轴为t，时刻0的瞬间频率是f0，时刻t1的瞬间频率是f1，method就是你可以产生波形的方法，phi就是相位，一般来说vertex_zero设置成缺省值就行。

来看一下：

 17、连续矩形周期信号采样变成离散信号：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
f=6000
nt=3
N=15
T=1/f
dt=T/N
n=np.linspace(0,50,51,dtype=int)
tn=n*dt
x=signal.square(2*f*np.pi*tn,duty=0.25)+1
fig=plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.step(tn,x,'k')
ax1.axis([0,nt*T,1.2*min(x),1.1*max(x)])
ax1.set_ylabel('x(t)')
ax2=fig.add_subplot(2,1,2)
ax2.stem(tn,x,'r')
ax2.axis([0,nt*T,1.2*min(x),1.1*max(x)])
ax2.set_ylabel('x(n)')
plt.show()

 18、随机函数，这个用numpy就可以直接生成：

#生成随机函数：
t=np.linspace(0,49,dtype=int)
N=len(t)
x=np.random.random(len(t))
fig2=plt.figure()
ax1=fig2.add_subplot(2,1,1)
ax1.plot(t,x)
ax1.set_xlabel('n')
ax1.set_ylabel('x(n)')
ax2=fig2.add_subplot(2,1,2)
ax2.stem(t,x)
ax2.set_xlabel('n')
ax2.set_ylabel('x(n)')
plt.show()

 19、三角波，用signal的sawtooth：

#生成三角波
fs=10000
t=np.linspace(0,1,fs,dtype=float)
x1=signal.sawtooth(1*np.pi*40*t,0)
x2=signal.sawtooth(1*np.pi*40*t,1)
fig3=plt.figure()
ax1=fig3.add_subplot(2,1,1)
ax1.plot(t,x1)
ax1.set_xlabel('n')
ax1.set_ylabel('x(n)')
ax1.axis([0,0.25,-1,1])
ax2=fig3.add_subplot(2,1,2)
ax2.plot(t,x2)
ax2.set_xlabel('n')
ax2.set_ylabel('x(n)')
ax2.axis([0,0.25,-1,1])
plt.show()

 20、sinc曲线:

t=np.linspace(-3*np.pi,4*np.pi,400,dtype=float)
plt.plot(t,np.sinc(t))
plt.show()

 但是，有一点我到现在还没想通，同样的sinc，用公式推到的和内置的，出来的效果就是不一样：

t=np.linspace(-10,10,10000,dtype=float)
x=np.random.randint(len(t))
y=np.sinc(t)
y1=np.sin(t)
# y2=np.divide(y1,t)
y2=y1/t
fig3=plt.figure()
ax1=fig3.add_subplot(3,1,1)
ax1.plot(t,y)
ax2=fig3.add_subplot(3,1,2)
ax2.plot(t,y1)
ax3=fig3.add_subplot(3,1,3)
ax3.plot(t,y2)
ax3.plot(t,y,'r*')
plt.show()
#好奇怪啊，为什么，已知的是sinc(t)=sin(t)/t，但是我用这种方法做出来的，两个信号的宽度不一样
from mayavi import mlab
from mpl_toolkits.mplot3d import Axes3D
fig=plt.figure()
z=np.linspace(-10,2500,10000,dtype=float).reshape(100,100)
ax=fig.add_subplot(1,1,1, projection='3d')
t2=t.reshape(100,100)
x2=np.sinc(t2)
ax.plot_surface(t2,x2,z,rstride=4,cstride=3,color='r',alpha=0.9)

 三维的我会展示另一幅，以为上面这个代码中的三维太难看了，

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 生成数据
x = y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X**2 + Y**2)
# Z = np.sin(R) / R
z=np.sinc(R)

# 绘制图像
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, z, cmap='coolwarm')
plt.show()

 就这样吧

21、生成非周期三角波

这个利用自己写的函数生成一个就好：

#非周期三角波信号：
t=np.linspace(-2,2,4000)
def triangle_wave(x,c,hc): #幅度为hc，宽度为c,斜度为hc/2c的三角波
    if x>=c/2:
        r=0.0
    elif x<=-c/2:
        r=0.0
    elif((x>-c/2)and(x<0)):
        r=2*x/c*hc+hc
    else:
        r=-2*x/c*hc+hc
    return r
x=np.array([triangle_wave(i,0.5,0.5) for i in t ])
plt.plot(t,x)
plt.show()

 22、高斯脉冲的实现：

#高斯正弦脉冲：signal.gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False,retenv=False):
tc=signal.gausspulse('cutoff',60e3,0.6,tpr=-40)
tG=np.linspace(-tc,tc,100000)
y=signal.gausspulse(tG,60e3,0.6)
plt.plot(tG,y)
plt.show()

 23、脉冲序列发生器：

这个就是上面的square函数和sawtooth：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
n=np.linspace(0,10,dtype=float)
h=signal.square(1*np.pi*40*n)
z=signal.sawtooth(1*np.pi*40*n)
fig=plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.plot(n,h)
ax1.set_title('square')
ax2=fig.add_subplot(2,1,2)
ax2.plot(n,z)
ax2.set_title('sawtooth')
plt.show()

24、产生非周期方波：

#产生非周期方波:matlab 中用的是rectpuls,python似乎还没有这个函数，但是可以自己实现，这是我复制的别人滴
def rect_wave(x,c,c0):     #起点为c0，宽度为c的矩形波
     if x>=(c+c0):
          r=0.0
     elif x<c0:
          r=0.0
     else:
          r=1
     return r

t=np.linspace(-3,3,6000)
y1=np.array([rect_wave(i,1.0,-2.0) for i in t])
plt.plot(t,y1)
plt.show()

25、连续时间信号的时域分析

（注意，到这里就很注重系统的概念了）

就是求解齐次方程和非齐次方程求解零输入响应和零状态响应：

（1）零状态响应：

#连续时间系统数值求解：matlab有提供函数lsim，python中应该有：scipy包里面有lsim函数：def lsim(system, U, T, X0=None, interp=True):和matlab中的几乎一模一样
from scipy import signal as signal
ts=0
te=5
dt=0.01
# 计算系统的零状态响应,系统是：y''+2y'+100y=10cos(2*pi*t)
# 这里，第一项是右端系数，第二项是左端系数，实际上，lti官方文档给的示例是4个二维矩阵，但我没有明白他们是要干什么。
sys=signal.lti([1],[1,2,200])
t=np.linspace(ts,te,500)
f=10*np.cos(2*np.pi*t)
T,yout,xout=signal.lsim2(sys,f,t)
plt.plot(t,yout)
plt.show()
# 和matlab书上展示的一模一样,展示了该系统的零状态响应

 （2）冲激响应和阶跃响应：

# 连续时间系统系统冲击响应和节约响应：这个matlab有impuls和step，python也当然是有滴，就在singal里面:
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
t=np.linspace(0,4,2000)
system=([1,32],[1,4,64])
t,h=signal.impulse(system,T=t,N=2000)
t1,g=signal.step(system,T=t,N=2000)
fig1=plt.figure()
ax1=fig1.add_subplot(2,1,1)
ax1.plot(t,h)
ax1.set_xlabel("t")
ax1.set_ylabel('h(t)')
ax1.set_title('impulse')
ax1.grid(visible=True)
ax1.axis([0,4,-1.5,3])
ax2=fig1.add_subplot(2,1,2)
ax2.plot(t1,g)
ax2.set_xlabel("t")
ax2.set_ylabel('g(t)')
ax2.set_title('step')
ax2.grid(visible=True)
ax2.axis([0,4,0,1])
plt.show()

 （3）各种信号的响应：

#计算一个特定系统：
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

b=[-0.46,-0.25,-0.12,-0.06]
a=[1,0.64,0.94,0.51,0.01]
system=(b,a)
t1=np.linspace(0,10,10000)
t2=np.linspace(-5,5,10000)
f1=np.zeros(len(t1))
def rect_wave(x,c,c0):     #起点为c0，宽度为c的矩形波
     if x>=(c+c0):
          r=0.0
     elif x<c0:
          r=0.0
     else:
          r=1
     return r
def heaviside(x):
     if x==0:
          r=0.5
     elif x>0:
          r=1
     elif x<0:
          r=0
     return r
def unit(t):
    r=np.where(t>0.0,1.0,0.0)
    return r
f1=np.array([heaviside(i)for i in t1])-np.array([heaviside(i) for i in t1])
# f2=np.array([rect_wave(i,1,-1) for i in t2])
f2=np.ones(len(t1))
f2[10:11]=0
f3=t1
f4=np.sin(t1)
T1,yout1,xout1=signal.lsim2(system,f1,t1)
T2,yout2,xout2=signal.lsim2(system,f2,t1)
T3,yout3,xout3=signal.lsim2(system,f3,t1)
T4,yout4,xout4=signal.lsim2(system,f4,t1)
#f1、f2这两个有些问题
fig2=plt.figure()
ax1=fig2.add_subplot(2,2,1)
ax1.plot(t1,yout1)
ax1.set_title('chongji')
ax2=fig2.add_subplot(2,2,2)
ax2.plot(t1,yout2)
ax2.set_title('jump')
ax3=fig2.add_subplot(2,2,3)
ax3.plot(t1,yout3)
ax3.set_title('xiepo')
ax4=fig2.add_subplot(2,2,4)
ax4.plot(t1,yout4)
ax4.set_title('sin')

 [0,0]这幅图有些问题，其他的都没啥问题。

（4）连续时间信号的卷积：

#连续时间信号卷积求解：
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
dt=0.01
t=np.linspace(-1,2.5,350)
t1=t-2
# matlab的heaviside函数是一个阶跃函数，当输入为0时返回0.5，当输入大于0时返回1，当输入小于0时返回0。它的定义是:
# heaviside(x) = 0.5, x = 0
# heaviside(x) = 1, x > 0
# heaviside(x) = 0, x < 0
def heaviside(x):
    if x==0:
          r=0.5
    elif x>0:
          r=1
    elif x<0:
          r=0
    return r
def unit(t):
    r=np.where(t>0.0,1.0,0.0)
    return r
f1=np.array([heaviside(i)for i in t])-np.array([heaviside(i) for i in t1])*0.5
f2=2*np.exp(-3*t)*np.array([heaviside(i)for i in t])
f=np.convolve(f1,f2)*dt
n=len(f)
tt=np.linspace(0,n,n)*dt-2
fig3=plt.figure()
ax1=fig3.add_subplot(3,2,1)
ax1.plot(t,f1)
ax1.axis([-1,2.5,0,1.2])
ax1.set_title('f1(t)')
ax1.grid(visible=True)
ax2=fig3.add_subplot(3,2,2)
ax2.plot(t,f2)
ax2.axis([-1,2.5,0,2])
ax2.set_title('f2(t)')
ax2.grid(visible=True)
ax3=fig3.add_subplot(3,1,2)
ax3.plot(tt,f)
ax3.axis([-2,5,0,1])
ax3.set_title('convolve')
ax3.grid(visible=True)
f4=signal.convolve(f1,f2)*dt
ax4=fig3.add_subplot(3,1,3)
ax4.plot(tt,f4)
ax4.axis([-2,5,0,1])
ax4.set_title('convolve')
ax4.grid(visible=True)
plt.show()

 这里第三张用的是numpy的卷积函数，第四张用的是signal的卷积函数，我就是想看看，两者有没有上面出入。

26、离散时间系统：

实现离散时间系统的实现连续时间系统是很简单的，用numpy就行，说白了，连续就是点很多，很密，就和微分一样，我分的越细，它逼近的就越好，离散也一样，你就认为在和连续时间信号两者范围相同的情况下，区区有限个点就行：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
#离散
t=np.linspace(-10,10,41,dtype=float)
#连续
t1=np.linspace(-10,10,10000000,dtype=float)
x=signal.square(np.pi*2*t)
x1=signal.square(np.pi*2*t1)
fig=plt.figure()
ax1=fig.add_subplot(2,1,1)
ax1.stem(t,x)
ax1.set_title('discrete')
ax2=fig.add_subplot(2,1,2)
ax2.plot(t1,x1)
ax2.set_title('coiled')
plt.show()

 可以看到，上下两图的区别。

27、离散时间系统的冲激和阶跃响应：

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
a=[1,-0.35,1.5]
b=[1,1]
system=(b,a)
t=np.linspace(0,21,21,dtype=int)
x=np.power(0.5,t)
T,yout,xout=signal.lsim2(system,T=t)
print(T.size,yout.size,xout.size)
fig1=plt.figure()
ax1=fig1.add_subplot(1,2,1)
ax1.stem(t,x)
ax1.set_title('imput Sequence')
ax1.grid(visible=True)
ax2=fig1.add_subplot(1,2,2)
ax2.stem(t,yout)
ax2.set_title('output Sequence')
ax2.grid(visible=True)
plt.show()

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
a=[1,6,4]
b=[1,3]
k=np.linspace(0,10,11,dtype=int)
system=signal.dlti(b,a)
T,yout=signal.dstep(system,t=k)
plt.stem(T,np.squeeze(yout))

 28、卷积和运算，这都不用说啥了，前面都说过了：

#离散时间信号的卷积和运算
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
def uDt(n):
    if n>=0:
        y=1
    else:
        y=0
    return y
nx=np.linspace(-1,5,7,dtype=int,endpoint=True)
nh=np.linspace(-2,10,13,dtype=int,endpoint=True)
nx2=nx-4
nh2=nh-9
x=np.array([uDt(i)for i in nx])-np.array([uDt(i)for i in nx2])
h=np.power(0.9,nh)
h1=np.array([uDt(i)for i in nh])-np.array([uDt(i)for i in nh2])
y=np.convolve(h,h1)
ny1=nx[0]+nh[0]
ny=ny1+np.linspace(0,(len(y)),len(y),dtype=int)
fig=plt.figure()
ax1=fig.add_subplot(3,1,1)
ax1.stem(nx,x)
ax1.grid(visible=True)
ax1.set_title('x(n)')
ax1.axis([-4,16,0,1.5])
ax2=fig.add_subplot(3,1,2)
ax2.stem(nh,h)
ax2.grid(visible=True)
ax2.set_title('h(n))')
ax2.axis([-4,16,0,1.5])
ax3=fig.add_subplot(3,1,3)
ax3.stem(ny,y)
ax3.set_title('y(n)')
ax3.grid(visible=True)
ax3.axis([-4,16,0,9])
plt.show()

#已知序列卷积求和：
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from scipy.stats import pearsonr
x=np.array([1,3,5,7,9,11,13,15,17,19])
y=np.array([1,1,1,1,2,2,2,2,2,2])
z=np.convolve(x,y)
xlength=np.linspace(0,len(x),len(x),dtype=int)
ylength=np.linspace(0,len(y),len(y),dtype=int)
zlength=np.linspace(0,len(z),len(z),dtype=int)
figure=plt.figure()
ax1=figure.add_subplot(3,1,1)
ax1.stem(xlength,x)
ax1.set_title('x(n)')
ax1.grid(visible=True)
ax1.axis([0,len(x),0,20])
ax2=figure.add_subplot(3,1,2)
ax2.stem(ylength,y)
ax2.set_title('y(n)')
ax2.grid(visible=True)
ax2.axis([0,len(y),0,2.2])
ax3=figure.add_subplot(3,1,3)
ax3.stem(zlength,z)
ax3.set_title('z(n)=x(n)*y(n)')
ax3.grid(visible=True)
ax3.axis([0,20,0,max(z)+10])
plt.show()
print(np.corrcoef(x,y))
#这个pearsonr就不用看来，他输出的结果是xy这两个序列的相关系数矩阵，这个在统计学里面会用得到，但现在似乎没有任何用，是我刚开始搞错了
corr_coef, p_value = pearsonr(x, y)
print(corr_coef,p_value)

 29、相关序列（自相关、互相关）：

这些前面都讲过了，没啥意思的，这里就不再说了

其实，这两天主要就在忙这个，因为，怎么说呢，我是可以看懂matlab的，但是你要让我去用matlab写代码，我是一百万个不情愿，可能是因为1、我的水平还很低级，2、本科的时候学的是面向对象的C++和Python，对面向对象的编程方式比较熟悉，所以我就选择了使用python来学习数字信号处理，Python很强大，而且都是开源的，对于我来说，Python用起来比Matlab顺手的多，当然，这也是个人原因，办公室的几个师兄师姐就觉得matlab比c++和python简单，所以他们Matlab用的多，但是，说白了，编程语言就是个简单工具，最重要的还是算法和思想。

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：/a/24674.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

开源云原生数仓引擎ByConity 存储计算分离架构和优势

开源云原生数仓引擎ByConity 存储计算分离架构和优势

供稿 | ByConity技术团队 出品 | CSDN 云计算 ByConity是一款字节跳动开源的云原生数仓引擎。它的一个重要优势是采用存储计算分离的架构&#xff0c;实现了读写分离和弹性扩缩容。这种架构确保读操作和写操作不会相互影响&#xff0c;使得计算资源和存储资源解耦&#xff0c;…
阅读更多...
基于html+css的图展示102

基于html+css的图展示102

准备项目 项目开发工具 Visual Studio Code 1.44.2 版本: 1.44.2 提交: ff915844119ce9485abfe8aa9076ec76b5300ddd 日期: 2020-04-16T16:36:23.138Z Electron: 7.1.11 Chrome: 78.0.3904.130 Node.js: 12.8.1 V8: 7.8.279.23-electron.0 OS: Windows_NT x64 10.0.19044 项目…
阅读更多...
BLE协议栈结构

BLE协议栈结构

// 开坑BLE协议栈 0 镇楼图 接下来会自下往上粗略分析各个层级&#xff0c;后续会有对各层的细致解读 1 CONTROLLER 1.1 PHY BLE使用ISM频段&#xff08;频率范围是2.400-2.4835 GHz&#xff09;。将整个频带分为40份&#xff0c;每份的带宽为2MHz&#xff0c;称作RF Chann…
阅读更多...
如何使用Python自动化测试工具Selenium进行网页自动化？

如何使用Python自动化测试工具Selenium进行网页自动化？

引言 Selenium是一个流行的Web自动化测试框架&#xff0c;它支持多种编程语言和浏览器&#xff0c;并提供了丰富的API和工具来模拟用户在浏览器中的行为。Selenium可以通过代码驱动浏览器自动化测试流程&#xff0c;包括页面导航、元素查找、数据填充、点击操作等。 与PyAuto…
阅读更多...
抖音账号矩阵系统源码开发之——视频发布功能开发

抖音账号矩阵系统源码开发之——视频发布功能开发

视频发布权限在账号矩阵系统研发之初&#xff0c;都是一个备受争议的功能&#xff0c;最早之前我们使用的视频发布权限名字是Video.creat, video.delete权限&#xff0c;但是该权限于2022年10月份做了权限的收回&#xff0c;后又在上架了一个能力叫发布内容至抖音&#xff1a;…
阅读更多...
PostGIS的10个最佳实践

PostGIS的10个最佳实践

PostGIS 是一个功能强大的开源空间数据库&#xff0c;可用于存储、查询和分析地理空间数据。 对于需要存储和分析大量地理空间数据的组织来说&#xff0c;这是一个流行的选择。 但是&#xff0c;正确使用 PostGIS 以充分利用它很重要。 在本文中&#xff0c;我们将讨论 10 个 …
阅读更多...
【2023年4月美赛加赛】Z题：The future of Olympics 25页完整论文

【2023年4月美赛加赛】Z题：The future of Olympics 25页完整论文

【2023年4月美赛加赛】Z题&#xff1a;The future of Olympics 25页完整论文 1 题目 背景 国际奥委会(IOC)正面临着夏季奥运会和冬季奥运会申办数量的减少**[1]**。在过去&#xff0c;举办奥运会的竞争非常激烈&#xff0c;声望也很高。然而&#xff0c;最近&#xff0c;主办…
阅读更多...
自定义注解和@Target、@Retention注解的使用

自定义注解和@Target、@Retention注解的使用

说明&#xff1a;注解可以理解为另一种形式的配置&#xff0c;可用于在类上、方法上等&#xff0c;标志是“”&#xff0c;如重写方法上的“Override”就是一种注解。这里我通过一个实例&#xff0c;来介绍自定义注解和java元注解&#xff08;Target、Retention&#xff09;的使…
阅读更多...
一分钟了解乐观锁、悲观锁、共享锁、排它锁、行锁、表锁以及使用场景

一分钟了解乐观锁、悲观锁、共享锁、排它锁、行锁、表锁以及使用场景

大家好&#xff0c;我是冰点&#xff0c;今天给大家带来&#xff0c;关于MySQL中的锁的使用。 我首先提个问题&#xff0c;大家知道什么是 乐观锁、悲观锁、共享锁&#xff0c;、排它锁、行锁、表锁&#xff0c;以及每种锁的使用场景吗&#xff1f; !! 背景&#xff1a;最近在各…
阅读更多...
（学习日记）2023.04.23

（学习日记）2023.04.23

写在前面&#xff1a; 由于时间的不足与学习的碎片化&#xff0c;写博客变得有些奢侈。 但是对于记录学习&#xff08;忘了以后能快速复习&#xff09;的渴望一天天变得强烈。 既然如此 不如以天为单位&#xff0c;以时间为顺序&#xff0c;仅仅将博客当做一个知识学习的目录&a…
阅读更多...
关于安卓以及微软用户chatgpt上一篇文章如今第五点无法正常进入更新解决方法以及附加本地部署

关于安卓以及微软用户chatgpt上一篇文章如今第五点无法正常进入更新解决方法以及附加本地部署

目录 一、问题出现&#xff1a; 1、问题&#xff1a; 原因&#xff1a; 二、解决办法&#xff08;本地部署chatgpt&#xff09; 1、解决&#xff08;国内网络使用真的chatgpt并非镜像&#xff09;一次部署终生使用 第一步&#xff1a; ​编辑第二步&#xff1a; 三、实现结…
阅读更多...
让你不再好奇怎么给小说配音

让你不再好奇怎么给小说配音

你是否曾经想象过&#xff0c;当你在读小说时&#xff0c;你可以听到人物的声音&#xff0c;感受到情感和气氛的变化&#xff1f;有声书的出现已经让这一切成为可能。然而&#xff0c;如何为小说创造生动的配音效果却是一个需要仔细考虑的问题。如果你还不知道怎么给小说配音的…
阅读更多...
智能计价器-第14届蓝桥杯省赛Scratch中级组真题第5题

智能计价器-第14届蓝桥杯省赛Scratch中级组真题第5题

[导读]&#xff1a;超平老师的《Scratch蓝桥杯真题解析100讲》已经全部完成&#xff0c;后续会不定期解读蓝桥杯真题&#xff0c;这是Scratch蓝桥杯真题解析第140讲。 智能计价器&#xff0c;本题是2023年5月7日举行的第14届蓝桥杯省赛Scratch图形化编程中级组真题第5题&#…
阅读更多...
信号处理与分析-卷积的性质与推导

信号处理与分析-卷积的性质与推导

目录 一、引言 二、信号分析中的卷积 1. 什么是卷积 2. 卷积的性质 3. 卷积的应用 三、离散卷积 1. 离散卷积的定义 2. 离散卷积的计算 3. 离散卷积的性质 四、连续卷积 五、卷积的实际应用 六、总结 一、引言 在信号处理中&#xff0c;卷积是一种非常重要的数学运…
阅读更多...
2023年05月小结

2023年05月小结

5月一晃即将过去&#xff0c;先看看5月份都写了哪些文章&#xff1a; 往期推荐 国内手机号如何注册Google账号 Python中的布尔类型以及布尔值介绍 Python中深拷贝与浅拷贝介绍 Python中的集合介绍 Python中的字典学习笔记 Python学习历程-元组的基本操作回顾 记录我的Python学习…
阅读更多...
STM32F103C8T6+2.4寸SPI TFT触摸屏代码+标准库 项目开发

STM32F103C8T6+2.4寸SPI TFT触摸屏代码+标准库 项目开发

目录 模块清单&#xff1a; 模块介绍&#xff1a; 1&#xff1a;STM32F103C8T6 2&#xff1a;2.4寸SPI TFT触摸屏 项目结果展示 2.4寸 TFT SPI显示触摸屏 2.4寸 SPI TFT 显示触摸屏代码下载链接&#xff1a; (1条消息) 2.4寸SPITFT显示触摸屏资源-CSDN文库 模块清单&#x…
阅读更多...
低资源方面级情感分析研究综述

低资源方面级情感分析研究综述

文章目录 前言1. 引言2. 问题定义、数据集和评价指标2.1 问题定义2.2 任务定义2.3 常用数据集 3. 方面级情感分析的方法3.1 **方面词抽取**3.1.1 基于无监督学习的方法3.1.1.1 基于规则的方面词抽取3.1.1.2 基于统计的方面词抽取 3.1.2 基于有监督浅层模型的方法3.1.3 基于有监…
阅读更多...
Golang每日一练(leetDay0075) 打家劫舍II、最短回文串

Golang每日一练(leetDay0075) 打家劫舍II、最短回文串

目录 213. 打家劫舍 II House Robber ii &#x1f31f;&#x1f31f; 214. 最短回文串 Shortest Palindrome &#x1f31f;&#x1f31f;&#x1f31f; &#x1f31f; 每日一练刷题专栏 &#x1f31f; Rust每日一练 专栏 Golang每日一练 专栏 Python每日一练 专栏 C/C…
阅读更多...
mysql加强小结 203446

mysql加强小结 203446

数据库三范式: 什么是范式 规则:想要设计一个好的关系,必须要满足一定的约束条件,有几个等级,一级比一级高 ​ 解决什么问题:让数据库设计更加简洁,结构更加清晰,否则容易造成数据冗余 数据库有哪些范式? ​ 数据库有七大范式,常用的只有三个范式 **第一范式:**业务上属…
阅读更多...
【1++的C++初阶】之模板

【1++的C++初阶】之模板

&#x1f44d;作者主页&#xff1a;进击的1 &#x1f929; 专栏链接&#xff1a;【1的C初阶】 文章目录 一&#xff0c;浅谈模板二&#xff0c;函数模板三&#xff0c;类模板 一&#xff0c;浅谈模板 在前面的文章【【1的C初阶】之C入门篇1】中我们对函数重载有了一定的认识&a…
阅读更多...
最新文章

Copyright @ 2022~2023