Python 信号与系统

信号与系统 电路与系统函数

scipy.signal信号处理的库

https://blog.csdn.net/chehec2010/article/details/117109246
freqs（b，a [，worN，plot]）计算模拟滤波器的频率响应。
freqs_zpk（z，p，k [，worN]）计算模拟滤波器的频率响应。
freqz（b [，a，worN，整体，情节，fs，…]）计算数字滤波器的频率响应。
freqz_zpk（z，p，k [，worN，整个，fs]）计算ZPK形式的数字滤波器的频率响应
butter（N，Wn [，btype，模拟量，输出，fs]）巴特沃思数字和模拟滤波器设计。

lti（系统）连续时间线性时不变系统基类。
StateSpace（ system，** kwargs）状态空间形式的线性时不变系统。
TransferFunction（* system，** kwargs）传递函数形式的线性时不变系统类。
ZerosPolesGain（* system，** kwargs）零，极点，增益形式的线性时不变系统类别。
lsim（系统，U，T [，X0，插值]）模拟连续时间线性系统的输出。
lsim2（系统[，U，T，X0]）使用ODE求解器模拟连续时间线性系统的输出scipy.integrate.odeint。
impulse（系统[，X0，T，N]）连续时间系统的脉冲响应。
impulse2（系统[，X0，T，N]）输入连续时间线性系统的脉冲响应。
step（系统[，X0，T，N]）连续时间系统的阶跃响应。
step2（系统[，X0，T，N]）连续时间系统的阶跃响应。
freqresp（系统[，w，n]）计算连续时间系统的频率响应。
bode（系统[，w，n]）

tf2zpk（b，a）从线性滤波器的分子，分母表示中返回零极点增益（z，p，k）表示形式。
tf2sos（b，a [，配对]）从传递函数表示中返回二阶部分
tf2ss（数字，巢穴）将函数转移到状态空间表示形式。
zpk2tf（z，p，k）从零和极返回多项式传递函数表示
zpk2sos（z，p，k [，配对]）从零，极点和系统增益返回二阶部分
zpk2ss（z，p，k）零极点增益表示到状态空间表示
ss2tf（A，B，C，D [，输入]）状态空间传递函数。

from scipy import signal
num = [1, 3, 3]
den = [1, 2, 1]
signal.TransferFunction(num, den)
TransferFunctionContinuous(array([1., 3., 3.]),array([1., 2., 1.]),dt: None)
signal.TransferFunction(num, den, dt=0.1)
TransferFunctionDiscrete(array([1., 3., 3.]),array([1., 2., 1.]),dt: 0.1)

  from scipy.signal import freqs, iirfilter
 b, a = iirfilter(4, [1, 10], 1, 60, analog=True, ftype='cheby1')
 w, h = freqs(b, a, worN=np.logspace(-1, 2, 1000))
  import matplotlib.pyplot as plt
 plt.semilogx(w, 20 * np.log10(abs(h)))
 plt.xlabel('Frequency')
 plt.ylabel('Amplitude response [dB]')
 plt.grid()
 plt.show()

.

# -*- coding: utf-8 -*-
from scipy import signal
import numpy as np
import matplotlib.pyplot as pl

pl.figure(figsize=(5,5))

b, a = signal.butter(6, 1.0, analog=1)
z,p,k = signal.tf2zpk(b, a)
pl.plot(np.real(p), np.imag(p), '^', label=u"6阶巴特沃斯极点")

b, a = signal.butter(7, 1.0, analog=1)
z,p,k = signal.tf2zpk(b, a)
pl.plot(np.real(p), np.imag(p), 's', label=u"7阶巴特沃斯极点")

pl.axis("equal")
pl.legend(loc="center right")
pl.show()

 b, a = signal.butter(2, 1.0, analog=1)
 np.real(b)
array([ 1.])
np.real(a)
array([ 1.        ,  1.41421356,  1.        ])

