使用Python作时频分析-1

import numpy as np
from matplotlib import pyplot as plt

模拟白噪声

使用均匀分布和标准正太分布函数，rand，randn，生成1000个 0到1的数据以模拟白噪声。

#生成一个长度为10000的，[0,1)区间内的均匀分布和正态分布的随机信号
Yu = np.random.rand(1000,1)
Yn = np.random.randn(1000,1)
#绘制信号
#在一个子图中绘制俩个信号的图像
plt.subplot(211)
plt.hold(True)
plt.plot(Yn,'b')
plt.plot(Yu,'r')
# plt.xticks([i*200 for i in range(6)])
plt.hold(False)

#绘制直方图
plt.subplot(223)
plt.hist(Yu,200)
plt.title("Distribution of uniform noise")
plt.subplot(224)
plt.hist(Yn,200)
plt.title("Distribution of random nosie")
plt.tight_layout()#调整子图间间距和边距，避免重叠

模拟粉色噪声

#生成白噪声
wn = np.random.randn(10000,1)
#计算白噪声声谱
wnX = np.fft.fft(wn) #使用傅里叶变换计算频谱
#计算粉色白噪声频谱，其中np.real是为了得到功率谱，*2是将功率*2，之后进行归一化
N = len(wnX) #获取功率谱长度
#使用IFFT将频谱转换为自相关函数
pn = np.fft.ifft(wnX*(np.linspace(-1,1,N)**2).reshape(N,1)*2)
pn = pn.real

#绘制图像信号
plt.subplot(221)
plt.hold(True)
plt.plot(wn,'b') #绘制白噪声图像
plt.plot(pn,'r') #绘制粉噪声图像
plt.xlabel("Time (a.u.)")
plt.ylabel("Amplitude (a.u.)")
plt.legend(["white","pink"],loc="upper")

#绘制散点图可视化白噪声和粉噪声的关系
plt.subplot(222)
plt.plot(wn,pn,'.')
plt.xlabel("Amplitude white noise")
plt.ylabel("Amplitude pink noise")
plt.title("Scatter Plot")

#绘制功率谱 白噪声和粉噪声
plt.subplot(212)
plt.hold(True)
plt.plot(np.abs(wnX),'r')
plt.plot(np.abs(np.fft.fft(pn)),'b')
plt.xlabel("Frequnecy (a.u.)")
plt.ylabel("Amplitude (a.u.)")
plt.hold(False)
plt.tight_layout()

正弦波

• 正弦波相加&正弦波加入白噪声
• 频率以及振幅变化的正弦波
• 高斯分布函数

正弦波叠加并且引入白噪声

#生成时间向量t 从 0-5 步长为0.0001
t = np.arange(0,5,0.0001)
#计算时间向量t的个数
n = len(t)
#设置振幅
a = [10,2,5,8]
#设置频率
f = [3,1,6,12]
#设置相位
p = [0,np.pi/4,-np.pi,np.pi/2]
#初始化每个时刻正弦波的值
swave = np.zeros(n)
#将每个正弦的值相加
for i in range(len(a)):
    swave += a[i]*np.sin(2*np.pi*f[i]*t+p[i])
#绘制叠加后的正弦波图象
plt.plot(t,swave)
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title("summed Sines")
plt.show()

• 引入白噪声

#加入白噪声
wn = np.random.randn(n)+np.mean(a)
swave += wn
plt.plot(t,swave)
plt.xlabel("Time (s)")
plt.ylabel("Aplitude")
plt.title("Add noise")
plt.show()

频率以及振幅变化的正弦波

#振幅变化的正弦波
#设置设置振幅和频率
a = np.array([10,2,5,8])
f = np.array([3,1,6,12])

t = np.arange(0,5,0.001)
n = len(t)
#将时间分块
block_size = int(n/len(a))
tchunks = np.array([i*block_size for i in range(len(a)+1)])

#初始化每一时刻正弦波的值
swave = np.zeros_like(t)
#循环生成正弦波信号
for i in range(len(a)):
    #计算每个时间块内的正弦波信号值
    signal = a[i]*np.sin(2*np.pi*f[i]*t[tchunks[i]:tchunks[i+1]])
    swave[tchunks[i]:tchunks[i+1]] +=signal
#绘制图像
plt.plot(t,swave)
plt.xlabel("Time (a.u.)")
plt.ylabel("Amplitude (a.u.)")
plt.title("Syntetic Snusidal Waveform")
plt.show()

# 绘制频率不断变化的正弦波
#定义频率数组
f = np.linspace(2,10,n)
#计算正弦波的值
swave = np.sin(2*np.pi*f*t)
#绘制图像
plt.plot(t,swave)
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title("Time-avrying frequency")
plt.show()

# 绘制振幅不断变化的正弦波
#定义频率数组
t = np.arange(-1,5,0.0001)
a = np.linspace(1,10,len(t))
#计算正弦波的值
swave = a*np.sin(2*np.pi*t*3)
#绘制图像
plt.plot(t,swave)
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title("Amplitude-avrying frequency")
plt.show()

高斯分布

t = np.arange(-1,5,0.0001) #定义时间变量-1到5，步长为0.001
s = [0.5,0.1] #标准差，用于控制高斯分布函数的宽度
a = [4,5] # 高斯分布函数的幅值
g1 = a[0]*np.exp(-(t)**2/(2*s[0]**2))
g2 = a[1]*np.exp(-(t-2)**2/(2*s[1]**2))
#绘制图像
plt.plot(t,g1)
plt.hold(True)
plt.plot(t,g2)
plt.legend({"G1","G2"})
plt.show()