深入解构magnitude_spectrum()

深入解构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 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 = 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 = 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,freqs

IQcomplex = 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)