def freqz(b, a=1, worN=None, whole=0, plot=None):
    b, a = map(atleast_1d, (b,a))
    if whole:
        lastpoint = 2*pi
    else:
        lastpoint = pi
    if worN is None:
        N = 512
        w = numpy.arange(0,lastpoint,lastpoint/N)
    elif isinstance(worN, types.IntType):
        N = worN
        w = numpy.arange(0,lastpoint,lastpoint/N)
    else:
        w = worN
    w = atleast_1d(w)
    zm1 = exp(-1j*w)
    h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
    if not plot is None:
        plot(w, h)
    return w, h

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
import scipy.signal as signal

def logfreqz(b, a, f0, f1, fs, N):
    """
    以对数频率坐标计算滤波器b,a的频率响应
    f0, f1: 计算频率响应的开始频率和结束频率
    fs: 取样频率
    """
    w0, w1 = np.log10(f0/fs*2*np.pi), np.log10(f1/fs*2*np.pi)
    # 不包括结束频率
    w = np.logspace(w0, w1, N, endpoint=False)
    zm1 = np.exp(-1j*w)
    h = np.polyval(b[::-1], zm1) / np.polyval(a[::-1], zm1)
    return w/2/np.pi*fs, h
    
for n in range(1, 6):    
    # 设计n阶的通频为0.1*4000 = 400Hz的高通滤波器
    b, a = signal.iirfilter(n, [0.1, 1])    
    f, h = logfreqz(b, a, 10.0, 4000.0, 8000.0, 400)
    gain = 20*np.log10(np.abs(h))
    pl.semilogx(f, gain, label="N=%s" % n)
    slope = (gain[100]-gain[10]) / (np.log2(f[100]) - np.log2(f[10]))
    print "N=%s, slope=%s dB" % (n, slope)
pl.ylim(-100, 20)
pl.xlabel(u"频率(Hz)")
pl.ylabel(u"增益(dB)")
pl.legend()
pl.show()

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
def stoz(s):
    """
    将s复平面映射到z复平面
    为了方便起见，假设取样周期T=1
    """
    return (2+s)/(2-s)
    
def make_vline(x):
    return x + 1j*np.linspace(-100.0,100.0,20000)
    
fig = pl.figure(figsize=(7,3))    
axs = pl.subplot(121)
axz = pl.subplot(122)
for x in np.arange(-3, 4, 1):
    s = make_vline(x)
    z = stoz(s)
    axs.plot(np.real(s), np.imag(s))
    axz.plot(np.real(z), np.imag(z))

axs.set_xlim(-4,4)
axz.axis("equal")
axz.set_ylim(-3,3)

pl.show()

# -*- coding: utf-8 -*-
# 本程序演示如何用多个正弦波合成三角波
import numpy as np
import pylab as pl

# 取FFT计算的结果freqs中的前n项进行合成，返回合成结果，计算loops个周期的波形
def fft_combine(freqs, n, loops=1):
    length = len(freqs) * loops
    data = np.zeros(length)
    index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi)
    for k, p in enumerate(freqs[:n]):
        if k != 0: p *= 2 # 除去直流成分之外，其余的系数都*2
        data += np.real(p) * np.cos(k*index) # 余弦成分的系数为实数部
        data -= np.imag(p) * np.sin(k*index) # 正弦成分的系数为负的虚数部
    return index, data    

# 产生size点取样的三角波，其周期为1
def triangle_wave(size):
    x = np.arange(0, 1, 1.0/size)
    y = np.where(x<0.5, x, 0)
    y = np.where(x>=0.5, 1-x, y)
    return x, y

fft_size = 256

# 计算三角波和其FFT
x, y = triangle_wave(fft_size)
fy = np.fft.fft(y) / fft_size

# 绘制三角波的FFT的前20项的振幅，由于不含下标为偶数的值均为0， 因此取
# log之后无穷小，无法绘图，用np.clip函数设置数组值的上下限，保证绘图正确
pl.figure()
pl.plot(np.clip(20*np.log10(np.abs(fy[:20])), -120, 120), "o")
pl.xlabel("frequency bin")
pl.ylabel("power(dB)")
pl.title("FFT result of triangle wave")

