基于Bokeh&Matplotlib的时域信号绘制频谱图的简单方案

        笔者在工作中遇到了想在本地获得频谱的问题，通常从示波器上down下来的波形文件剪一剪送进matlib里面仿真一下scope就可以了。

        但是其一不是每台电脑都安装了matlib，这软件的启动也很慢，而且不适应频繁的生成频谱（就是操作有一点点麻烦），导致工作交流，现场使用都不是很方便；其二信号的源数据不一定就是十进制的时域信号（电压幅值），甚至可能不完整（存在拼接或者多余位）需要修剪并转换成快速傅里叶变换所需的数据格式。

        所以笔者最终还是决定设计一个专用的简单方案来处理频谱的生成问题

1、信号的模拟

        实验采用的是Three Tone的信号，分别是xMHz、4x.xMHz、5x.xMHz，基频24xxMHz。

        笔者采用的是perl5脚本来生成，也是顺手选的语言:

        perl没有🥧的概念，所以需要计算一个pi出来，这里用到的是atan2：

my $pi = atan2 0, -1;

        当然，下面的也行：

print 4*atan2(1,1);

         分辨率为1000GHz，即0.001ns（看着唬人，其实更方便理解），采样率为160MHz，生成点100M个，采样点160K个，输出到wfm（其实就是txt）文件（仅供参考）：

#!/usr/bin/perl

print " generation======================================\n";

use File::Find;
$first_tone  = $ARGV[0];
$second_tone = $ARGV[1];
$third_tone  = $ARGV[2];

open(wave_info ,">", "target_waveform_down_conversion.wfm") or die "Cannot open file: $!";

$y=null;
my $pi = atan2 0, -1;
for (my $i = 0; $i <= 99999999; $i++) {
#fs 100g 0.01ns
    $y =  sin( 2 * $pi * 2491300000 * $i * 1e-11) + sin( 2 * $pi * 2492200000 * $i * 1e-11) + sin( 2 * $pi * 2450000000 * $i * 1e-11);
#down_conversion
    $y =  $y - sin( 2 * $pi * 2442000000 * $i * 1e-11);
#fs 80MHz 12.5ns = 1250T; 160MHz 6.25ns = 625T
    if ($i % 625 == 0)
    {
      $y =  sin( 2 * $pi * 49300000 * $i * 1e-12) + sin( 2 * $pi * 50200000 * $i * 1e-12) + sin( 2 * $pi * 8000000 * $i * 1e-12);
      #$y = sin( 2 * $pi * 8000000 * $i * 1e-11);
      print wave_info "$y\n";
      print "\n$y----".$i."-------------$pi";
    }
}

close(wave_info) || die"cannot close file";

        最终效果如下（肯定从0开始嘛，应该没错了）：

        这样一个时域的Three Tone信号就记录下来了

2、频谱图的生成

        时域信号转频域，最重要的就是FFT，具体的理论再此省略，生产中也比较少自建FFT函数的，笔者也是直接使用api。

        至于程序语言，在开发时候考虑到可能使用各种各样的库（其实很简单，并没用多少），包括GUI、读文件、计算等，笔者最终选择了Python，更灵活一些，笔者也更熟悉。

        1、读文件

folder_path = "./"

file_names = glob.glob(folder_path + '*.wfm') 
print(file_names[0])

key_word = []
x_axis = []
i=0

with open(file_names[0], 'r', encoding='utf-8') as f:
    for line in f:
        x_axis.append(i)
        key_word.append(line)
        i=i+1

         2、FFT 

#fft计算后要除掉二分之一个采样点数
fft_v = np.abs(np.fft.fft(key_word))/80000

        3、 Matplotlib：使用 Python 进行可视化

        Matplotlib 是一个综合库，用于在 Python 中创建静态、动画和交互式可视化。Matplotlib 让简单的事情变得简单，让困难的事情成为可能。

链接：Matplotlib 中文网

        python使用的话，pip就可以

# pip安装
$ pip install matplotlib

        时域是将key_word数组直接填入plt中，但是频域首先需要FFT，然后将幅值转换成dbw，其实这也就是一步放大作用，不然我们模拟的信号过于标准，FFT后如果用watt来表示纵坐标，连一点噪声都看不到：

from matplotlib import pyplot as plt

#time domain
x=0
for key in key_word:
    y = float(key)
    plt.plot(x, y, c='k', ls='-',  marker='*', mfc='r', mec='b', mew=1)
    x = x+1
# plt.legend([file_names+"_"], loc="best")
plt.show()
plt.clf()

