• 数字信号处理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、离散时间信号序列的表示：

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

     2、采样定理：

    1. #用10hz的采样频率采样
    2. import matplotlib.pyplot as plt
    3. import numpy as np
    4. from scipy import signal
    5. dt=0.01
    6. n=np.linspace(0,89,90,dtype=int)
    7. t=n*dt
    8. f=10
    9. x=np.sin(3*np.pi*f*t+0.5)#原始信号
    10. dt=0.1
    11. n=np.linspace(0,9,10,dtype=int)
    12. t1=n*dt
    13. x1=np.sin(3*f*np.pi*t1+0.5)# 采样后的信号
    14. fig1=plt.figure()
    15. ax1=fig1.add_subplot(3,1,1)
    16. ax1.plot(t,x)
    17. ax1.set_title('origional signal')
    18. ax2=fig1.add_subplot(3,1,2)
    19. ax2.plot(t,x)
    20. ax2.plot(t1,x1,'*')
    21. ax2.set_title('Sampling process')
    22. ax3=fig1.add_subplot(3,1,3)
    23. ax3.plot(t1,x1)
    24. ax3.set_title('Sampling signal')
    25. plt.show()

     3、简单离散信号序列：

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy import signal
    4. n=50
    5. x=np.zeros(n)
    6. x[1]=1
    7. xn=np.linspace(0,n-1,n,dtype=int)
    8. fig1=plt.figure()
    9. ax1=fig1.add_subplot(2,1,1)
    10. ax1.stem(xn,x)
    11. x[1]=0
    12. x[10]=1
    13. ax2=fig1.add_subplot(2,1,2)
    14. ax2.stem(xn,x)
    15. plt.show()

     4、单位阶跃函数：

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

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy import signal
    4. def u(n):
    5.     if n>=0:
    6.         r=1
    7.     else:
    8.         r=0
    9.     return r
    10. x=np.linspace(-10,10,21,dtype=int)
    11. y=np.array([u(i)for i in x])
    12. plt.stem(x,y)
    13. plt.show()

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

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

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

    5、正弦信号序列：

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

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

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

    6、实指数序列：

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

    x(n)=a^nu(n)

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy import signal
    4. import sympy as sy
    5. n=np.linspace(-10,10,21,dtype=int)
    6. y1=np.power(1.6,n)
    7. y2=np.power(-1.6,n)
    8. y3=np.power(0.9,n)
    9. y4=np.power(-0.9,n)
    10. f=plt.figure()
    11. ax1=f.add_subplot(2,2,1)
    12. ax1.stem(n,y1)
    13. ax2=f.add_subplot(2,2,2)
    14. ax2.stem(n,y2)
    15. ax3=f.add_subplot(2,2,3)
    16. ax3.stem(n,y3)
    17. ax4=f.add_subplot(2,2,4)
    18. ax4.stem(n,y4)
    19. ax1.set_title('1.6^n')
    20. ax2.set_title('(-1.6)^n')
    21. ax3.set_title('0.9^n')
    22. ax4.set_title('(-0.9)^n')
    23. plt.show()

     7、复指数序列：

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

    8、序列值累加和乘积：

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

    9、序列反转、移位：

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

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

    1. #序列位移：
    2. import matplotlib.pyplot as plt
    3. import numpy as np
    4. nx=np.linspace(-2,5,8,dtype=int)
    5. x=np.array([9,8,7,6,5,5,5,5])
    6. y=x
    7. ny1=nx+3
    8. ny2=nx-2
    9. fig=plt.figure()
    10. ax1=fig.add_subplot(2,1,1)
    11. ax1.stem(nx,x,'.')
    12. ax1.axis([-5,9,-1,10])
    13. ax1.grid(visible=True)
    14. ax1.set_xlabel('n')
    15. ax1.set_ylabel('x(n)')
    16. ax1.set_title('origional sequence')
    17. ax2=fig.add_subplot(2,2,3)
    18. ax2.stem(ny1,y,'.')
    19. ax2.axis([-5,9,-1,10])
    20. ax2.grid(visible=True)
    21. ax2.set_xlabel('n')
    22. ax2.set_ylabel('y1(n)')
    23. ax2.set_title('y1=(n+3)')
    24. ax3=fig.add_subplot(2,2,4)
    25. ax3.stem(ny2,y,'.')
    26. ax3.axis([-5,9,-1,10])
    27. ax3.grid(visible=True)
    28. ax3.set_xlabel('n')
    29. ax3.set_ylabel('y2(n)')
    30. ax3.set_title('y2=(n-2)')
    31. plt.show()

     10、信号的尺度变换：

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

    1. t = np.linspace(-4, 4, 8000, dtype=float)
    2. T = 2
    3. f = np.zeros((8000),dtype=float)
    4. t1=2*t
    5. f1=np.zeros((8000),dtype=float)
    6. for i in range(len(t)):
    7.     if ((-1 <= t[i]) & (t[i] <= 1)).any():
    8.         f[i] = 1
    9. for i in range(len(t1)):
    10.     if ((-0.5 <= t1[i]) & (t1[i] <= 0.5)).any():
    11.         f1[i] = 1
    12. fig=plt.figure()
    13. ax1=fig.add_subplot(2,1,1)
    14. ax1.plot(t,f)
    15. ax1.axis([-4,4,-0.5,1.5])
    16. ax2=fig.add_subplot(2,1,2)
    17. ax2.plot(t1,f1)
    18. ax2.axis([-4,4,-0.5,1.5])
    19. plt.show()

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

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

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

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

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

    还是上面的

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

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

    13、信号的微积分：

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

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. import sympy as sp
    4. from scipy import signal
    5. import sympy.plotting as syp
    6. def heaviside(x):
    7.     if x==0:
    8.           r=0.5
    9.     elif x>0:
    10.           r=1
    11.     elif x<0:
    12.           r=0
    13.     return r
    14. #本来要用的，但我们没有用heaviside函数
    15. t=sp.symbols('t')
    16. f=sp.Function('f')(t)
    17. f=t*t+2*t-1
    18. f1=f.diff(t)#一阶导
    19. f2=f.diff(t,t)#二阶导
    20. f3=f.diff(t,t,t)#三阶导
    21. syp.plot(f,f1,f2,f3,(t,-1,2))

     14、积分：

    这里用的是：

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

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

    15、卷积运算和相关计算

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

    1. #离散时间信号的卷积和运算
    2. import matplotlib.pyplot as plt
    3. import numpy as np
    4. from scipy import signal
    5. def uDt(n):
    6.     if n>=0:
    7.         y=1
    8.     else:
    9.         y=0
    10.     return y
    11. nx=np.linspace(-1,5,7,dtype=int,endpoint=True)
    12. nh=np.linspace(-2,10,13,dtype=int,endpoint=True)
    13. nx2=nx-4
    14. nh2=nh-9
    15. x=np.array([uDt(i)for i in nx])-np.array([uDt(i)for i in nx2])
    16. h=np.power(0.9,nh)
    17. h1=np.array([uDt(i)for i in nh])-np.array([uDt(i)for i in nh2])
    18. y=np.convolve(h,h1)
    19. ny1=nx[0]+nh[0]
    20. ny=ny1+np.linspace(0,(len(y)),len(y),dtype=int)
    21. fig=plt.figure()
    22. ax1=fig.add_subplot(3,1,1)
    23. ax1.stem(nx,x)
    24. ax1.grid(visible=True)
    25. ax1.set_title('x(n)')
    26. ax1.axis([-4,16,0,1.5])
    27. ax2=fig.add_subplot(3,1,2)
    28. ax2.stem(nh,h)
    29. ax2.grid(visible=True)
    30. ax2.set_title('h(n))')
    31. ax2.axis([-4,16,0,1.5])
    32. ax3=fig.add_subplot(3,1,3)
    33. ax3.stem(ny,y)
    34. ax3.set_title('y(n)')
    35. ax3.grid(visible=True)
    36. ax3.axis([-4,16,0,9])
    37. plt.show()

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

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

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

    16、产生信号波形：

    1. from scipy import signal as signal
    2. t=np.linspace(0,1,200,dtype=float)
    3.  
    4. # scipy.signal.chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True)
    5.  
    6. h=signal.chirp(t,0,1,120,method='linear',phi=np.pi/3)#linear线性\quadratic二次扫描、logarithmic对数扫描（这时候f0、f1均不能为零）
    7. plt.plot(t,h)
    8. plt.grid(visible=True)
    9. plt.show()

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

    来看一下：

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

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy import signal
    4. f=6000
    5. nt=3
    6. N=15
    7. T=1/f
    8. dt=T/N
    9. n=np.linspace(0,50,51,dtype=int)
    10. tn=n*dt
    11. x=signal.square(2*f*np.pi*tn,duty=0.25)+1
    12. fig=plt.figure()
    13. ax1=fig.add_subplot(2,1,1)
    14. ax1.step(tn,x,'k')
    15. ax1.axis([0,nt*T,1.2*min(x),1.1*max(x)])
    16. ax1.set_ylabel('x(t)')
    17. ax2=fig.add_subplot(2,1,2)
    18. ax2.stem(tn,x,'r')
    19. ax2.axis([0,nt*T,1.2*min(x),1.1*max(x)])
    20. ax2.set_ylabel('x(n)')
    21. plt.show()

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

    1. #生成随机函数：
    2. t=np.linspace(0,49,dtype=int)
    3. N=len(t)
    4. x=np.random.random(len(t))
    5. fig2=plt.figure()
    6. ax1=fig2.add_subplot(2,1,1)
    7. ax1.plot(t,x)
    8. ax1.set_xlabel('n')
    9. ax1.set_ylabel('x(n)')
    10. ax2=fig2.add_subplot(2,1,2)
    11. ax2.stem(t,x)
    12. ax2.set_xlabel('n')
    13. ax2.set_ylabel('x(n)')
    14. plt.show()

     19、三角波，用signal的sawtooth：

    1. #生成三角波
    2. fs=10000
    3. t=np.linspace(0,1,fs,dtype=float)
    4. x1=signal.sawtooth(1*np.pi*40*t,0)
    5. x2=signal.sawtooth(1*np.pi*40*t,1)
    6. fig3=plt.figure()
    7. ax1=fig3.add_subplot(2,1,1)
    8. ax1.plot(t,x1)
    9. ax1.set_xlabel('n')
    10. ax1.set_ylabel('x(n)')
    11. ax1.axis([0,0.25,-1,1])
    12. ax2=fig3.add_subplot(2,1,2)
    13. ax2.plot(t,x2)
    14. ax2.set_xlabel('n')
    15. ax2.set_ylabel('x(n)')
    16. ax2.axis([0,0.25,-1,1])
    17. plt.show()

     20、sinc曲线:

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

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

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

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

    1. import numpy as np
    2. import matplotlib.pyplot as plt
    3. from mpl_toolkits.mplot3d import Axes3D
    4.  
    5. # 生成数据
    6. x = y = np.linspace(-10, 10, 100)
    7. X, Y = np.meshgrid(x, y)
    8. R = np.sqrt(X**2 + Y**2)
    9. # Z = np.sin(R) / R
    10. z=np.sinc(R)
    11.  
    12. # 绘制图像
    13. fig = plt.figure()
    14. ax = fig.add_subplot(111, projection='3d')
    15. ax.plot_surface(X, Y, z, cmap='coolwarm')
    16. plt.show()

     就这样吧

    21、生成非周期三角波

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

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

     22、高斯脉冲的实现：

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

     23、脉冲序列发生器：

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

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy import signal
    4. n=np.linspace(0,10,dtype=float)
    5. h=signal.square(1*np.pi*40*n)
    6. z=signal.sawtooth(1*np.pi*40*n)
    7. fig=plt.figure()
    8. ax1=fig.add_subplot(2,1,1)
    9. ax1.plot(n,h)
    10. ax1.set_title('square')
    11. ax2=fig.add_subplot(2,1,2)
    12. ax2.plot(n,z)
    13. ax2.set_title('sawtooth')
    14. plt.show()

    24、产生非周期方波：

    1. #产生非周期方波:matlab 中用的是rectpuls,python似乎还没有这个函数，但是可以自己实现，这是我复制的别人滴
    2. def rect_wave(x,c,c0):     #起点为c0，宽度为c的矩形波
    3.      if x>=(c+c0):
    4.           r=0.0
    5.      elif x
    6.           r=0.0
    7.      else:
    8.           r=1
    9.      return r
    10.  
    11. t=np.linspace(-3,3,6000)
    12. y1=np.array([rect_wave(i,1.0,-2.0) for i in t])
    13. plt.plot(t,y1)
    14. plt.show()

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

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

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

    （1）零状态响应：

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

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

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

     （3）各种信号的响应：

    1. #计算一个特定系统：
    2. import matplotlib.pyplot as plt
    3. import numpy as np
    4. from scipy import signal
    5.  
    6. b=[-0.46,-0.25,-0.12,-0.06]
    7. a=[1,0.64,0.94,0.51,0.01]
    8. system=(b,a)
    9. t1=np.linspace(0,10,10000)
    10. t2=np.linspace(-5,5,10000)
    11. f1=np.zeros(len(t1))
    12. def rect_wave(x,c,c0):     #起点为c0，宽度为c的矩形波
    13.      if x>=(c+c0):
    14.           r=0.0
    15.      elif x
    16.           r=0.0
    17.      else:
    18.           r=1
    19.      return r
    20. def heaviside(x):
    21.      if x==0:
    22.           r=0.5
    23.      elif x>0:
    24.           r=1
    25.      elif x<0:
    26.           r=0
    27.      return r
    28. def unit(t):
    29.     r=np.where(t>0.0,1.0,0.0)
    30.     return r
    31. f1=np.array([heaviside(i)for i in t1])-np.array([heaviside(i) for i in t1])
    32. # f2=np.array([rect_wave(i,1,-1) for i in t2])
    33. f2=np.ones(len(t1))
    34. f2[10:11]=0
    35. f3=t1
    36. f4=np.sin(t1)
    37. T1,yout1,xout1=signal.lsim2(system,f1,t1)
    38. T2,yout2,xout2=signal.lsim2(system,f2,t1)
    39. T3,yout3,xout3=signal.lsim2(system,f3,t1)
    40. T4,yout4,xout4=signal.lsim2(system,f4,t1)
    41. #f1、f2这两个有些问题
    42. fig2=plt.figure()
    43. ax1=fig2.add_subplot(2,2,1)
    44. ax1.plot(t1,yout1)
    45. ax1.set_title('chongji')
    46. ax2=fig2.add_subplot(2,2,2)
    47. ax2.plot(t1,yout2)
    48. ax2.set_title('jump')
    49. ax3=fig2.add_subplot(2,2,3)
    50. ax3.plot(t1,yout3)
    51. ax3.set_title('xiepo')
    52. ax4=fig2.add_subplot(2,2,4)
    53. ax4.plot(t1,yout4)
    54. ax4.set_title('sin')

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

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

    1. #连续时间信号卷积求解：
    2. import matplotlib.pyplot as plt
    3. import numpy as np
    4. from scipy import signal
    5. dt=0.01
    6. t=np.linspace(-1,2.5,350)
    7. t1=t-2
    8. # matlab的heaviside函数是一个阶跃函数，当输入为0时返回0.5，当输入大于0时返回1，当输入小于0时返回0。它的定义是:
    9. # heaviside(x) = 0.5, x = 0
    10. # heaviside(x) = 1, x > 0
    11. # heaviside(x) = 0, x < 0
    12. def heaviside(x):
    13.     if x==0:
    14.           r=0.5
    15.     elif x>0:
    16.           r=1
    17.     elif x<0:
    18.           r=0
    19.     return r
    20. def unit(t):
    21.     r=np.where(t>0.0,1.0,0.0)
    22.     return r
    23. f1=np.array([heaviside(i)for i in t])-np.array([heaviside(i) for i in t1])*0.5
    24. f2=2*np.exp(-3*t)*np.array([heaviside(i)for i in t])
    25. f=np.convolve(f1,f2)*dt
    26. n=len(f)
    27. tt=np.linspace(0,n,n)*dt-2
    28. fig3=plt.figure()
    29. ax1=fig3.add_subplot(3,2,1)
    30. ax1.plot(t,f1)
    31. ax1.axis([-1,2.5,0,1.2])
    32. ax1.set_title('f1(t)')
    33. ax1.grid(visible=True)
    34. ax2=fig3.add_subplot(3,2,2)
    35. ax2.plot(t,f2)
    36. ax2.axis([-1,2.5,0,2])
    37. ax2.set_title('f2(t)')
    38. ax2.grid(visible=True)
    39. ax3=fig3.add_subplot(3,1,2)
    40. ax3.plot(tt,f)
    41. ax3.axis([-2,5,0,1])
    42. ax3.set_title('convolve')
    43. ax3.grid(visible=True)
    44. f4=signal.convolve(f1,f2)*dt
    45. ax4=fig3.add_subplot(3,1,3)
    46. ax4.plot(tt,f4)
    47. ax4.axis([-2,5,0,1])
    48. ax4.set_title('convolve')
    49. ax4.grid(visible=True)
    50. plt.show()

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

    26、离散时间系统：

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

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy import signal
    4. #离散
    5. t=np.linspace(-10,10,41,dtype=float)
    6. #连续
    7. t1=np.linspace(-10,10,10000000,dtype=float)
    8. x=signal.square(np.pi*2*t)
    9. x1=signal.square(np.pi*2*t1)
    10. fig=plt.figure()
    11. ax1=fig.add_subplot(2,1,1)
    12. ax1.stem(t,x)
    13. ax1.set_title('discrete')
    14. ax2=fig.add_subplot(2,1,2)
    15. ax2.plot(t1,x1)
    16. ax2.set_title('coiled')
    17. plt.show()

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

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

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy import signal
    4. a=[1,-0.35,1.5]
    5. b=[1,1]
    6. system=(b,a)
    7. t=np.linspace(0,21,21,dtype=int)
    8. x=np.power(0.5,t)
    9. T,yout,xout=signal.lsim2(system,T=t)
    10. print(T.size,yout.size,xout.size)
    11. fig1=plt.figure()
    12. ax1=fig1.add_subplot(1,2,1)
    13. ax1.stem(t,x)
    14. ax1.set_title('imput Sequence')
    15. ax1.grid(visible=True)
    16. ax2=fig1.add_subplot(1,2,2)
    17. ax2.stem(t,yout)
    18. ax2.set_title('output Sequence')
    19. ax2.grid(visible=True)
    20. plt.show()

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

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

    1. #离散时间信号的卷积和运算
    2. import matplotlib.pyplot as plt
    3. import numpy as np
    4. from scipy import signal
    5. def uDt(n):
    6.     if n>=0:
    7.         y=1
    8.     else:
    9.         y=0
    10.     return y
    11. nx=np.linspace(-1,5,7,dtype=int,endpoint=True)
    12. nh=np.linspace(-2,10,13,dtype=int,endpoint=True)
    13. nx2=nx-4
    14. nh2=nh-9
    15. x=np.array([uDt(i)for i in nx])-np.array([uDt(i)for i in nx2])
    16. h=np.power(0.9,nh)
    17. h1=np.array([uDt(i)for i in nh])-np.array([uDt(i)for i in nh2])
    18. y=np.convolve(h,h1)
    19. ny1=nx[0]+nh[0]
    20. ny=ny1+np.linspace(0,(len(y)),len(y),dtype=int)
    21. fig=plt.figure()
    22. ax1=fig.add_subplot(3,1,1)
    23. ax1.stem(nx,x)
    24. ax1.grid(visible=True)
    25. ax1.set_title('x(n)')
    26. ax1.axis([-4,16,0,1.5])
    27. ax2=fig.add_subplot(3,1,2)
    28. ax2.stem(nh,h)
    29. ax2.grid(visible=True)
    30. ax2.set_title('h(n))')
    31. ax2.axis([-4,16,0,1.5])
    32. ax3=fig.add_subplot(3,1,3)
    33. ax3.stem(ny,y)
    34. ax3.set_title('y(n)')
    35. ax3.grid(visible=True)
    36. ax3.axis([-4,16,0,9])
    37. plt.show()

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

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

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

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

  • 相关阅读:
    关于android：Retrofit-@Body参数不能与表单或多部分编码一起使用
    phpstorm配置php运行环境
    Linux（b站视频兄弟连）自学笔记第十三章——Linux系统管理
    TCP相关面试题
    Spring SSM整合步骤
    第十章 单调栈 part03 84. 柱状图中最大的矩形
    C++——模板进阶
    vue项目类微信聊天页面，输入法弹出，ios的标题会整体上移问题
    【场外衍生品系列】雪球结构定价研究
    Oracle 存储过程游标遍历查询信息，并打印输出异常信息
  • 原文地址：https://blog.csdn.net/faltas/article/details/130910386