世界的本质是旋转(9) 利用声波传输mFSK多音调频文本

在上一篇文章中，主要介绍的是使用专用的业余无线电设备传输调相波形。

对于买不起SDR设备的学生来说，可以使用这篇文章介绍的思路，使用声卡的数据线传输IQ路的基带数据。线路输入的好处：

不经过空气的媒介，波形的本质一直保持在电磁波，噪声极少，相位线性度很高，没有什么失真。
左右声道独立传输，同步采集，恰好可以用来表示IQ路的波形。
线路输入的采样率通常可以很高，有的声卡支持到192000，也就是192KHz。

但是并不是所有的计算机都有线路输入口。目前大量的微型主板、笔记本电脑采用的是四芯耳麦口，如果强行把麦的连线和对方耳机通道贯通，阻抗也不匹配。那如果采用喇叭直放呢？笔者本来以为只是做一个实信号到0频段复信号的搬移、滤波即可，但一通实验下来，非常抓狂，发现一旦电磁波换能为机械波，其信道特性就变了：

十几块钱的麦克风、粗制滥造的USB小喇叭的频率、相位响应千奇百怪。
一般麦克风和音箱即使能够字面上支持极高的采样率，比如192K，但其实会发现大部分产品对12KHz以上的频率衰减就很大了，几乎不能用。
操作系统里设置的均衡器要清空，否则各种“大厅”、“低音炮”模式会加重失真。

考虑到上述原因，通过麦克风接收到的声波的相位、幅度已经不太靠谱了，我们这次使用更为直接的多音调频的方式传输数据。即使在这种模式下，也是更换了3种麦克风，最终用一个比较贵的麦克风（60多元）配合 Lenovo 笔记本自带的扬声器实现了测试。10元的扬声器和麦克风收到的音频频偏很大，且回音、黏连干扰很严重，背景还有很大的干扰（类似一种嗡嗡的声音）。

1. 实验目的

生成一个多音波形，用做基于麦克风、小音箱的点对点数据传递。

2. 波形设计

协议栈沿用上几篇的BPSK连续收发代码，只对调制进行设计如下：

层级	协议
网络层	字符或者IP
链路层	pseudo HDLC
物理层-纠错	连续VIT序列
物理层-波形	mFSK

2.1 mFSK多音位置选取

在归一化采样率（声卡是实信号）下，有效频带为 $[0,\pi]$ ， $[\pi,2\pi]$ 是实信号的镜像频率。我们利用其中一半的资源，做 $2^m$ 分频,其中m范围为 [1,4]，对应2,4,8,16fsk。以4 fsk为例，m=2，其归一化频率位置如下：

mFSK

在以采样率归一化的频率下，每个实样点相位前进 $\pi$ 就是采样窗口能够表达的最高频率，比如采样率为24000Hz，对应 0~ $2\pi$ , 折半后为12000Hz，即12KHz。
在生成样点时，控制相位值每个样点前进0.25-0.75 $\pi$ 之间的某个角度，即生成了上图的各个频率位置。

相位推进的代码如下：

	//symbols 里是fsk符号const double c_pi = 3.14159265355;static double phy = 0; //当前样点相位const int msk  = 2;  //2^2 = 4fskconst double freq_span =  0.5 * c_pi /((1<<msk)-1);const int N = 128;//可根据调制指数反推for (int sym :symbols){//sym from 0 to 3, 归一化频率是 nfdouble nf = c_pi * 0.25 + sym * freq_span;//声卡单声道实波形,每个符号N个样点for (int i=0;i<N;++i){//加窗抑制旁瓣double win = sin(c_pi * i/(N-1)) * exp(-4.0*(i-N/2)/N*(i-N/2)/N);//输出样点int v = 16384 * cos(phy) * win;append_wav(v);//推进相位phy += nf ;}}}return signal;

可以注意到的是：

每个符号持续的样点数可以设置，这里是128个点。样点数越长，速率越慢，但频谱越趋近梳状。比如极端情况样点设置为24000点，那在24K采样率下，1秒才传输1个符号。样点数越少，传输越快，但切换造成的频谱展宽更宽，直到完全重叠。描述这种频谱利用率的参数叫做调制指数。由于我们先确定了频率的位置、间隔，所以可以选取一个调制指数，计算N的值。

bad
但由于低成本扬声器、话筒本身的恶劣信道特性，带来的非线性失真太严重，不能用电磁波的fsk来对待声波。这个值在 16fsk下，取到96已经是极限了，非常挑剔扬声器和话筒。此时，若以24K采样率，传输16fsk,其速率为：

$\delta=24000/96 Bd= 250Bd=1000bps$

如果希望双工通信，则进行频率交叉或者上下分片，速率继续折半。
加窗是为了在符号交叠区域幅度归零，抑制各个符号之间的频谱旁瓣干扰。

2.2 mFSK 映射取值

为保证相邻符号只变化1bit，减少误码率，我们可以使用简单的算法构造一个距离矩阵，找出以二进制相似性为约束的映射取值：

