深入解构magnitude_spectrum()

2019/08/25 22:56
阅读数 36

深入解构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值作频谱图如下:

img

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)

img

展开阅读全文
打赏
0
0 收藏
分享
加载中
更多评论
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部