Python分析系统的时域特性和频率域特性

在不使用matlab的情况下,可以选择用python来实现自动控制理论有关系统打时域分析和频率域分析等,安装的package是python-control,在windows的控制台(cmd)或者linux终端下输入pip install control 即可,注意,如果同时安装了2.73.x(3.4或者3.5或者3.6 版本,使用pip 命令打时候需要指定版本号,如pip2 install control 或者pip3.4 install control ,当然,常用打科学计算用的package也要安装,numpy,scipy,sympy,matplotlib,pandas 等。

下面是自己练习时写的代码,写在此作记录和分享用,因为函数语法和matlab相差无几,这里就没有写太多的注释了,有需要打话可以去python-control打官网查看相关文档。

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 16 17:48:02 2016

@author: kindy
"""

from control import *
from scipy import signal as sgl
from matplotlib import pyplot as plt
import numpy as np

T=np.mgrid[0:8:0.02]
U1=T
U2=T**2


sys1 = tf([1],[0.5,1]) # 
sys2 = tf([2],[1,2,4]) # 

# Step Response
def step_resp():
    sout1,stime1 = step(sys1)
    sout2,stime2 = step(sys2)
    plt.plot(stime1,sout1,'b',linewidth=0.5)
    plt.plot(stime2,sout2,'b',linewidth=0.5)
    plt.xlabel("Time")
    plt.ylabel("Amplitude")
    plt.title("Step Resopnse",fontsize=12)
    #plt.legend()
    plt.show()

# Impulse Response
def impulse_resp():

    iout1,itime1 = impulse(sys1)
    iout2,itime2 = impulse(sys2)
    plt.plot(itime1,iout1,'m',linewidth=0.8)
    plt.plot(itime2,iout2,'r',linewidth=0.8)
    plt.show()

#impulse_resp()

# 任意输入信号的输出,lsim
def lsim_plot():
    yout1,Time1, xout1 = lsim(sys1, U1, T)
    yout2,Time2, xout2 = lsim(sys2, U1, T)
    plt.plot(Time1, yout1, 'b', linewidth=0.7)
    plt.plot(Time2, yout2, 'b', linewidth=0.7)
    plt.show()

#lsim_plot()

# 波特图
def bode_plot():
    bode(sys1)
    bode(sys2)

#bode_plot()

# Nyquist图
def nyquist_plot():
    nyquist(sys1)
    nyquist(sys2)

#nyquist_plot()  


# 根轨迹
def root_locus():
    rlocus(sys1)
    rlocus(sys2)

root_locus()


下面是运行打一些结果图:
这里写图片描述

这里写图片描述

这里写图片描述

这里写图片描述

时域分析是信号处理和控制系统中常用的一种分析方法,可以用来分析信号的时域特性,如幅度、相位、周期、频率等。在Python中,可以使用numpy、scipy、matplotlib等库进行时域分析。 1. 采样和时间序列生成 在Python中,可以使用numpy库生成时间序列,例如: ```python import numpy as np import matplotlib.pyplot as plt # 采样频率 fs = 1000 # 采样时间 t = np.arange(0, 1, 1/fs) # 生成正弦波信号 f = 10 # 信号频率 x = np.sin(2*np.pi*f*t) # 绘制信号图像 plt.plot(t, x) plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.show() ``` 上述代码中,np.arange()函数生成了一个时间序列,从0到1,步长为1/fs,即采样周期。然后通过np.sin()函数生成了一个频率为10Hz的正弦波信号,并使用matplotlib库绘制了信号的时域波形图。 2. 时域特性分析Python中,可以使用scipy库中的函数进行时域特性分析,如: - 平均值和方差 ```python import numpy as np from scipy.stats import describe # 生成随机信号 x = np.random.rand(100) # 计算信号的平均值和方差 mean, var, skew, kurt = describe(x) print('Mean:', mean) print('Variance:', var) ``` - 峰值和峰峰值 ```python import numpy as np # 生成信号 x = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1]) # 计算信号的峰值和峰峰值 peak = np.max(x) peak_to_peak = np.max(x) - np.min(x) print('Peak:', peak) print('Peak to peak:', peak_to_peak) ``` - 自相关函数和互相关函数 ```python import numpy as np from scipy.signal import correlate # 生成信号 x = np.array([1, 2, 3, 4, 5]) y = np.array([1, 0, 1, 0, 1]) # 计算信号的自相关函数和互相关函数 corr_x = correlate(x, x, mode='same') corr_xy = correlate(x, y, mode='same') print('Auto-correlation:', corr_x) print('Cross-correlation:', corr_xy) ``` 上述代码中,使用scipy库中的correlate()函数计算信号的自相关函数和互相关函数。其中,mode参数表示相关函数的计算方式,具体有:'valid'、'same'和'full'三种模式。 3. 时域滤波 在Python中,可以使用scipy库中的函数进行时域滤波,如: - 移动平均滤波 ```python import numpy as np from scipy.signal import lfilter # 生成信号 x = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1]) # 定义滤波器系数 N = 3 b = np.ones(N) / N # 进行滤波 y = lfilter(b, 1, x) print('Filtered signal:', y) ``` 上述代码中,使用lfilter()函数进行了移动平均滤波。其中,b参数为滤波器系数,1为滤波器的阶数,x为原始信号。 - 中值滤波 ```python import numpy as np from scipy.signal import medfilt # 生成信号 x = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1]) # 进行中值滤波 y = medfilt(x) print('Filtered signal:', y) ``` 上述代码中,使用medfilt()函数进行了中值滤波,无需定义滤波器系数,直接输入原始信号即可。 4. 傅里叶变换 在Python中,可以使用numpy库和scipy库中的函数进行傅里叶变换,如: ```python import numpy as np import matplotlib.pyplot as plt # 采样频率 fs = 1000 # 采样时间 t = np.arange(0, 1, 1/fs) # 生成正弦波信号 f = 10 # 信号频率 x = np.sin(2*np.pi*f*t) # 进行傅里叶变换 X = np.fft.fft(x) # 计算频率分辨率 df = fs / len(x) # 生成频率轴 freqs = np.arange(0, fs, df) # 绘制频谱图 plt.plot(freqs[:len(x)//2], np.abs(X[:len(x)//2])/len(x)) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude') plt.show() ``` 上述代码中,使用numpy库中的fft()函数进行了傅里叶变换,并使用matplotlib库绘制了信号的频谱图。其中,np.abs()函数计算了信号的幅度谱,len()函数用于计算信号长度,freqs为频率轴。 5. 窗函数 在进行傅里叶变换时,可以使用窗函数对信号进行加窗处理,以减少频谱泄漏等问题。在Python中,可以使用numpy库中的函数生成窗函数,如: ```python import numpy as np import matplotlib.pyplot as plt # 采样频率 fs = 1000 # 采样时间 t = np.arange(0, 1, 1/fs) # 生成正弦波信号 f = 10 # 信号频率 x = np.sin(2*np.pi*f*t) # 定义窗函数 win = np.hanning(len(x)) # 进行傅里叶变换 X = np.fft.fft(x*win) # 计算频率分辨率 df = fs / len(x) # 生成频率轴 freqs = np.arange(0, fs, df) # 绘制频谱图 plt.plot(freqs[:len(x)//2], np.abs(X[:len(x)//2])/len(x)) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude') plt.show() ``` 上述代码中,使用numpy库中的hanning()函数生成了一个汉宁窗,并将其应用于信号加窗处理。然后,使用fft()函数进行傅里叶变换,并使用matplotlib库绘制了信号的频谱图。
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值