fft_v = np.abs(np.fft.fft(key_word))/80000
yy = np.multiply(20,np.log10(fft_v)) #dbm


#freq domain
freqs = np.fft.fftfreq(len(key_word), 1/1600000000)

plt.plot(abs(freqs), yy)
plt.xlabel('Hz')

plt.show()
plt.clf()

运行效果：

放大后： 

 频域：

放大后：

        从结果来看，频谱是没问题的，三个峰代表着笔者模拟的三个输入信号。

        4、 Bokeh：更快渲染

        实际运行中可以发现，Matplotlib产出时域图像时候要很久，判断是渲染速度非常之慢导致的，如果只是几十几百个点还好，但是时域波形通常都是上万个。本次模拟波形十六万个点需要足足三十秒才能渲染出来，而且放大缩小等操作也响应的很差，所以笔者决定换一个。

        无论是自身硬件还是Matplotlib本身的问题，既然渲染的慢，那就换个渲染方式，也就是Bokeh。

        Bokeh是一个Python库，用于为现代web浏览器创建交互式可视化。它可以帮助您构建美丽的图形，从简单的绘图到具有流数据集的复杂仪表板。使用Bokeh，您可以创建JavaScript驱动的可视化，而无需自己编写任何JavaScript。 （官方如是介绍到）

        笔者相中其采用浏览器进行渲染的方式，相信Chorem的力量！

Partial code：

from bokeh.plotting import figure, output_file, show
from bokeh.layouts import row, column

#time domain
p = figure(title="Simple Example", x_axis_label="X 轴", y_axis_label="Y 轴", y_range=(-4, 4),)
p.line(x_axis, key_word, legend_label="Test", line_width=2)
show(p)        

#freq domain
p1 = figure(title="freq domain", x_axis_label="Hz 轴", y_axis_label="Y 轴")
# p1.line(abs(freqs), np.absolute(yy), legend_label="Test", line_width=2)
p1.line(freqs, np.real(yy) ,legend_label="Testtest", line_width=2)
show(p1)

时域效果 ：

放大后：

 

 频域效果：

放大后

 

         改用bokeh来做图形化大概提升了五到六倍的速度，大概五秒即可出图。且放大拖动也很丝滑流畅。

         5、 小结

        两种可视化在笔者使用过程中都各有优缺点，但是matplotlib在时域表现得过差，所以上万点数的还是依靠Bokeh更好。如果仅仅看频域，点数少，速度不那么重要的时候，matplotlib可随鼠标取点，获取横纵坐标值会更方便一点。

3、"纵坐标"的转换 

fft_v = np.abs(np.fft.fft(key_word))/80000
yy = np.multiply(20,np.log10(fft_v))
print("voltage", fft_v[xxx])
print("voltagr->dbw", yy[xxx])
watt = np.multiply(fft_v,fft_v) #watt
milliwatt = np.multiply(np.multiply(fft_v, fft_v), 1000)
print("watt", watt[xxx])
print("milliwatt", milliwatt[xxx])
print("watt->dbm", np.multiply(10,np.log10(np.multiply(watt[xxx], 1000))))
print("milliwatt->dbm", np.multiply(10,np.log10(milliwatt[xxx])))

4、特殊类型数据的分割处理

        笔者也会遇到二进制数据的时候，前后也有多余位，我更习惯用正则的方式来处理从文件中顺序读出来的line数据；

        二进制格式数据通常也有最高位判断正负值的功能（振幅总归有正有负），也需要注意；

with open(file_names[0], 'r', encoding='utf-8') as f:
    for line in f:
        mt = re.match(r'\#\#\#read out start.*\=(\d+)', line)
        if mt != None:
            # print(mt.group(1))
            I_bit = int(mt.group(1)[24:34],2)/1024
            Q_bit = int(mt.group(1)[34:44],2)/1024

            #judge postive&negative 
            I_bit = int(mt.group(1)[25:34],2)/1024 * ((-1) ** int(mt.group(1)[24]))
            Q_bit = int(mt.group(1)[35:44],2)/1024 * ((-1) ** int(mt.group(1)[34]))

            print(int(mt.group(1)[24]))
            print(int(mt.group(1)[34]))
            print((-1) ** 0)
            print(-1 ** 0)
            print(-1 ** int(mt.group(1)[24]))
            print(-1 ** int(mt.group(1)[34]))

            print(I_bit) 
            print(Q_bit)

            key_word_I.append(I_bit)
            key_word_Q.append(Q_bit) 

 

方案留存，之后更新吧，当个技术储备

原创文章，欢迎转载，请注明来源，未经书面允许，请勿用于商业用途。