# 绘制原始的三角波和用正弦波逐级合成的结果，使用取样点为x轴坐标
pl.figure()
pl.plot(y, label="original triangle", linewidth=2)
for i in [0,1,3,5,7,9]:
    index, data = fft_combine(fy, i+1, 2)  # 计算两个周期的合成波形
    pl.plot(data, label = "N=%s" % i)
pl.legend()
pl.title("partial Fourier series of triangle wave")
pl.show()

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl

sampling_rate = 8000
fft_size = 512
t = np.arange(0, 1.0, 1.0/sampling_rate)
x = np.sin(2*np.pi*156.25*t)  + 2*np.sin(2*np.pi*234.375*t)
xs = x[:fft_size]
xf = np.fft.rfft(xs)/fft_size
freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)
xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))
pl.figure(figsize=(8,4))
pl.subplot(211)
pl.plot(t[:fft_size], xs)
pl.xlabel(u"时间(秒)")
pl.title(u"156.25Hz和234.375Hz的波形和频谱")
pl.subplot(212)
pl.plot(freqs, xfp)
pl.xlabel(u"频率(Hz)")
pl.subplots_adjust(hspace=0.4)
pl.show()

import numpy as np
import scipy.signal as signal
import pylab as pl

def average_fft(x, fft_size):
    n = len(x) // fft_size * fft_size
    tmp = x[:n].reshape(-1, fft_size)
    tmp *= signal.hann(fft_size, sym=0)
    xf = np.abs(np.fft.rfft(tmp)/fft_size)
    avgf = np.average(xf, axis=0)
    return 20*np.log10(avgf)

 from scipy import signal
 import matplotlib.pyplot as plt

 b, a = signal.butter(4, 100, 'low', analog=True)
 w, h = signal.freqs(b, a)                        # 由分子分母的系数求解频率响应，w为频率，h为对应的响应。
 plt.semilogx(w, 20 * np.log10(abs(h))) 	#绘制幅频响应，频率轴取对数，幅度轴转换成dB。
 plt.title('Butterworth filter frequency response')
 plt.xlabel('Frequency [radians / second]')
 plt.ylabel('Amplitude [dB]')
 plt.margins(0, 0.1)
 plt.grid(which='both', axis='both')
 plt.axvline(100, color='green') # cutoff frequency
 plt.show()

表达式

$$
G(s)={4\over s^3+2s^2+10s+8}
$$

$$G(s)=\frac{4}{s^3+2s^2+10s+8}$$

符号算法

表示正弦信号

import numpy as np
from sympy import plot,sin,Symbol
t=Symbol('t')   #定义符号变量t
y=sin(np.pi/4*t)
plot(y)

代替control库的符号运算
在sympy中是否有一个函数可以将多项式转换为系数和参数的列表？

from sympy import * 
x, xa, xb, a0, a1, a2, a3 , a4 = symbols(('x', 'xa', 'xb', 'a0', 'a1', 'a2', 'a3', 'a4'))
c0 = y.subs(x, xa)
c1 = y.subs(x, xb)
ydiff = diff(y, x)
c2 = ydiff.subs(x, xa)
c3 = ydiff.subs(x, xb)
c4 = ydiff.subs(x, (xa + xb) * 0.5)
c0ArgList = list(c0.args)
c3ArgList = list(c3.args)


x = Symbol('x')
P = x**2 - 3*x + 5
A = Matrix([ [1,3], [-1,2] ])
Poly(P).subs(A)
Poly(Matrix([ [ 1, 3], [-1, 2]])**2 - 3*Matrix([ [ 1, 3],        
    [-1, 2]]) + 5, Matrix([ [ 1, 3], [-1, 2]]), domain='ZZ')

要获取所有系数，请使用all_coeffs方法：

Poly(1+x**2*2+5*x**4).all_coeffs()
[5, 0, 2, 0, 1]
def poly_matrix(P, A):
    coeffs = Poly(P).all_coeffs()[::-1]
    res = zeros(A.rows)
    for i in range(len(coeffs)):
        res += coeffs[i]*(A**i)
    return res

如何在SymPy中的系数列表中创建多项式？

例如，给定一个列表[1, -2, 1]，获得Poly(x**2 - 2*x + 1)
使用Poly.from_list构造多项式：

 x = sympy.Symbol('x')
