0. 本文目的
- 主要辅助信号与系统课程的学习与讲解。
- 一些零碎的记录，以后再酌情细化。
1. 利用sympy库
基本方式是直接定义函数，之后调用自带的函数做处理、计算、变换，以及调用plot函数画图。
- 可以描述连续函数波形。可以描述冲激信号（DiracDelta）、阶跃信号（Heaviside）、符号函数（sign）等多种函数的波形。
- 可以波形变换（subs）、求导、定积分（integrate）、卷积（好像没看到，但可以用定积分定义）、傅里叶变换（fourier_transform）、得到复函数的幅度相位等、以及微分方程求通解等。
- 可以将表达式转为numpy（lambdify），这个感觉很方便。
一些坑：
- 如何描述离散信号，这个还没来得及看，其实有个离散的模块，好像也有FFT之类的东西，挺全的。但感觉怎么也不会比numpy更方便，暂时不太打算看了。
- 自带的plot方法画图，感觉很弱，目前没有搞定的几个点：
- 如何画虚线和网格，没看到官方的例子。
- 如何嵌入到matplotlib，网上有一些教程，有点麻烦，还没来得及看，但感觉不值得。因为可以将sympy表达式转为numpy序列之后再画图。
- sympy里的sinc函数实际是Sa函数（不是Sa(pi*t)），文档里到是也说了；numpy当中的sinc就是真正的Sa函数。
- 尝试利用定积分函数integrate和符号替换（subs）定义连续信号的卷积，但很受限制，主要问题：
- 很慢，感觉上这东西是“数学”工具（吗？），不能做工程工具。
- 一些复杂函数积分的时候会卡死，懒得总结规律了。有个奇葩的例子：
x = Heaviside(t)-Heaviside(t-0.1)
h = sinc(pi * 2 * 25 * t)
#下一步做卷积有点慢，plot画图的时候卡死
x = Heaviside(t)-Heaviside(t-0.1*pi)
h = sinc(2 * 25 * t)
#一切正常，还有点快
- 会出现一些很复杂的积分结果，无法化简，但有时候化简得就特别完美。但特别慢。
- 积分上下限而可能需要根据计算结果调整：1，计算方波卷积sin区间的时候，上下限用无穷，可以得到一个很简化的表达式，基本是最简形式。如果自行设定一个足够大的上下限，虽然图像正确，但表达式超级复杂，simplify也简化不了；2，如果计算方波卷积方波，则给一个小的上下限比较好。积分限给正负无穷，则表达式无法化简、画图卡死。
- fourier_transforms对复杂函数也很难处理，情况和上面积分、卷积的情况类似，比如基本的sin、cos信号的频谱等（可能是因为包频谱还有冲激形式）。
2. 利用numpy
总的来说挺强大的，啥都能干。
一些笔记：
- 描述连续信号的时候会觉得有点奇怪——只说描述方式，不考虑底层原理。sympy描述的是一个连续的表达式，但numpy代码描述的是一个序列（这个实际也算不上坑）。
- 做FFT的时候，出现过相位谱混乱的情况，问题和解决方案：
#产生横坐标序列，如果采用np.arange，长度需要加入一个sample_interval，否则相位图会混乱，可能和arange的特性有关
t=np.arange(0,10+sample_interval,sample_interval)
#如果采用linspace则没有问题，注意最后一个参数是采样个数
t=np.linspace(0,10,int(10/sample_interval))
但总的来说，FFT相位谱的效果，如果仅利用基本的傅里叶变换原理不是很好解释。
- 对于Sa函数（截短），有两种定义方式，即sint/t和直接使用sinc函数：
t1=np.linspace(-10,10,int(20/sample_interval))
np.seterr(divide='ignore',invalid='ignore')#忽略除以0的警告
ht1 =np.sin( 2 * np.pi * 25 * t1)/(2 * np.pi * 25 * t1)
ht1[0]=1 #特殊点直接赋值了
ht2 = np.sinc( 2 * 25 * t1)
#看一下相减不为零的个数，对比序列长度（6875 和20480）
ht3=(ht1-ht2)
print(np.count_nonzero(ht3),len(ht3))
这两个定义方式出来的序列是不一样的，可能和np.pi的精度有关，但用起来效果基本相同。
sint/t的顶部放大（sinc的顶部是平的）：
- 模拟连续信号卷积，卷积结果存在幅度问题，实际就是一个连续时间卷积的数值近似方法问题，这里算是换了个方式推导一下。
按照连续卷积的定义，卷积结果的最大值应该是2，参见sympy的结果：
原因：
numpy（以及scipy）的卷积是离散卷积，卷积结果的最大值和序列的个数有关。例如：
[1,1]*[1,1]=[1,2,1]
[1,1,1]*[1,1,1]=[1,2,3,2,1]
注意上述两个卷积的形状是一样的，但高度不同。
用numpy模拟连续信号时，模拟信号的长度和实际数字序列的长度不是一回事：
t=np.arange(-10,10,delta)
t=np.linspace(-10,10,sample_amout) #sample_amout = 20 * delta
模拟信号的长度是20，数字序列的长度是10*delta=sample_amout。区间一定的情况下，序列个数越多，则离散卷积的幅度越大，因此要进行修正。修正方式：
模拟卷积结果就是：convolve 乘以 delta，公式如下：
或者convolve结果 乘以两个信号的模拟长度 ，再除以两个信号的序列个数。
-
扩展问题，如果卷积时两个信号的抽样率不同，则卷积结果可能有误，例如：
定义两个u(t)-u(t-1)的方波，理论上卷积应该得到一个从0到2的三角波。
但如果第二个方波的采样率是前者的一倍，则可以理解为是一个窄方波卷积一个宽方波，此时会得到一个梯形。
image.png
类似的问题结合理论分析应该还能找到其他的吧，但并不难解决。
画信号频谱图有两种方式：
- 直接调用matplotlib的 magnitude_spectrum和phase_spectrum，输入时域的信号（signal）即可作图，最好结合实际给一个Fs抽样间隔参数（比如：1/信号长度），别的参数都是可选项，能省大概5、6条语句吧，但貌似坐标轴的label不能调，可能是没找对地方吧。
plt.magnitude_spectrum(signal, Fs=sample_freq )
plt.phase_spectrum(signal, Fs = sample_freq)
- 手动FFT，再用matplotlib画折线图，需要自己求幅度相位、自己搞定坐标、自己做归一化。
好处是自己画图可以做更多的美化，但自己做FFT的
时候需要注意一个小问题：
fft_data= fft(et)
fft_amp= np.abs(fft_data) * deltaT # 双边幅度谱
xAxis = fftfreq(len(fft_data1), deltaT )#fft的双边频域坐标
fft_amp1 = fftshift(fft_amp1)#修改坐标范围，方便画图
xAxis = fftshift(xAxis ) #修改坐标范围，方便画图
fftfreq和fft语句，输出的数据和坐标顺序为：[0到正坐标最大值，负坐标最大值到0]，即不是从小到大排列的。者在画图的时候会出现一些奇怪的线，从正坐标最大值一直连到负的坐标最大值，
例如：
x = [1,2,3,4,5,-5,-4,-3,-2,-1] #x轴
y = [1,2,3,4,5, 6, 7, 8, 9, 10] #数据
出来的图为：
解决这个问题，需要执行fftshift语句，将数据的坐标顺序改为从小到达，相当于改为：
x = [1,2,3,4,5,-5,-4,-3,-2,-1]
y = [1,2,3,4,0, 0, 7, 8, 9, 10]
3. 利用Sympy再转numpy
利用sympy定义连续函数表达式，可以处理根据需要做函数变换、卷积、处理等，得到解析解，然后再转到numpy序列，再画图。
n = np.arange(-3, 3,0.01)
ft = Heaviside(t+1)-Heaviside(t-1) #sympy表达式
ftn = lambdify(t, ft, "numpy") #转为可执行函数，并且是numpy类型的函数
r = ftn(n) #np序列
或者
n = np.arange(-3, 3,0.01)
r = np.array([i*(Heaviside(i)-Heaviside(i-1)) for i in n])
4. 利用SCIPY
一是能做FFT这些，但numpy也能做，简单使用都差不多。
二是处理s域或z域的系统函数，能描述函数、求零极点（tf2zpk）、状态方程和特征值（tf2ss、eig、eigh）、求频响特性（freqs，freqz）
三是能做2D卷积，图像处理等。numpy下暂时没找到现成的二维卷积方法。
5. control库
主要是系统函数的描述（tf）、求零极点（pzmap）、频响特性（bode_plot）。方便之处是直接可以画图，比如z域零极点图可以一句话出单位圆和零极点。
