2019/08/25 22:56

# 深入解构magnitude_spectrum()

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

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

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

window：窗函数，默认汉宁窗

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

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

Fc：中心频率

spectrum: 频谱值，纵坐标

freqs：频率值，横坐标

def load_data():
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

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)


0
0 收藏

0 评论
0 收藏
0