void init_mapCodes(int map_codes[],int map_syms[])
{bool bChecked[256]={true};map_codes[0]= 0;map_syms[0] = 0;for (int i=1;i<256;++i)bChecked[i] = false;fprintf(stderr,"Init MAP : 0 ");for (int g=1;g<(1<<msk);++g){std::map<int,std::set<int> > map_dis;for (int i=1;i<(1<<msk);++i){int h = i ^ map_codes[g-1];if (!bChecked[i])map_dis[qPopulationCount((quint32)h)].insert(i);}if (!map_dis.size())break;std::set<int> set_good = map_dis.begin()->second;int test = g-2;if (set_good.size()>1 && test >=0){std::map<int,std::set<int> > map_test;for (int p:set_good){int v = 0;for (int t = 0; t<= test;++t){int h = p ^ map_codes[test - t];v += qPopulationCount((quint32)h);}map_test[v].insert(p);}set_good = map_test.begin()->second;}bChecked[*set_good.begin()] = true;map_codes[g] = *set_good.begin();map_syms[*set_good.begin()] = g;fprintf(stderr,"%X ",map_codes[g]);}fprintf(stderr,"\nMAP_CODES=");for (int i=0 ;i<(1<<msk) ;++i){fprintf(stderr,"%02X ",map_codes[i]);}fprintf(stderr,"\nMAP_SYMS=");for (int i=0 ;i<(1<<msk) ;++i){fprintf(stderr,"%02X ",map_syms[i]);}fprintf(stderr,"\n");for (int i=0 ;i<(1<<msk) ;++i){fprintf(stderr,"%02X ",map_codes[i]);for (int j=0;j<(1<<msk) ;++j){int h = map_codes[i] ^ map_codes[j];int v = qPopulationCount((quint32)h);fprintf(stderr,"%2d ",v);}fprintf(stderr,"\n");}
}

以 8fsk\16fsk为例，各个符号之间的距离矩阵如下：

   00 01 03 02 06 04 05 07 
00  0  1  2  1  2  1  2  3 
01  1  0  1  2  3  2  1  2 
03  2  1  0  1  2  3  2  1 
02  1  2  1  0  1  2  3  2 
06  2  3  2  1  0  1  2  1 
04  1  2  3  2  1  0  1  2 
05  2  1  2  3  2  1  0  1 
07  3  2  1  2  1  2  1  0

   00 01 03 02 06 04 05 07 0F 0B 09 0D 0C 08 0A 0E 
00  0  1  2  1  2  1  2  3  4  3  2  3  2  1  2  3 
01  1  0  1  2  3  2  1  2  3  2  1  2  3  2  3  4 
03  2  1  0  1  2  3  2  1  2  1  2  3  4  3  2  3 
02  1  2  1  0  1  2  3  2  3  2  3  4  3  2  1  2 
06  2  3  2  1  0  1  2  1  2  3  4  3  2  3  2  1 
04  1  2  3  2  1  0  1  2  3  4  3  2  1  2  3  2 
05  2  1  2  3  2  1  0  1  2  3  2  1  2  3  4  3 
07  3  2  1  2  1  2  1  0  1  2  3  2  3  4  3  2 
0F  4  3  2  3  2  3  2  1  0  1  2  1  2  3  2  1 
0B  3  2  1  2  3  4  3  2  1  0  1  2  3  2  1  2 
09  2  1  2  3  4  3  2  3  2  1  0  1  2  1  2  3 
0D  3  2  3  4  3  2  1  2  1  2  1  0  1  2  3  2 
0C  2  3  4  3  2  1  2  3  2  3  2  1  0  1  2  1 
08  1  2  3  2  3  2  3  4  3  2  1  2  1  0  1  2 
0A  2  3  2  1  2  3  4  3  2  1  2  3  2  1  0  1 
0E  3  4  3  2  1  2  3  2  1  2  3  2  1  2  1  0

上述算法保证相邻的符号只差1比特，次相邻的差2比特。这样即使左右错位，也只差1比特，用VIT比较容易纠错。

3. 数据的接收

对于电磁波实信号处理，要测量其相位，可变到0频（复数信号），直接计算角度变化，就是频率。这种计算很方便，参考《50行实现C语言FM收音机-Taskbus Stdio封装器在SDR课程中的应用》一文。

但是对声波，由于扬声器、麦克风本身的特点，我们准备采用另一种更为鲁棒的方式，即短时FFT变换来进行接收、解调。

3.1 时域波形

使用Audacity可以从声卡录制，通过时域波形可以看到明显的频率选择。对某些低频的增益大，高频部分不仅更弱，相位也有抖动。

这种情况下，可以通过两步来完成跟踪解调。

同步红色窗口的位置，使其对准每个符号。
判决符号，生成比特流。

3.2 短时FFT偏移同步

使用N=符号长度的窗口截取波形，使用短时fft可以对符号窗口呢的数据进行变换，获得瞬时的频谱。

