前言
今天要介绍的是一个C语言编写的随机数库Radiom，它需要RTLSDR硬件设备支持。目前仅能在Linux平台使用，不过可以很方便地修改，使其能够在Windows和macOS上生成。其实一个月以前就在做这个工作了，毕竟那一段时间比较闲，想要找一点事情做。不过基本完工之后呢，又渐渐地忙了起来，所以这里其实是写的一个月前的事情。
两年以前，我就写过一篇《软件定义无线电导论》，简单介绍了一些廉价的 RealTeK SDR 设备的玩法。一个月前经常纠结于外卖点什么的问题，经常使用RANDOM.ORG的服务。那么也算是心血来潮，我想，为什么不可以自己做一个相似的东西呢？想到手边有一个RTL2832的电视棒，觉得可行，于是说干就干了。
技术验证
librtlsdr库是一切的基石，各种各样的应用都需要它来支撑，它自身也被许多人制作了各种各样语言的binding。然而不那么美妙的是，它本身的文档工作做得不好，难以下手。我选择从dump1090的源代码切入，虽然它是一份数千行的C代码，但是真正调用到librtlsdr的地方并不是那么多，稍加阅读可以总结出一部分API的用法。
那么简单学习之后，我已经能够使用librtlsdr从电波中读出字节流了。但是原始的无线电信号完全不能满足随机性要求，怎么办呢？在GitHub上，我找到了一个思路与这个接近的项目：REDOUBLER，他的做法是构造一个ChaCha20，然后把SDR的字节流送进去，出来的结果随机性和均匀性都比较好了。是否真的如此呢？我对于这个方法进行了实现，为了检测随机性，这里采用wx-csy/Randomness-Test的随机性测试工具，这个工具的实现遵从了国家密码管理局于2009年2月发布的《随机性检测规范》。不过该工具代码需要进行一番修改才能正确地在非Windows平台上运行，且需要编写一些代码，将二进制数据转换成一个bool的std::vector才行。这都不难。
我随意截取了 1 MiB 二进制数据进行随机性检测，结果如下：
T E S T R E S U L T
Name P-value
-----------------------------------------------------------------
Monobit frequency test 14.9934% P
Frequency test within a block 28.8003% P
Poker test 100.0000% P
Poker test 81.9964% P
Serial test (1) 17.1196% P
Serial test (2) 22.7412% P
Serial test (1) 73.9796% P
Serial test (2) 62.1269% P
Runs test 22.7953% P
Runs distribution test 65.6587% P
Test for the longest run of ones in a block 8.4909% P
Binary derivative test 69.5659% P
Binary derivative test 34.1147% P
Autocorrelation test 22.7545% P
Autocorrelation test 38.2001% P
Autocorrelation test 67.9145% P
Autocorrelation test 91.4763% P
Binary matrix rank test 3.3140% P
Cumulative test 27.1074% P
Approximate entropy test 34.9364% P
Approximate entropy test 66.6382% P
Linear complexity test 87.3137% P
Maurer's universal test 84.1940% P
Discrete Fourier Transform test 47.7300% P可以看到，所有的测试全部通过了。
基本结构
技术验证通过后，进入正式的编码环节。总的来说这个库并不复杂，只要能提供一个随机引擎，能让用户读出其中的字节流，或者提供一下转换成随机数的封装就好了。整个库由4个部分组成：engine引擎、fifobuf缓冲队列（用于暂存一部分数据，避免反复IO）、whitener均匀处理（ChaCha20实现），和封装过的各种分布。
前三个部分在开始时采用的面向对象思路实现，事实证明这种想法是有偏差的。实现完毕之后进行了一次大重构——显然一个SDR设备不可以对应多个引擎（不能同时打开多次），那构造一个单例，然后把指针传来传去的有什么意思呢？
各种分布
均匀分布
考虑到本身“白化”过的数据已经比较均匀，假如要得到一个int32_t的正整数，我们可以直接取4个字节，然后把最高位置0。这样得到的是在0至$2^{31}-1$上均匀分布的随机数。在cppreference.com上提供了一个分布可视化的思路，按照这个思路将均匀分布算法实际得到的随机数可视化一下得到的结果如下（取100,000个随机数）：
1 *************************************************
2 **************************************************
3 *************************************************
4 **************************************************
5 *************************************************
6 *************************************************
7 *************************************************
8 *************************************************
9 **************************************************
10 *************************************************结果还行。
怎么得到在$[0,1]$上分布的小数呢？嗯，得到一个整数然后除以它的范围。
伯努利分布类
在C++的<random>中，伯努利分布、二项分布、负二项分布、几何分布都被归为伯努利分布一类。这几种分布的生成，特点就是简单，可以说是直接模拟。只要我们能够得到均匀分布的小数，做到都不难，就不详细说了。
正态分布类
正态分布
这里采用了Box-Muller变换法。取100,000个数，均值0，标准差5，可视化效果如下：
-13 *
-12 **
-11 ***
-10 *****
-9 ********
-8 **********
-7 ***************
-6 *******************
-5 ************************
-4 *****************************
-3 *********************************
-2 ************************************
-1 ***************************************
0 ****************************************
1 ***************************************
2 ************************************
3 *********************************
4 ****************************
5 ************************
6 *******************
7 **************
8 ***********
9 ********
10 *****
11 ***
12 **
13 *对数正态分布
有了正态分布，要获得对数正态分布就简单了，可以直接变换。取一$r$服从$N(0,1)$，那么参数$(\mu, \sigma)$的对数正态分布为：
$$\exp(r\ln(\sigma) + \ln(\mu)) $$
需要注意的是$r\sim N(0,1)$，$r$不是一定落在$[-1,1]$上的。需要取一个$[-1,1]$上的$r$才行。
泊松分布类
泊松分布
可将$[0,1]$上的均匀分布转换到泊松分布：
int32_t radiom_poisson(double mean)
{
int32_t r;
double limit = exp(-mean);
double p = radiom_uniform_double();
for(r = 0; p >= limit; ++r)
p *= radiom_uniform_double();
return r;
}指数分布
取一$r\sim U(0,1)$，则：
$$\displaystyle\frac{\ln(1-r)}{-\lambda} $$
威布尔分布
取一$r\sim U(0,1)$，则有比例参数$\lambda$、形状参数$k$的威布尔分布：
$$\lambda\cdot(-\ln(r))^{k^{-1}} $$
基本用法
安装没什么好说的，就是一个普通的CMake工程。你可以这样做：
git clone https://github.com/Bokjan/Radiom
cd Radiom
mkdir build
cd build
cmake ..
make
sudo make install
sudo ldconfig -v | grep radiom举例
下面的例子中，我们将硬件频率设置到 457.7 MHz，并且输出10个在$[1, 10]$均匀分布的整数。
#include <stdio.h>
#include <radiom.h>
int main(int argc, char *argv[])
{
int cnt = 10;
radiom_init();
radiom_engine_set_freq(457700000); /* Hz */
radiom_engine_start();
while(cnt-- > 0)
printf("%d\n", radiom_uniform_int32_in(1, 10));
radiom_exit();
}在生成随机数之前，需要调用radiom_init、radiom_engine_start来初始化、启动引擎。若要自定义频率，则应该在radiom_engine_start之前调用radiom_engine_set_freq，参数是频率，单位是赫兹。使用者只需要显式引入radiom.h头文件即可。
API 参考
引擎部分
void radiom_init(void); 初始化Radiom环境，应当在任何其他操作前执行一次。
void radiom_exit(void); 退出Radiom环境。
int radiom_engine_start(void); 启动引擎。返回RADIOM_ERROR为失败，RADIOM_SUCCESS为成功。
int radiom_engine_stop(void); 停止引擎。调用radiom_exit时会自动停止。
int radiom_engine_getbytes(void *dest, size_t count); 将count个字节的随机字节存放到dest处，返回值意义同上。
int radiom_engine_entropy_bytes(void *dest, size_t count); 将count个字节的原始无线电数据存放到dest处，返回值意义同上。
int radiom_engine_refresh_buffer(void); 强制刷新已缓存的无线电数据，返回值意义同上。
void radiom_engine_set_freq(int frequency); 设置SDR硬件频率，单位是赫兹，应在引擎启动前进行。
概率分布部分
均匀分布
int32_t radiom_uniform_int32(void);
int32_t radiom_uniform_int32_in(int32_t min, int32_t max);
int64_t radiom_uniform_int64(void);
int64_t radiom_uniform_int64_in(int64_t min, int64_t max);
double radiom_uniform_double(void);提供5种均匀分布接口。带参数的整数类型是左闭右闭区间，不带参数的整数类型是$[0,\text{MAX}]$上的非负数。浮点型的范围是$[0,1]$。
伯努利分布
int radiom_bernoulli(double probability);
int32_t radiom_binomial(int32_t n, double probability);
int32_t radiom_negative_binomial(int32_t n, double probability);
int32_t radiom_geometric(double probability);共提供4种，即伯努利分布、二项分布、负二项分布、几何分布，参数意义很明显，不详细描述。其中radiom_bernoulli的返回值只会是0或1。
正态分布
double radiom_normal(double mean, double stddev);
double radiom_lognormal(double mean, double stddev);提供正态分布和对数正态分布，生成方法已在上面描述。参数是均值和标准差（standard deviation）。
泊松分布
int32_t radiom_poisson(double mean);
double radiom_exponential(double mean);
double radiom_weibull(double shape, double scale);提供泊松分布、指数分布和威布尔分布，生成方法已在上面描述。泊松分布和指数分布的参数均是均值$\lambda$；威布尔分布有两个参数：形状（shape）参数$k$、比例（scale）参数$\lambda$。
用途
本有做类似于RANDOM.ORG的服务的想法。其实这件事情是很麻烦的，刨去其他方面不说，就受硬件设备限制，不可能跑在公有云上；做内网穿透呢，又麻烦又不稳，只好作罢。事实上我试着结合Scaffold做了一个随机双色球投注的JSON接口：
#include <set>
#include <cstdio>
#include <cstdlib>
#include <radiom.h>
#include <scaffold/Request.hpp>
#include <scaffold/Response.hpp>
#include <scaffold/Scaffold.hpp>
#include <rapidjson/writer.h>
#include <rapidjson/document.h>
#include <rapidjson/stringbuffer.h>
void InitializeRadiom(void)
{
radiom_init();
radiom_engine_set_freq(1210500000); // 121.05 MHz, Guangzhou Center (ATC)
if(radiom_engine_start() != RADIOM_SUCCESS)
{
fputs("ENGINE START FAILED!", stderr);
exit(-1);
}
}
void ShuangSeQiu(Request &req, Response &res)
{
using namespace rapidjson;
int count = 1;
if(req.query.find("count") != req.query.end())
sscanf(req.query["count"].c_str(), "%d", &count);
if(count > 100 || count < 1)
count = 1;
StringBuffer s;
Writer<StringBuffer> w(s);
w.StartArray();
while(count--)
{
std::set<int> s;
w.StartObject();
w.Key("red");
while(s.size() < 6)
s.insert(radiom_uniform_int32_in(1, 33));
w.StartArray();
for(const auto i : s)
w.Uint(i);
w.EndArray();
w.Key("blue");
w.Uint(radiom_uniform_int32_in(1, 16));
w.EndObject();
}
w.EndArray();
res.type("json");
res.send(s.GetString());
}
int main(void)
{
InitializeRadiom();
auto app = Scaffold();
app.get("/lottery/shuangseqiu", ShuangSeQiu);
app.listen(8000);
radiom_exit();
}Anyway, it works.
话说回来，事实上目前还没有什么实质性的应用；也由于时间问题，像样的文档仍然在路上。写下此文，至少对几百行代码有个交代。
