深入解构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 IQdata IQdata = np.fromfile('usrp_data_with_telosb.dat',dtype="float32",count=1024*2) # merge these data into complex form IQcomplex = map(complex,IQdata[0::2],IQdata[1::2]) # change these data into complex type IQcomplex = np.array(list(IQcomplex),dtype=complex) return IQcomplex 三种方法逐步解析
# 法一：暴力调包plt.magnitude_spectrum(IQcomplex,Fs=10,Fc=2435,sides='twosided',scale='dB')
# 法二：采用mlab调包fft_size = 1024spec, freqs = mlab.magnitude_spectrum(IQcomplex[0:fft_size)], Fs=10, sides='twosided')Fc = 2435freqs += Fc Z = 20.*np.log10(spec)plt.plot(freqs,Z)
# 法三：用numpy解决，从底层FFT看起def cal_mag(IQcomplex,fft_size=1024):# fft_size = 1024 IQcomplex_1024 = np.asarray(IQcomplex[0:fft_size]) # 加窗，不然频谱现象不明显 result, windowVal = mlab.apply_window(IQcomplex_1024,window=mlab.window_hanning,axis=0,return_window=True) # 进行fft result = np.fft.fft(result, n=fft_size) # 参数1：FFT点数，参数2：采样周期，即1/采样频率 freqs = np.fft.fftfreq(fft_size, 0.1) # 以下代码的作用貌似是调整双边的位置 freqcenter = fft_size//2 freqs = np.concatenate((freqs[freqcenter:],freqs[:freqcenter])) result = np.concatenate((result[freqcenter:],result[:freqcenter]),0) # 貌似是归一化 result = np.abs(result)/np.abs(windowVal).sum() # 加上中心频率 Fc = 2435 freqs += 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 = 20for 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) 