深入解构magnitude_spectrum()

二话不说，现附上官方文档链接matplotlib.axes.Axes.magnitude_spectrum

magnitude_spectrum(self, x, Fs=None, Fc=None, window=None, pad_to=None, sides=None, scale=None, , data=None*, **kwargs)

作用：绘制幅度谱

功能：根据x计算幅度谱值，可由pad_to参数指定补零长度，指定截断采用的窗函数

参数：

x: 一维实数或复数向量或序列

Fs：标量值，采样率，默认是2

window：窗函数，默认汉宁窗

sides：枚举值，{'default', 'onesided', 'twosided'}，实数数据默认单边，复数默认双边

pad_to：整数，补零个数，增加频谱密度

scale： 枚举值，{'default', 'linear', 'dB'}，默认linear，而dB amplitude 算法(20 * log10)

Fc：中心频率

返回值：

spectrum: 频谱值，纵坐标

freqs：频率值，横坐标

以下进行数据读取及处理

def load_data():# 加载从USRP采样的IQ信号，load 1024 IQdataIQdata = np.fromfile('usrp_data_with_telosb.dat',dtype="float32",count=1024*2)# merge these data into complex formIQcomplex = map(complex,IQdata[0::2],IQdata[1::2])# change these data into complex typeIQcomplex = np.array(list(IQcomplex),dtype=complex)return IQcomplex

三种方法逐步解析

# 法一：暴力调包
plt.magnitude_spectrum(IQcomplex,Fs=10,Fc=2435,sides='twosided',scale='dB')

# 法二：采用mlab调包
fft_size = 1024
spec, freqs = mlab.magnitude_spectrum(IQcomplex[0:fft_size)], Fs=10, sides='twosided')
Fc = 2435
freqs += Fc
Z = 20.*np.log10(spec)
plt.plot(freqs,Z)

# 法三：用numpy解决，从底层FFT看起
def cal_mag(IQcomplex,fft_size=1024):
#     fft_size = 1024IQcomplex_1024 = np.asarray(IQcomplex[0:fft_size])# 加窗，不然频谱现象不明显result, windowVal = mlab.apply_window(IQcomplex_1024,window=mlab.window_hanning,axis=0,return_window=True)# 进行fftresult = np.fft.fft(result, n=fft_size)# 参数1：FFT点数，参数2：采样周期，即1/采样频率freqs = np.fft.fftfreq(fft_size, 0.1)# 以下代码的作用貌似是调整双边的位置freqcenter = fft_size//2freqs = np.concatenate((freqs[freqcenter:],freqs[:freqcenter]))result = np.concatenate((result[freqcenter:],result[:freqcenter]),0)# 貌似是归一化result = np.abs(result)/np.abs(windowVal).sum()# 加上中心频率Fc = 2435freqs += Fc # 取对数坐标spec = 20.*np.log10(result)return spec,freqsIQcomplex = load_data()
mag, freqs = cal_mag(IQcomplex)
plt.plot(freqs,mag)

USRP采样所得的IQ值作频谱图如下：

IQcomplex = load_data()
plt.figure(figsize=(6,10))
num = 20
for n in range(10):plt.subplot(5,2,n+1)IQcomplex_tmp = IQcomplex[(n+num)*1024:(n+1+num)*1024]mag, freqs = cal_mag(IQcomplex_tmp)plt.plot(freqs,mag)
#     plt.title('{}'%n)