如果窗口没有对准符号位置，则FFT后，相邻两个符号的位置都有能量。

如果窗口对准了符号的位置，则只有当前符号的频谱能量最高。

如果已经测量当前时刻各个fsk频点位置上的能量m，并按照大小排序，最强的放在第一个，次强的第二个，以此类推，则可以定义一个测度：

$\eta=\cfrac{m_1}{m_2}$

作为跟踪最佳取样位置的依据。下图显示的是符号变换的过程，在最佳位置上，只有1个显著的峰值，在过渡位置，有多个。

在这里插入图片描述
下面的代码通过fft峰值跟踪，锁定最佳的判决窗口偏移。

const double c_pi = 3.14159265355;
int msk  = 3;
//double * in;
fftw_complex *out = 0, *in = 0;
fftw_plan p = 0;
int fftsize = 512;
int sym_len = 128;
int fft_win = fftsize / 2;
int sym_order = 1;
int sym_freqct = 2;
int pos_ffttest[256];std::vector<unsigned char> demod_msk(const vector<unsigned char> &packagedta)
{static unsigned long long clk = 0;static int map_codes[256]={-1},map_syms[256]={-1};static short last_syms[1024];static int offset = 0;static double max_curr_v = 0;static double max_left_v = 0;static double max_right_v = 0;static long long nextJudge = 0;static double fftv[4096];if (map_codes[0]<0){init_mapCodes(map_codes,map_syms);}std::vector<unsigned char> results;//声卡单声道实波形int nSyms = packagedta.size()/2;short * pSym = (short *) packagedta.data();bool fftout = false;for (int i=0;i<nSyms;++i){last_syms[clk % sym_len] = pSym[i];++clk;if (clk <sym_len*160+33)continue;//Judgeif ((clk - offset) % sym_len==0){max_curr_v = 0;int max_curr_s = 0;double secd = 0;//FFTfor (int j = 0; j<fftsize;++j){if (j<sym_len)in[j][0] = last_syms[(clk - sym_len -1 + j) % sym_len];elsein[j][0] = 0;in[j][1] = 0;}fftw_execute(p);for (int j=0;j<sym_freqct;++j){double absv = sqrt(out[pos_ffttest[j]][0]*out[pos_ffttest[j]][0] + out[pos_ffttest[j]][1]*out[pos_ffttest[j]][1]);if (max_curr_v < absv){secd = max_curr_v;max_curr_v = absv;max_curr_s = j;}else if (secd < absv)secd = absv;}if (sub_fft && !fftout ){for (int j=0;j<fft_win;++j){double absv = sqrt(out[j][0]*out[j][0] + out[j][1]*out[j][1]);fftv[j] = absv;}TASKBUS::push_subject(sub_fft,instance,(unsigned int)fft_win*sizeof(double),(unsigned char *)fftv);fftout = true;}int syb = map_codes[max_curr_s];for (int j=0;j<sym_order;++j)results.push_back((syb>>(sym_order-1-j))&0x01);max_curr_v = max_curr_v / (secd+0.0001);}else if ((clk - offset - 4) % sym_len==0){max_left_v = 0;double secd = 0;//FFTfor (int j = 0; j<fftsize;++j){if (j<sym_len)in[j][0] = last_syms[(clk - sym_len -1 + j) % sym_len];elsein[j][0] = 0;in[j][1] = 0;}fftw_execute(p);for (int j=0;j<sym_freqct;++j){double absv = sqrt(out[pos_ffttest[j]][0]*out[pos_ffttest[j]][0] + out[pos_ffttest[j]][1]*out[pos_ffttest[j]][1]);if (max_left_v < absv){secd = max_left_v;max_left_v = absv;}else  if (secd < absv)secd = absv;}max_left_v = max_left_v / (secd+0.0001);}else if ((clk - offset + 4) % sym_len==0){max_right_v = 0;double secd = 0;//FFTfor (int j = 0; j<fftsize;++j){if (j<sym_len)in[j][0] = last_syms[(clk - sym_len -1 + j) % sym_len];elsein[j][0] = 0;in[j][1] = 0;}fftw_execute(p);for (int j=0;j<sym_freqct;++j){double absv = sqrt(out[pos_ffttest[j]][0]*out[pos_ffttest[j]][0] + out[pos_ffttest[j]][1]*out[pos_ffttest[j]][1]);if (max_right_v < absv){secd = max_right_v;max_right_v = absv;}else  if (secd < absv)secd = absv;}		if (nextJudge ==0)nextJudge = clk + sym_len +15;max_right_v = max_right_v / (secd+0.0001);}//钟差跟踪if (clk > nextJudge && nextJudge >0){//adjust offset;if (max_curr_v < max_right_v)--offset;else if (max_curr_v < max_left_v)++offset;nextJudge = 0;}}static int cct = 0;if (++cct % 50==0)fprintf(stderr,"clk=%lld,offset=%d\n",clk,offset);return results;
}