1 引言
语音信号传输的应用广泛，以语音信号为载体实现隐蔽通信具有不易被察觉的特点，在信息安全传输领域具有重要应用价值[1-2]。基于掩蔽阈值分析的隐藏算法能够跟随语音载体非平稳特征自适应调整嵌入位置和强度，在满足信息嵌入后信噪比要求的前提下，能够使信号在听觉意义上失真降为最小，具有较好的不可感知性。文献[3]通过小波包分解和离散余弦变换(discrete cosine transform，DCT)估算掩蔽阈值，自适应选取信息嵌入位置和嵌入强度，与人耳频域掩蔽特性关联较好，但临界频带划分不够准确，掩蔽阈值计算复杂，信息提取时的嵌入强度盲检测存在较大误差，算法自身不稳定；文献[4]采用DCT变换计算掩蔽阈值，确定嵌入强度，设置能量补偿点，从而实现了信息隐藏前后能量的一致，算法稳定且鲁棒性强，但临界带划分和掩蔽阈值计算不够准确，固定分帧与载体非平稳特征关联不好，不可感知性较弱；文献[5]在离散小波变换(discrete wavelet transform，DWT)域计算掩蔽阈值，基于先前的载密数据计算当前的嵌入强度，嵌入强度盲检测准确，鲁棒性强，但嵌入强度计算复杂，掩蔽阈值计算不够准确，与载体特征关联性有待提升。
本文提出一种采用动态分帧和基于掩蔽阈值分析的自适应嵌入方法。仿真实验及性能分析表明：本文方法在保证隐藏容量达到一定要求的同时，具有运算量小、鲁棒性较强、不可感知性较好的特点，可以稳定可靠地应用于以语音为载体的隐蔽通信系统中。
2 基本原理
人耳存在掩蔽效应，充分利用这一听觉感知特性，可以使嵌入信息不可感知[6-8]。本文方法基于心理声学模型Ⅱ，采用DWT和DCT变换精确计算掩蔽阈值，自适应确定嵌入强度。由于鲁棒性与量化索引调制(quantization index modulation，QIM)嵌入强度大小直接相关[9-10]，本文采用动态分帧，将信息隐藏在载体语音能零比较大帧的低频段，在保证不可感知性前提下最大限度提高嵌入强度，并有效保证鲁棒性。盲提取要求稳定可靠实现嵌入强度盲检测，本文采用能量补偿方法，使信息嵌入前后各临界频带能量保持不变，消除信息嵌入过程对嵌入强度盲检测的影响。
2.1 动态分帧与快速重定位
由人耳掩蔽特性和鲁棒性要求可知，能零比较大的低频段更适于隐藏信息。然而语音信号是非平稳随机信号，固定分帧嵌入信息的方法与语音信号的非平稳随机特性关联不好[11-12]，本文采用动态分帧选取合适的信息嵌入帧。
语音信号能零比能够检测语音端点，区分有声段和无声段，但对数据平移不敏感，提取端难以再次准确定位隐藏信息[12]。本文根据载体音频能量比分布情况，将载体音频区分为能零比较大的数据隐藏段(有声段)和能零比较小的非隐藏段(无声段)，秘密信息隐藏在载体音频能零比较大的位置。图1为动态分帧结构示意图，数据隐藏段由起始帧、信息嵌入帧和结束帧构成。其中，对有声段语音做小波变换后，采用QIM法在小波系数中插入m序列标识起始帧。结束帧通过对3层尺度系数进行能零比检测来确定，用于标识当前数据隐藏段的终端，并做容错处理防止误检测。在起始帧和结束帧之间的信息嵌入，以1帧为单位进行逐帧嵌入秘密信息。
图1 动态分帧结构示意图
Fig.1 The dynamic framing structure
提取端利用小波变换的时频局部特性，在小波域快速搜索匹配m序列实现精准重定位。Haar小波变换[13]3层系数只与原信号8个样点有关，当原信号平移8个样点时，3层系数产生1个样点的平移。在3层小波系数中搜索同步码时，可利用这种特殊的平移性，只全部计算前8次逐点平移的小波系数，后续小波系数可依次通过对应的前8次小波系数向左平移1位，再计算最后8个样点对应的小波系数做为最后1位而得到，这样可以极大提高搜索效率。
m序列自相关函数为
(1)
式(1)中，p为m序列的周期。设嵌入的m序列为m(i)，由QIM法检测到的m序列为m′(i)，对m(i)和m′(i)做自相关检测，如果有
(2)
就认为检测到同步码。其中，z是一个适当的阈值参数，z的取值因根据信道质量、同步码嵌入强度确定。仿真实验表明，门限值一般取z/p=0.9比较合适。取值太高，误码会导致漏检测。取值太低，小波三层系数的稳定性会导致误检测，误同步。
2.2 基于心理声学模型Ⅱ的掩蔽阈值分析
基于心理声学模型Ⅱ计算掩蔽阈值更为精确，当信息嵌入引入噪声能量谱密度低于掩蔽阈值一定值时，就能保证隐藏信息不被感知。语音信号能量集中在低频段，为提高算法鲁棒性，本文只将信息隐藏在载体语音的低频段，采用DWT和DCT变换在载体低频段划分临界频带[14-15]，计算各临界频带的掩蔽阈值。
2.2.1 临界频带划分
人耳听觉范围为0.020～22.05 kHz，根据听觉特性可以划分为25个临界频带，临界频带序号称为临界频带域或Bark域，表1给出前9个临界频带频率。
表1 前9个临界频带频率
Table1 The first nine critical frequency band
Bark域频率/Hz中心低端高端带宽1 5020 10080215010020010032502003001004350300400100545040051011065705106301207700630770140884077092015091 0009201 080160
对于采样率fs=8 kHz的语音信号，能量主要集中在低频段，且人耳对低频信号不敏感，因此只选取1～9临界频带隐藏信息，忽略其他临界频带[9]。动态分帧时帧长为L，每帧信号先做3层哈尔(Haar)小波变换，3层尺度系数c3(i)和3层小波系数d3(i)对应频带为0～500 Hz和500～1 000 Hz；再分别对c3(i)和d3(i)做DCT变化，各得到L/8个DCT系数Xc(k)和Xd(k)，用DCT系数等间隔平分低频带宽，其频率分辨率均为fs/2L(Hz)；再将Xc(k)和Xd(k)所代表的频点对应划分到各Bark域，实现临界频带划分。本文取L=504，则频率分辨率约为7.94 Hz，对照表1中各临界频带的上下边界，可得各临界频带包含的DCT系数范围和数量如表2所示。
表2 各临界频带包含的DCT系数范围和数量
Table 2 The DCT coefficients range and quantity
for each critical frequency band
Bark域DCT系数范围数量1Xc(1)～Xc(12)122Xc(13)～Xc(26)143Xc(27)～Xc(38)124Xc(39)～Xc(51)135Xc(52)～Xc(63)Xd(1)～Xd(2)146Xd(3)～Xd(17)157Xd(18)～Xd(35)188Xd(36)～Xd(53)189Xd(54)～Xd(63)10
2.2.2 掩蔽阈值计算
本文基于心理声学模型Ⅱ计算掩蔽阈值，计算方法如下：
步骤1：当前帧内所有临界频带对第j个临界频带叠加产生的掩蔽效应由各临界频带能量与掩蔽扩展函数卷积进行定量估算。
① 计算各临界频带能量Bj。设Xj(i)为Xc(k)和Xd(k)映射到第j个临界频带的系数，则第j个临界频带的能量谱密度为：
(3)
第j个临界频带的能量为
(4)
式中， jl和jh分别为第j个临界频带内最低频点和最高频点对应的DCT系数。
② 计算掩蔽扩展函数[5]值SFj(dB)。
(5)
③ 计算扩展后的掩蔽阈值Cj。
(6)
式(6)中，*表示做卷积运算。
步骤2：掩蔽阈值的大小与音频信号的谱平坦度和音调特性有关，纯音掩蔽噪声和噪声掩蔽纯音的掩蔽阈值不同。实际信号往往既不是纯音也不是噪音，因此采用音调系数对掩蔽阈值进行修正，通过估算第j个临界频带的谱平坦度和音调系数对掩蔽阈值进行定量修正。
① 计算临界频带的谱平坦度SFM(dB)：
SFM=10lg(G/A)
(7)
式(7)中，G为临界频带能量谱密度的几何平均值，A为临界频带能量谱密度的算数平均值。
② 计算音调系数α。
α=min(SFM/-60，1)
(8)
式中，当音频信号为纯音时α=1，当音频信号为白噪声时α=0，而实际音频信号的音调系数往往介于[0，1]之间。
③ 计算掩蔽阈值Tzj(dB)。
由音调系数得到第j个临界频带掩蔽阈值修正值[1,5] 为
Oj=α(14.5+j)+(1-α)5.5
(9)
修正后的临界频带j的掩蔽阈值为
Tj=10lgCj-Oj/10
(10)
考虑估算得到的掩蔽阈值修正值Tj不应小于该临界频带的绝对掩蔽阈值Tq( f )，临界频带j的最终掩蔽阈值为：
Tzj=max(Tj，Tq(f))
(11)
式中，绝对掩蔽阈值Tq( f )(dBSPL)取值如下[5]：
(12)
由于Tj表示能量值，而Tq( f )表示声强级，二者单位不统一，因此不能直接比较二者大小。由声强级和绝对掩蔽阈值定义可知，当确定0 dB声强级对应的绝对能量谱密度值后，就可以确定各频点处的绝对能量谱密度值。本文设定DCT域0 dB声强级对应的绝对能量谱密度值为S，为便于运算，取各临界频带中点处绝对掩蔽阈值为该临界频带绝对掩蔽阈值Tq( j )，将前9个Bark域中心频率代入计算得Tq(j)，对应的绝对能量谱密度值为S( j )。如表3所示，其中S值通过实验设定，最终与仿值归一化数值实现统一。
表3 绝对掩蔽阈值
Table 3 The absolute masking threshold
Bark域中心频率/Hz中心频率绝对掩蔽阈值Tq( j )(SPL/dB)绝对能量谱密度值S( j )15039.976 19 945.1S215016.587 745.579 5S325011.009 912.618 0S43508.395 36.910 8S54506.845 34.836 5S65705.632 73.658 2S77004.729 72.971 5S88404.013 12.519 5S91 0003.369 12.172 3S
临界频带j的最终掩蔽阈值为
Tz(j)=max[Tj，lj×Tq(j)]=max[Tj，lj×S(j)]
(13)
式(13)中，lj为第j个临界频带DCT系数数量。
2.3 能量补偿
由掩蔽阈值的计算过程可知，扩展后的掩蔽阈值Cj和谱平坦度SFM都与载体信号能量谱密度相关，而信息嵌入势必影响载体信号的能量谱密度，这会导致提取端掩蔽阈值与发送端掩蔽阈值产生偏差，进而导致盲检测嵌入强度产生误差，盲提取误码率增大。
为提高盲提取的准确性和算法稳定性，本文采用能量补偿法，基本原理为：将各临界频带DCT系数分为2组，设为g1和g2。g1用于隐藏信息，g2用于能量补偿，通过调整g1组嵌入算法和g2组DCT系数大小，使得各临界带信息嵌入前后能量不变，即：
(14)
式中，Xj(i)′为Xj(i)嵌入和调整后的数值。
令
(15)
则有
(16)
为便于计算，对g2组补偿系数做等比例缩放，则有
(17)
由上式可知，应有k≥0且k尽量接近1。为达到这一目的，首先应合理划分g1和g2组DCT系数数量，本文方法中当临界带DCT系数少于15个时，g2组选取2个DCT系数；否则，g2组选取3个DCT系数。其次，由k的表达式可知，应使g1组嵌入前后能量差最小，g2组能量和不能过小，以防止调整方位过大。为使能量差最小，量化过程应尽量满足：
(18)
为实现上述目的，本文在QIM法嵌入信息过程中采取以下措施：
① 当0≤k≤2时，不做处理；
② 当k>2时，表明嵌入后系数平方和过小，将嵌入后变小的系数逐个增大一个嵌入强度Δj，直至0≤k≤2；
③ 当k<2时，表明嵌入后系数平方和过大，将嵌入后变大的系数逐个减小一个嵌入强度值Δj直到刚好k>0为止。
同时，g2组系数最小应与初始嵌入强度
同一量级，不能过小，一是防止k值对调整过于敏感，二是防止g1组系数幅值全部小于嵌入强度需调整而无法调整的情况发生。
本文方法中对g2组系数做如下处理：
(19)
能量补偿能够消除信息嵌入对各临界频带能量Bj、算数平均值A的影响，但不能消除对几何均值G带来的影响，导致音调系数α在信息嵌入前后发生变动。但本文方法经过能零比检测处理，信息嵌入段载体信号谱平坦度变化不是十分剧烈，仿真实验也表明，掩蔽阈值修正值大体范围为[5，12]，因此，本文将Oj设定为8 dB。
2.4 嵌入强度选取
本文信息嵌入采用QIM调制的方法，量化步长的选取关系方法的不可感知性和鲁棒性，根据文献[3]的结论和实验测试，如果信息嵌入后的音频噪声掩蔽比NMR≤-5 dB，信息隐藏后具有较好的不可感知性。设QIM嵌入法量化步长或嵌入强度为Δj，则g1组单点信息嵌入引起的最大量化误差为Δj/2，进过能量补偿调整后，单点最大量化误差为Δj；g1组系数幅值越大，调整比例k越趋近于1；对于调整过的g2组系数，k值极少大于2，一般在[0，2]之间，因此可取g2组系数单点最大量化误差亦为Δj。单点最大噪声能量为
(20)
又有
NMR(j)=En-Tz(j)/lj
(21)
当量化噪声不可感知时，有
10lgNMR(j)≤-5
(22)
则当前帧第j个临界频带的嵌入强度
(23)
2.5 量化索引调制
设待隐藏的二进制信息为x(i)，采用QIM法嵌入信息的公式为
① x(i)=1时，
Xj1(i)′=
(24)
② x(i)=0时，
Xj1(i)′=
(25)
式中，Xj1(i)、Xj1(i)′为嵌入前后的DCT系数，Δj为当前帧第j个临界频带的嵌入强度，⎣·」表示向下取整。
设盲提取的二进制信息为x′(i)，采用QIM法盲提取隐藏信息的公式为
(26)
式中，Xj1(i)″为提取端当前帧各临界频带g1组DCT系数，Δj″为盲检测得到的各临界带的嵌入强度。
由QIM法嵌入公式可知，信息嵌入前后各DCT系数量化误差最大值为Δj/2。由于能量补偿算法对部分嵌入后的DCT系数进行了调整，导致部分DCT系数信息嵌入前后量化误差最大值为Δj，但调整幅度均为+Δj或-Δj，且调整前后DCT系数正负值不发生变化，故不会对提取过程造成任何影响。由QIM法提取公式可知，提取时的噪声容限为±Δj″/4，只要最终检测误差不超过噪声容限，就可以正确的恢复出隐藏信息，否则就可能产生误码；提取过程中的Δj″为盲检测得到，提取过程不需要原始载体，故能够实现盲提取。
3 信息嵌入提取过程
3.1 信息自适应嵌入
为了保证鲁棒性，本文所提方法在DCT域基于QIM在载体语音的低频段进行信息的嵌入和提取过程。在嵌入端，利用能零比检测和m序列自相关性，动态选取和标记适于信息隐藏的语音段，进而结合DWT和DCT准确划分临界频带并基于心理声学模型Ⅱ计算掩蔽阈值，实现嵌入强度的自适应调整。自适应嵌入流程如图2所示。
图2 自适应嵌入流程框图
Fig.2 The adaptive embedding flow chart
具体步骤为：
① 动态分帧，经过DWT、DCT变换计算初始掩蔽阈值，得到初始嵌入强度
② 根据预调整g2组系数；
③ 再次计算掩蔽阈值和嵌入强度Δj；
④ 采用QIM法在各临界带g1组系数中依次嵌入隐藏信息；
⑤ 计算各临界带k(j)值，依据k(j)值调整g1组相应系数量化误差；
⑥ 根据k(j)缩放g2组系数，使信息嵌入前后各临界带能量不变，实现能量补偿目的；
⑦ 经过IDCT、IDWT得到信息隐藏后数据，对应替换原数据，实现信息隐藏。
3.2 信息盲提取
在提取端，利用小波变换的时频局部特性，在小波域采用容错处理的方法实现快速准确重定位，并在各临界带采用能量补偿方法消除嵌入信息对嵌入强度盲检测准确度的影响，进而采用QIM方法提取信息。盲提取流程如图3所示。
图3 盲提取流程框图
Fig.3 The blind extraction flow chart
具体步骤为：
① 利用小波变换时频局部特性和m序列自相关检测实现快速重定位；
② 对嵌入帧做DWT和DCT变换；
③ 按照发送端方法划分临界频带，计算掩蔽阈值，确定嵌入强度Δj″；
④ 从各临界带g1组采用QIM法顺序盲提取隐藏信息；
⑤ 顺序组合得到所有隐藏信息，实现信息盲提取。
4 仿真实验及性能分析
本文基于心理声学模型Ⅱ分析掩蔽阈值，通过动态分帧、精确划分临界频带、计算掩蔽阈值，充分利用人耳掩蔽效应，自适应调整嵌入位置和嵌入强度。下面对该方法的不可感知性、隐藏容量、抗干扰能力和复杂度进行测试分析。性能测试载体语音采用采集的3段8 kHz采样率、16位采样精度、时长5 s(40 000点)单声道语音信号(样本一为女声、样本二为男中音、样本三为男低音)，隐藏信息为任意的一段随机二进制序列，同步帧和嵌入帧帧长均为504点，m序列周期p=63，0 dB声强级对应的绝对能量谱密度值S=3×10-6。
4.1 不可感知性分析
采用平均意见得分(average opinion score，MOS)和信噪比(signal to noise ratio，SNR)对不可感知性进行综合表征和评判。鉴于自适应方法的不可感知性不随载体音频改变而改变，测试对载体音频的影响不做考虑。
表4给出了测试样本的MOS和SNR值。由于嵌入信息是伪随机序列，SNR值存在小幅波动，表4中的数值为10次测量平均值。
表4 不可感知性测试
Table 4 The imperceptible experiments
MOS/分SNR/dB样本15.0约18.16样本25.0约16.65样本35.0约17.94文献[3]5.0约30文献[4]4.5约17.51文献[5]4.5约20
对表4数据进行分析可知：
① 本文方法。在保证MOS值足够大的前提下，即使SNR值较小(时域差值大，0.1量级)，不可感知性依然较好；SNR值小意味着嵌入强度大，有利于算法提高鲁棒性；算法MOS值不随载体音频不同而变化，自适应特性良好；
② 文献[3]方法。为了提高信噪比，选取嵌入系数时选取了较小的嵌入强度(时域差值0.01量级)，降低了嵌入容量，进而降低了鲁棒性，且小波包变换与临界带对应不是很准，不可感知性与本文方法相当；
③ 文献[4-5]方法。由于掩蔽阈值计算不够准确，在同样低信噪比的情况下不可感知性低于本文方法。
样本一信息嵌入前后时域如图4所示。
图4 样本一信息嵌入前后的时域
Fig.4 The time-domain comparison before and
after information embedding in sample 1
4.2 隐藏容量分析
语音信号是典型的非平稳随机信号，且不同于连续音频信号，有声段无声段特征区分明显，固定隐藏容量显然不能适应其特点，隐藏容量自适应调整具有更高应用价值。本文方法采取动态分帧选取嵌入位置，隐藏容量受载体语音自身能零比特征影响。表5给出测试样本的嵌入容量，分析可知：
表5 隐藏容量测试表
Table 5 The hiding capacity experiments
隐藏容量/bps样本11 218样本21 407样本31 449文献[3]约600文献[4]848.08文献[5]344.53
① 文献[3]自适应选取嵌入系数，嵌入容量动态变化；
② 文献[4-5]采用固定分帧和固定嵌入位置的方法，嵌入容量不随载体变化；
③ 本文方法隐藏容量与载体能零比分布情况相关，在不同测试样本下的隐藏容量不同，但都明显大于文献[3-5]之值。
4.3 抗干扰能力分析
不同于音频水印，音频信息隐藏受到的干扰主要为采样精度和计算精度引入的噪声和信道传输过程中误码导致的噪声，本文采用高斯白噪声模拟干扰，采用误比特率(bit error rate，BER)衡量抗扰能力。设置式(2)中z/p=0.9时实验数据统计性能表现较好。表6给出该条件下测试样本在不同干扰强度及不进行能量补偿时的BER值。
表6 抗干扰能力测试表
Table 6 The resistance experiments against interference
样本1样本2样本3文献[3]无扰有补偿无补偿02.82%02.43%01.79%-2.66%加噪σ=0.000 1σ=0.000 5σ=0.001σ=0.00203.14%6.4%9.98%04.36%9.18%14.97%03.4%7.38%9.86%2.81%7.16%11.57%15.15%
对表6数据进行分析可知：
① 无干扰情况下本文方法BER值始终为0、无干扰不进行能量补偿时BER值不为0，本文采用能量补偿确实提高了算法自身稳定性；
② 本文方法在受扰情况下BER值变化平缓，不存在明显的边界现象，说明本文方法有助于改善鲁棒性；
③ 文献[3]算法没有考虑信息隐藏对载体音频的影响，算法自身不稳定，且为了单纯提高信噪比而选取了较小嵌入强度，导致抗扰能力弱于本文方法；
此外，文献[4-5]采用了矢量量化方法，以降低隐藏容量换取强鲁棒性，属于强鲁棒性水印算法。在确定了具体的信息隐藏方法后，抗扰能力往往直接受嵌入强度影响，适当降低不可感知性要求，提高绝对掩蔽阈值参考值S，可以明显改善低干扰情况下的BER值。
4.4 计算复杂度分析
一个长度为N的离散数列的小波变换，使用Mallat算法的计算复杂度为O(N)；完全小波包分解的计算复杂度与DCT变换的计算复杂度相当，均为O(NlogN)；部分小波包分解的计算复杂度介于O(N)和O(NlogN)之间。
为保证频谱分析具有足够的分辨率，帧长N必须具有足够的长度。文献[3]采用部分小波包分解；文献[4]直接采用DCT变换；本文方法在动态分帧的基础上先进行3层DWT变化，使数列长度变为N/8后，再进行DCT变换。从变换过程计算量分析，文献[3]<本文方法<文献[4]算法，但当N值较小时，DCT变换和部分小波包分解运算量差别不大。因此本文方法具有和文献[4]相同的频率分辨率，优于文献[3]的频率分辨率；本文方法计算复杂度与文献[3]相当，优于文献[4]算法。
5 结论
本文方法在计算复杂度未有明显提升的情况下，掩蔽阈值计算准确，在较低的信噪比下具有较高的MOS值，不可感知性好，且不可感知性不受载体音频影响；充分利用了人耳掩蔽特性，隐藏容量较大；在信息嵌入端，采用动态分帧和自适应嵌入，提高了嵌入强度的统计平均值，在保证一定不感知性的前提下有利于抗干扰性能的提升；在信息提取端，容错处理和能量补偿措施有助于实现稳定的盲提取，保证了方法的稳定性。
[1] Sun R，Jiang X G，Yang X M.Audio watermarking algorithm based on wavelet packet and psychoacoustic model[J].Journal of System Simulation，2007，19(10)：2248-2251.
[2] Hua G，Huang J，Shi Y Q，et al.Twenty years of digital audio watermarking-a comprehensive review[J].Signal Processing，2016，128(11)：222-242.
[3] Zhang T，Zhang C X，Gao X Y.Adaptive audio watermarking algorithm in the hybrid domain[J].Journal of Signal Processing，2017，33(6)：828-835.
[4] Hu H T，Hsu L Y.Robust，transparent and high-capacity audio watermarking in DCT domain[J].Signal Processing，2015，109(C)：226-235.
[5] Hu H T，Hsu L Y.A DWT-Based rational dither modulation scheme for effective blind audio watermarking[J].Circuits，Systems，and Signal Processing，2016，35(2)：553-572.
[6] Painter T，Spanias A.Perceptual coding of digital audio[J].Proceedings of the IEEE，2000，88(4)：451-515.
[7] He X，Michael S.Scordilis.An enhanced psychoacoustic model based on the discrete wavelet packet transform[J].Journal of the Franklin Institute，2006，343(7)：738-755.
[8] Jiang W，Huang X，Quan Y.Audio watermarking algorithm against synchronization attacks using global characteristics and adaptive frame division[J].Signal Processing，2019，162(9)：153-160.
[9] Terchi Y，Bouguezel S.A blind audio watermarking technique based on a parametric quantization index modulation[J].Multimedia Tools and Applications，2018，77(19)：25681-25708.
[10] Wang L N，Zhang H G.Information hiding technology and application(the second edition)[M].Wuhan:Wuhan University Press，2009.
[11] Luo Yifan，Peng Dezhong，et al.Dual-domain audio watermarking algorithm based on flexible segmentation and adaptive embedding[J].IEEE Access，2019(1)：10533-10545.
[12] Song Z Y.Application of MATLAB in speech signal analysis and synthesis[M].Beijing:Beijing University of Aeronautics and Astronautics Press，2013.
[13] 张悦，张旭，丁一夫.Era-glonass紧急呼叫系统关键性能测试技术研究[J].重庆理工大学学报(自然科学)，2019，33(4)：41-48.
Zhang Y，Zhang X，Ding Y F.Research on key performance test technology of era-glonass emergency call system[J].Journal of Chongqing University of Technology(Natural Science)，2019，33(4)：41-48.
[14] Rekik S，Driss Guerchi.Speech steganography using wavelet and Fourier transforms[J].EURASIP Journal on Audio，Speech，and Music Processing，2012，2012(1)：20.
[15] Hu H T，Lee T T.Hybrid blind audio watermarking for proprietary protection，tamper proofing，and self-recovery[J].IEEE Access，2019，7(1)：180395-180408.