0. 本文目的
- 主要輔助信號(hào)與系統(tǒng)課程的學(xué)習(xí)與講解持偏。
- 一些零碎的記錄，以后再酌情細(xì)化氨肌。
1. 利用sympy庫(kù)
基本方式是直接定義函數(shù)鸿秆，之后調(diào)用自帶的函數(shù)做處理、計(jì)算怎囚、變換卿叽，以及調(diào)用plot函數(shù)畫(huà)圖。
- 可以描述連續(xù)函數(shù)波形恳守「矫保可以描述沖激信號(hào)（DiracDelta）、階躍信號(hào)（Heaviside）井誉、符號(hào)函數(shù)（sign）等多種函數(shù)的波形。
- 可以波形變換（subs）整胃、求導(dǎo)颗圣、定積分（integrate）、卷積（好像沒(méi)看到屁使，但可以用定積分定義）在岂、傅里葉變換（fourier_transform）、得到復(fù)函數(shù)的幅度相位等蛮寂、以及微分方程求通解等蔽午。
- 可以將表達(dá)式轉(zhuǎn)為numpy（lambdify），這個(gè)感覺(jué)很方便酬蹋。
一些坑：
- 如何描述離散信號(hào)及老，這個(gè)還沒(méi)來(lái)得及看抽莱，其實(shí)有個(gè)離散的模塊，好像也有FFT之類(lèi)的東西骄恶，挺全的食铐。但感覺(jué)怎么也不會(huì)比numpy更方便，暫時(shí)不太打算看了僧鲁。
- 自帶的plot方法畫(huà)圖虐呻，感覺(jué)很弱，目前沒(méi)有搞定的幾個(gè)點(diǎn)：
- 如何畫(huà)虛線和網(wǎng)格寞秃，沒(méi)看到官方的例子斟叼。
- 如何嵌入到matplotlib，網(wǎng)上有一些教程春寿，有點(diǎn)麻煩朗涩，還沒(méi)來(lái)得及看，但感覺(jué)不值得堂淡。因?yàn)榭梢詫ympy表達(dá)式轉(zhuǎn)為numpy序列之后再畫(huà)圖馋缅。
- sympy里的sinc函數(shù)實(shí)際是Sa函數(shù)（不是Sa(pi*t)），文檔里到是也說(shuō)了绢淀；numpy當(dāng)中的sinc就是真正的Sa函數(shù)萤悴。
- 嘗試?yán)枚ǚe分函數(shù)integrate和符號(hào)替換（subs）定義連續(xù)信號(hào)的卷積，但很受限制皆的，主要問(wèn)題：
- 很慢覆履，感覺(jué)上這東西是“數(shù)學(xué)”工具（嗎？）费薄，不能做工程工具硝全。
- 一些復(fù)雜函數(shù)積分的時(shí)候會(huì)卡死，懶得總結(jié)規(guī)律了楞抡。有個(gè)奇葩的例子：
x = Heaviside(t)-Heaviside(t-0.1)
h = sinc(pi * 2 * 25 * t)
#下一步做卷積有點(diǎn)慢伟众，plot畫(huà)圖的時(shí)候卡死
x = Heaviside(t)-Heaviside(t-0.1*pi)
h = sinc(2 * 25 * t)
#一切正常，還有點(diǎn)快
- 會(huì)出現(xiàn)一些很復(fù)雜的積分結(jié)果召廷，無(wú)法化簡(jiǎn)凳厢，但有時(shí)候化簡(jiǎn)得就特別完美。但特別慢竞慢。
- 積分上下限而可能需要根據(jù)計(jì)算結(jié)果調(diào)整：1先紫，計(jì)算方波卷積sin區(qū)間的時(shí)候，上下限用無(wú)窮筹煮，可以得到一個(gè)很簡(jiǎn)化的表達(dá)式遮精，基本是最簡(jiǎn)形式。如果自行設(shè)定一個(gè)足夠大的上下限败潦，雖然圖像正確本冲，但表達(dá)式超級(jí)復(fù)雜准脂，simplify也簡(jiǎn)化不了；2，如果計(jì)算方波卷積方波，則給一個(gè)小的上下限比較好沸停。積分限給正負(fù)無(wú)窮晕翠，則表達(dá)式無(wú)法化簡(jiǎn)、畫(huà)圖卡死。
- fourier_transforms對(duì)復(fù)雜函數(shù)也很難處理，情況和上面積分、卷積的情況類(lèi)似院塞，比如基本的sin、cos信號(hào)的頻譜等（可能是因?yàn)榘l譜還有沖激形式）性昭。
2. 利用numpy
總的來(lái)說(shuō)挺強(qiáng)大的拦止，啥都能干。
一些筆記：
- 描述連續(xù)信號(hào)的時(shí)候會(huì)覺(jué)得有點(diǎn)奇怪——只說(shuō)描述方式糜颠，不考慮底層原理汹族。sympy描述的是一個(gè)連續(xù)的表達(dá)式，但numpy代碼描述的是一個(gè)序列（這個(gè)實(shí)際也算不上坑）其兴。
- 做FFT的時(shí)候顶瞒，出現(xiàn)過(guò)相位譜混亂的情況，問(wèn)題和解決方案：
#產(chǎn)生橫坐標(biāo)序列元旬，如果采用np.arange榴徐，長(zhǎng)度需要加入一個(gè)sample_interval，否則相位圖會(huì)混亂匀归，可能和arange的特性有關(guān)
t=np.arange(0,10+sample_interval,sample_interval)
#如果采用linspace則沒(méi)有問(wèn)題坑资，注意最后一個(gè)參數(shù)是采樣個(gè)數(shù)
t=np.linspace(0,10,int(10/sample_interval))
但總的來(lái)說(shuō)，F(xiàn)FT相位譜的效果穆端，如果僅利用基本的傅里葉變換原理不是很好解釋袱贮。
- 對(duì)于Sa函數(shù)（截短），有兩種定義方式体啰，即sint/t和直接使用sinc函數(shù)：
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 #特殊點(diǎn)直接賦值了
ht2 = np.sinc( 2 * 25 * t1)
#看一下相減不為零的個(gè)數(shù)攒巍，對(duì)比序列長(zhǎng)度（6875 和20480）
ht3=(ht1-ht2)
print(np.count_nonzero(ht3),len(ht3))
這兩個(gè)定義方式出來(lái)的序列是不一樣的，可能和np.pi的精度有關(guān)狡赐，但用起來(lái)效果基本相同。
sint/t的頂部放大（sinc的頂部是平的）：
- 模擬連續(xù)信號(hào)卷積钦幔，卷積結(jié)果存在幅度問(wèn)題枕屉，實(shí)際就是一個(gè)連續(xù)時(shí)間卷積的數(shù)值近似方法問(wèn)題，這里算是換了個(gè)方式推導(dǎo)一下鲤氢。
按照連續(xù)卷積的定義搀擂，卷積結(jié)果的最大值應(yīng)該是2西潘，參見(jiàn)sympy的結(jié)果：
原因：
numpy（以及scipy）的卷積是離散卷積，卷積結(jié)果的最大值和序列的個(gè)數(shù)有關(guān)哨颂。例如：
[1,1]*[1,1]=[1,2,1]
[1,1,1]*[1,1,1]=[1,2,3,2,1]
注意上述兩個(gè)卷積的形狀是一樣的喷市，但高度不同。
用numpy模擬連續(xù)信號(hào)時(shí)威恼，模擬信號(hào)的長(zhǎng)度和實(shí)際數(shù)字序列的長(zhǎng)度不是一回事：
t=np.arange(-10,10,delta)
t=np.linspace(-10,10,sample_amout) #sample_amout = 20 * delta
模擬信號(hào)的長(zhǎng)度是20品姓，數(shù)字序列的長(zhǎng)度是10*delta=sample_amout。區(qū)間一定的情況下箫措，序列個(gè)數(shù)越多腹备，則離散卷積的幅度越大，因此要進(jìn)行修正斤蔓。修正方式：
模擬卷積結(jié)果就是：convolve 乘以 delta植酥，公式如下：
或者convolve結(jié)果 乘以兩個(gè)信號(hào)的模擬長(zhǎng)度 ，再除以兩個(gè)信號(hào)的序列個(gè)數(shù)弦牡。
-
擴(kuò)展問(wèn)題友驮，如果卷積時(shí)兩個(gè)信號(hào)的抽樣率不同，則卷積結(jié)果可能有誤驾锰，例如：
定義兩個(gè)u(t)-u(t-1)的方波卸留，理論上卷積應(yīng)該得到一個(gè)從0到2的三角波。
但如果第二個(gè)方波的采樣率是前者的一倍稻据，則可以理解為是一個(gè)窄方波卷積一個(gè)寬方波艾猜，此時(shí)會(huì)得到一個(gè)梯形。
image.png
類(lèi)似的問(wèn)題結(jié)合理論分析應(yīng)該還能找到其他的吧捻悯，但并不難解決匆赃。
畫(huà)信號(hào)頻譜圖有兩種方式：
- 直接調(diào)用matplotlib的 magnitude_spectrum和phase_spectrum，輸入時(shí)域的信號(hào)（signal）即可作圖今缚，最好結(jié)合實(shí)際給一個(gè)Fs抽樣間隔參數(shù)（比如：1/信號(hào)長(zhǎng)度）算柳，別的參數(shù)都是可選項(xiàng)，能省大概5姓言、6條語(yǔ)句吧瞬项，但貌似坐標(biāo)軸的label不能調(diào)，可能是沒(méi)找對(duì)地方吧何荚。
plt.magnitude_spectrum(signal, Fs=sample_freq )
plt.phase_spectrum(signal, Fs = sample_freq)
- 手動(dòng)FFT囱淋，再用matplotlib畫(huà)折線圖，需要自己求幅度相位餐塘、自己搞定坐標(biāo)妥衣、自己做歸一化。
好處是自己畫(huà)圖可以做更多的美化，但自己做FFT的
時(shí)候需要注意一個(gè)小問(wèn)題：
fft_data= fft(et)
fft_amp= np.abs(fft_data) * deltaT # 雙邊幅度譜
xAxis = fftfreq(len(fft_data1), deltaT )#fft的雙邊頻域坐標(biāo)
fft_amp1 = fftshift(fft_amp1)#修改坐標(biāo)范圍税手，方便畫(huà)圖
xAxis = fftshift(xAxis ) #修改坐標(biāo)范圍蜂筹，方便畫(huà)圖
fftfreq和fft語(yǔ)句，輸出的數(shù)據(jù)和坐標(biāo)順序?yàn)椋篬0到正坐標(biāo)最大值芦倒，負(fù)坐標(biāo)最大值到0]艺挪，即不是從小到大排列的。者在畫(huà)圖的時(shí)候會(huì)出現(xiàn)一些奇怪的線兵扬，從正坐標(biāo)最大值一直連到負(fù)的坐標(biāo)最大值麻裳，
例如：
x = [1,2,3,4,5,-5,-4,-3,-2,-1] #x軸
y = [1,2,3,4,5, 6, 7, 8, 9, 10] #數(shù)據(jù)
出來(lái)的圖為：
解決這個(gè)問(wèn)題，需要執(zhí)行fftshift語(yǔ)句周霉，將數(shù)據(jù)的坐標(biāo)順序改為從小到達(dá)掂器，相當(dāng)于改為：
x = [1,2,3,4,5,-5,-4,-3,-2,-1]
y = [1,2,3,4,0, 0, 7, 8, 9, 10]
3. 利用Sympy再轉(zhuǎn)numpy
利用sympy定義連續(xù)函數(shù)表達(dá)式，可以處理根據(jù)需要做函數(shù)變換俱箱、卷積国瓮、處理等，得到解析解狞谱，然后再轉(zhuǎn)到numpy序列乃摹，再畫(huà)圖。
n = np.arange(-3, 3,0.01)
ft = Heaviside(t+1)-Heaviside(t-1) #sympy表達(dá)式
ftn = lambdify(t, ft, "numpy") #轉(zhuǎn)為可執(zhí)行函數(shù)跟衅，并且是numpy類(lèi)型的函數(shù)
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也能做，簡(jiǎn)單使用都差不多伶跷。
二是處理s域或z域的系統(tǒng)函數(shù)掰读，能描述函數(shù)、求零極點(diǎn)（tf2zpk）叭莫、狀態(tài)方程和特征值（tf2ss蹈集、eig、eigh）雇初、求頻響特性（freqs拢肆，freqz）
三是能做2D卷積，圖像處理等靖诗。numpy下暫時(shí)沒(méi)找到現(xiàn)成的二維卷積方法郭怪。
5. control庫(kù)
主要是系統(tǒng)函數(shù)的描述（tf）、求零極點(diǎn)（pzmap）刊橘、頻響特性（bode_plot）鄙才。方便之處是直接可以畫(huà)圖，比如z域零極點(diǎn)圖可以一句話出單位圓和零極點(diǎn)促绵。