sympy.Poly.from_list([1, -2, 1], gens=x)
Poly(x**2 - 2*x + 1, x, domain='ZZ')
{x**m[0]:p.coeff_monomial(m) for m in p.monoms()}

from sympy import *
import pprint
x = symbols("x")
y = 4*x**5 + 10*x**2 + x**3 + 7*x + 3
p = Poly(y, x)#多项式子项系数
pprint.pprint(p.coeffs())#多项式子项次数
pprint.pprint(p.monoms())#多项式次数
pprint.pprint(p.degree())
y =(4*x**5 + 10*x**2)/(x**3 + 7*x + 3)
p = fraction(y)
num= Poly(p[0],x).all_coeffs()
dem= Poly(p[1],x).all_coeffs()
A=Poly.from_list(num,gens=x)
B=Poly.from_list(dem,gens=x)
Poly(y).as_expr()

在MATLAB中，复指数函数的调用格式为：

exp((a+j*w)*t)
在Python中的函数表示为：
exp((complex(a,w))*t)

在MATLAB中，矩形脉冲信号可用rectpuls函数产生

，其调用格式为：
y=rectpuls(t,width)
该函数生成幅度为1,宽度为width，以t=0为对称中心的矩形脉冲信号。
Python中用一个分段函数表示矩形脉冲信号。

import numpy as np
import matplotlib.pyplot as plt

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

x=np.linspace(-2,4,1000)
y=np.array([rect_wave(t,2.0,-1.0) for t in x])
plt.ylim(-0.2,1.2)
plt.plot(x,y)
plt.show()

在MATLAB中，

阶跃信号

用“t>=0”产生，调用格式为：
ft=(t>=0)
在Python中可以用where函数绘制其波形，调用格式为：
where(condition，[x,y])
该函数的返回结果是根据前面的条件判断输出x还是y。

import numpy as np
import matplotlib.pyplot as plt
#定义阶跃信号
def unit(t):
     r=np.where(t>0.0,1.0,0.0)
     return r
t=np.linspace(-1.0,3.0,1000)
plt.ylim(-1.0,3.0)
plt.plot(t,unit(t))
plt.show()

from scipy.signal import lti,step2,impulse2
import matplotlib.pyplot as plt

f,((ax1,ax2,ax3),(ax4,ax5,ax6),(ax7,ax8,ax9),(ax10,ax11,ax12)) = plt.subplots(4,3,sharex='col',sharey='row') # 开启subplots模式
for i in range(1,13,1):
    exec("s%s=lti([%d],[1,2,5])"%(i,i))
    exec("t%s,y%s=step2(s%s)"%(i,i,i))
    exec("ax%s.plot(t%s,y%s,'r',label='s%s Step Response',linewidth=0.5)"%(i,i,i,i))
    print (i)
plt.show()

Python 求解系统响应

scipy模块下的子模块signal，函数用step2（）和impulse2（）来分别求阶跃输出和脉冲输出

​

from scipy.signal import lti,step2,impulse2
import matplotlib.pyplot as plt

s1=lti([3],[1,2,10])    # 以分子分母的最高次幂降序的系数构建传递函数，s1=3/(s^2+2s+10）
s2=lti([1],[1,0.4,1])   # s2=1/(s^2+0.4s+1)
s3=lti([5],[1,2,5])     # s3=5/(s^2+2s+5)

t1,y1=step2(s1)         # 计算阶跃输出，y1是Step response of system.
t2,y2=step2(s2)
t3,y3=step2(s3)
t11,y11=impulse2(s1)
t22,y22=impulse2(s2)
t33,y33=impulse2(s3)

