import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
'''
思路:1)先利用fft计算得出其幅频值2)在利用rfft计算得出其幅频值,看1)和2)那个能还原出信号的原始幅值
'''
# 生成一个示例信号
np.random.seed(0)
fs = 512 # 采样频率
N = 512 # 信号长度
time = np.arange(N) / fs
sig = 10 * np.sin(2 * np.pi * 50 * time) + 5 * np.random.randn(N)############ 1)利用fft计算得出其幅频值 ############
fft_signal = fft(sig)
fft_p = np.abs(fft_signal/N)
# 计算单边频谱
fft_p2 = fft_p[:N//2 + 1]
fft_p2[1:-1] = 2*fft_p2[1:-1]
fft_freqs = np.arange(0, N//2 + 1)*(fs / N)############ 2)利用rfft计算得出其幅频值 ############
rfft_signal = rfft(sig)
rfft_p = np.abs(rfft_signal/N)
rfft_p2 = rfft_p
rfft_p2[1:-1] = 2*rfft_p2[1:-1]plt.figure(1)
plt.subplot(211)
plt.plot(fft_freqs, fft_p2)
plt.title("fft method")
plt.subplot(212)
plt.plot(fft_freqs, rfft_p2)
plt.title("rfft method")plt.show()
结论:利用fft方法和rfft方法皆可以反映出信号的真实幅值大小
############ 3)利用welch计算得出其功率谱 ############
f_welch, Pxx_welch = welch(sig, fs, nperseg=N)############ 4)利用rfft计算得出其功率谱 ############
win = signal.get_window('hann', N) # 获取窗函数
det_sig = detrend(sig, type='constant') # 去除线性趋势
sig_win = win*det_sig # 加窗操作
rfft_win = rfft(sig_win) # 快速傅里叶变换,获取单边频谱rfft_abs = np.conjugate(rfft_win)*rfft_win # 模的平方(此处的rfft_abs为复数)
# rfft_abs = abs(rfft_win)**2scale = 1.0 / (fs * (win*win).sum()) # 窗函数的缩放因子
Pxx_rfft = rfft_abs * scale # 乘于窗函数的缩放因子,保持加窗后的能量不变
Pxx_rfft[1:-1] = 2 * Pxx_rfft[1:-1] # 获取双边频谱值
Pxx_rfft = Pxx_rfft.real # 获取实部plt.figure(2)
plt.subplot(211)
plt.plot(f_welch, Pxx_welch)
plt.title('welch method')plt.subplot(212)
plt.plot(f_welch, Pxx_rfft)
plt.title('rfft method')plt.show()
结论:利用welch方法实现的功率谱和利用rfft方法实现的功率谱结果一致。只所以welch方法未能在数据上反映出真是数值,是因为没有除于样本点个数进行归一化。