f,((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2,3,sharex='col',sharey='row') # 开启subplots模式
ax1.plot(t1,y1,'r',label='s1 Step Response',linewidth=0.5)
ax1.set_title('s1 Step Response',fontsize=9)
ax2.plot(t2,y2,'g',label='s2 Step Response',linewidth=0.5)
ax2.set_title('s2 Step Response',fontsize=9)
ax3.plot(t3,y3,'b',label='s3 Step Response',linewidth=0.5)
ax3.set_title('s3 Step Response',fontsize=9)

ax4.plot(t11,y11,'m',label='s1 Impulse Response',linewidth=0.5)
ax4.set_title('s1 Impulse Response',fontsize=9)
ax5.plot(t22,y22,'y',label='s2 Impulse Response',linewidth=0.5)
ax5.set_title('s2 Impulse Response',fontsize=9)
ax6.plot(t33,y33,'k',label='s3 Impulse Response',linewidth=0.5)
ax6.set_title('s3 Impulse Response',fontsize=9)

plt.show()

import numpy as np
import control as ctl
import matplotlib.pyplot as plt

def step_plot(s):
    t,y=ctl.step_response(s)
    plt.plot(t,y,'b',linewidth=0.6)
    plt.title('Step Response',fontsize=9)
    plt.xlabel('Time(seconds)',fontsize=9)
    plt.ylabel('Amplitude',fontsize=9)
    plt.show()

def impulse_plot(s):
    t,y=ctl.impulse_response(s)
    plt.plot(t,y,'b',linewidth=0.6)
    plt.title('Impulse Response',fontsize=9)
    plt.xlabel('Time(seconds)',fontsize=9)
    plt.ylabel('Amplitude',fontsize=9)
    plt.show()


s=ctl.tf([4],[1,2,10,8])
step_plot(s)

伯德图

from scipy import signal
import matplotlib.pyplot as plt

s1 = signal.lti([2.557e-9,1], [4.999e-15,1.268e-8, 1])
w, mag, phase = s1.bode()
print(w)
plt.figure()
plt.semilogx(w, mag)    # Bode magnitude plot
plt.figure()
plt.semilogx(w, phase)  # Bode phase plot
plt.show()

求解非线性方程组。

from math import cos
from scipy import optimize
def f(x):
    d = 140
    l = 156
    a, r = x.tolist()
    return [
        cos(a) - 1 + (d*d)/(2*r*r),
        l - r * a
    ]
result = optimize.fsolve(f,[1, 1])
print result
print f(result)

离散概率分布

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
x = range(1, 7)
p = (0.4, 0.2, 0.1, 0.1, 0.1, 0.1)
dice = stats.rv_discrete(values=(x, p))  #创建随机变量
np.random.seed(42)
samples = dice.rvs(size=(2000, 50))
sample_mean = np.mean(samples, axis=1)
_, bins, step = plt.hist(sample_mean, bins=100, normed=True, histtype="step",label=u"直方图统计")
'''
_: 直方图向量，是否归一化由参数normed设定

bins: 返回各个bin的区间范围

step: 返回每个bin里面包含的数据，是一个list
'''
kde = stats.kde.gaussian_kde(sample_mean)
x = np.linspace(bins[0], bins[-1], 100)
plt.plot(x, kde(x), label=u"核密度估计")
mean, std = stats.norm.fit(sample_mean)
plt.plot(x, stats.norm(mean, std).pdf(x),alpha=0.8, label=u"正态分布拟合")
plt.legend()
plt.savefig("C:\\figure1.png")

插值积分信号处理

import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
#决定插值曲线的数据点有11个，插值之后的数据点有101个
y = np.sin(x)
xnew = np.linspace(0, 10, 101)
plt.plot(x, y, 'ro')
for kind in ['nearest', 'zero', 'slinear', 'quadratic']:
'''
zero、nearest:阶梯插值
slinear、linear:线性插值
quadratic、cubic:二阶和三阶B样条插值

'''
    f = interpolate.interp1d(x, y, kind=kind)#interpld可以计
    #算x的取值范围内的任意点的函数值
    ynew = f(xnew)
    plt.plot(xnew, ynew, label=str(kind))
plt.legend()
#plt.show()
plt.savefig("c:\\figure1.png")

2. 滤波器
   基于scipy模块使用Python实现简单滤波处理，包括内容有1.低通滤波，2.高通滤波，3.带通滤波，4.带阻滤波器
   scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘ba’)
   sos
   scipy.signal.filtfilt
   scipy.signal.remez函数设计equiripple高通滤波器