回声消除 Speex AEC 论文和代码解析

参考

"On Adjusting the Learning Rate in Frequency Domain Echo Cancellation With Double-Talk"

"Multidelay Block Frequency Domain Adaptive Filter"

"Echo Canceler with Two Echo Path Models "

注：由于包含很多LaTeX公式，可能需要多次刷新或者使用google浏览器

大纲

包括小节：1 信号模型、2 迭代公式推导、3 最优学习率、4 量的估计、5 双讲分析、6 代码解析

speexdsp AEC 的三个主要设计：

1）Multi-delay frequency domain adaptive filter，好处是：a）用来进行低延迟自适应，b）可以为每个频点单独使用学习率

2）Two path model，好处是：a）双滤波器切换，避免双讲时发散

3）Optimal learning rate under noise：a）噪声环境下能更好地控制学习率，b）避免双讲发散，对近端语音快速响应

下面对“On Adjusting the Learning Rate in Frequency Domain Echo Cancellation With Double-Talk”这篇论文（也就是上面的第三点）进行解析，并挑一些代码中的技术点进行梳理。

AEC系统图

噪声环境下最优NLMS学习率

1 信号模型

$\begin{aligned}&d(n) : 目标信号\\&\hat{y}(n): 估计的回声信号\\&\hat{w}: 估计的滤波器权重\\&w: 理想滤波器权重\\&x: 远端信号\\&[\cdot]^*: 表示求复数共轭\end{aligned}$

则，系统的误差公式如下：

$\begin{aligned}e(n) &= d(n) - \hat{y}(n) \\&= d(n) - \sum_{k=0}^{N-1}\hat{w}_k(n) x(n-k) \quad (1)\end{aligned}$

基于 NLMS 的迭代公式如下：

$\begin{aligned}\hat{w}_k(n+1) &=\hat{w}_k(n) + \mu\frac{e(n)x^*(n-k)}{\sum_{i=0}^{N-1} |x(n-i)|^2} \quad (2) \\ \end{aligned}$

代入（1）式有：

$\begin{aligned}\hat{w}_k(n+1) &=\hat{w}_k(n) + \mu\frac{[ d(n) - \sum_{k=0}^{N-1}\hat{w}_k(n) x(n-k) ]x^*(n-k)}{\sum_{i=0}^{N-1} |x(n-i)|^2} \quad (3) \\ \end{aligned}$

记 $\delta_k(n) = \hat{w}_k(n) - w_k(n)$ 为第 $n$ 次迭代后，滤波器权重第 $k$ 个系数的估计误差；并且，对于理想滤波器有： $d(n) = v(n) + \sum_{k=0}^{N-1} w_k(n) x(n-k)$ 。

将 $\delta_k(n),d(n)$ 的表达式代入（3）式中，并考虑到理想滤波器系数（声学路径）并不随着时间发生较大变化，可以认为 $w(n) = w(n+1)$ ，得到：

$\delta_k(n+1) = \delta_k(n) + \mu\frac{[v(n)-\sum_{i=0}^{N-1} \delta_i(n)x(n-i)]x^*(n-k)}{\sum_{i=0}^{N-1}|x(n-i)|^2} \quad (4)$

接着定义在第 n 次迭代时，滤波器失配度：

$\Lambda(n) = \sum_{k=0}^{N-1}\delta_k^*(n) \delta_k(n)$

2 滤波器失配度的迭代公式推导

由失配度的定义和迭代公式（4）可得，在第 n+1 次迭代时，滤波器的失配度为：

$\Lambda(n+1) = \sum_{k=0}^{N-1}\left | \delta_k(n) + \mu \frac{[v(n)-\sum_{i=0}^{N-1}\delta_i(n)x(n-i)]x^*(n-k)}{\sum_{i=0}^{N-1}|x(n-i)|^2} \right | \quad (5)$

若我们作一个强假设：远端信号 $x(n)$ 和近端语音 $v(n)$ 是互不相关的白噪声，假设给定 $x(n),\delta_i(n)$ 的条件下 $v(n)$ 的条件分布为 $p_{v|x(n),\Lambda(n)}$ ，则对 $\Lambda(n+1)$ 关于分布 $p_{v|x(n),\Lambda(n)}$ 求期望得：

$\mathbb E\{ \Lambda (n+1) | \Lambda(n), x(n) \} = \Lambda(n)\left[ 1 - \frac{2\mu}{N} + \frac{\mu^2}{N} + \frac{\mu^2\sigma_v^2}{\Lambda(n)\sum_{i=0}^{N-1}|x(n-i)|^2 }\right] \quad (6)$

其中： $\sigma_v^2$ 是近端语音信号得方差 $\sigma_v^2 = \mathbb E\{|v(n)|^2\}$ 。

我们需要对（6）式做一些推导：（为了简洁，将时间 n 略去）

先直接对绝对值的平方进行展开：

$\begin{aligned}\mathbb E\{\Lambda(n+1)|\Lambda(n), x\} &= \mathbb E\{\sum_k \left |\delta_k + \mu\frac{(v - \sum_i \delta_i x_i)x_k^*}{\sum_i|x_i|^2} \right |^2\} \\&= \mathbb E \{ \underbrace{ \sum_k |\delta_k|^2}_{(E1)} + \underbrace{\sum_k \delta_k^*\mu\frac{(v-\sum_i\delta_i x_i)x_k^*}{\sum_i |x_i|^2}}_{(E2)} \\&\quad +\underbrace{\sum_k\delta_k\mu\frac{(v^*-\sum_i\delta_i^*x_i^*)x_k}{\sum_i|x_i|^2}}_{(E3)} \\&\quad + \underbrace{\sum_k\mu^2\frac{|v-\sum_i\delta_ix_i|^2\cdot |x_k|^2}{(\sum_i|x_i|^2)^2}}_{(E4)}\}\end{aligned}$

其中： $(E1)$ 显然就是时刻 n 的滤波器失配度 $\Lambda(n)$ ， $(E2),(E3)$ 中有关 $v(n)$ 的一次项式的部分，在其取期望后为 0，即：

$\mathbb E\{ \sum_k\frac{\delta_k^* x_k^*\mu v}{\sum_i|x_i|^2}\} = 0,\quad \mathbb E\{ \sum_k\frac{\delta_k x_k\mu^* v^*}{\sum_i|x_i|^2}\} = 0$

余下的部分，重新整理得：

$\begin{aligned}\mathbb E\{ \Lambda(n+1)| \Lambda(n), x \} &= \mathbb E \{ \Lambda (n) - \underbrace{ \frac{2\mu\sum_k\delta_k^* x_k^*\sum_i\delta_ix_i}{\sum_i |x_i|^2} }_{(E5)} \\&\quad +\underbrace{\frac{\mu^2(|v|^2-v^*\sum_i\delta_ix_i-v\sum_i\delta_i^*x_i^* + \sum_i\delta_ix_i\sum_k \delta_k^* x_k^*)}{\sum_i|x_i|^2}}_{(E6)}\}\end{aligned}$

再一次，在上式 $(E6)$ 中，对 $v$ 得一次项取期望后为 0；

再对上式中 $(E5)$ 子式进行处理，其分子如下化简：

$\begin{aligned}&\mathbb E\{ \sum_k\sum_i\delta_k^*x_k^*\delta_i x_i \}\\&= \sum_k \mathbb E\{ \delta_k^*\delta_kx_k^*x_k \} + 交叉项\\&=\sum_k |\delta_k|^2\sigma_x^2 \\&= \sigma_x^2\Lambda(n)\end{aligned}$

其中 “交叉项” 为 0，因为假设（远端）语音是白噪，因此不相关； $\sigma_x^2$ 为远端语音信号得方差。

再对 $(E5)$ 中的分母改写为：

$\sum_k|x_i|^2 = N(\frac{1}{N}\sum_i |x_i|^2) \approx N\hat{\sigma}_x^2$

又 $\sigma_x^2 \approx \hat{\sigma}_x^2$ ，全部代入 $(E5)$ 子式得： $(E5) = \frac{2\mu\Lambda(n)}{N}$

同理，对子式 $(E6)$ 使用同样得处理可得：

$\mathbb E\{\Lambda(n+1)|\Lambda(n),x\} = \Lambda(n) + \frac{2\mu\Lambda(n)}{N} + \frac{\mu^2\sigma_v^2}{\sum_i|x_i|^2} + \frac{\mu^2\Lambda(n)}{N} \quad (*)$

提取公因子 $\Lambda(n)$ ，就得到了论文中的公式（6）。

3 最优学习率

接下来，我们寻找在新的一轮迭代中，使得滤波器失配度 $\Lambda(n+1)$ 达到最小的学习率；由于上式 $(*)$ 是一个关于学习率 $\mu$ 的二次函数，因而容易得到最小值点，求导并令其为 0 得：

$\frac{-2}{N} + \frac{2\mu}{N} + \frac{2\mu\sigma_v^2}{\Lambda(n)\sum_{i=0}^{N-1}|x(n-i)|^2} = 0 \quad (7)$

解方程得：

$\mu_{opt}(n) = \frac{1}{1+\frac{\sigma_v^2}{\Lambda(n)(1/N)\sum_i|x(n-i)|^2}} \quad (8)$

Observation：当没有近端语音时，即 $\sigma_v^2 = 0$ ，此时（8）式退化成 NLMS 对应的学习率。

虽然我们得到了噪声环境中的最优学习率（8）式，但是其中的失配度无法计算，因此我们需要进一步将其变为可估计的量。

3.1 得到可估计的学习率公式

继续分析（8）式中的 $\Lambda(n)\frac{1}{N}\sum_i|x(n-i)|^2$ 这一项，其中左半边 $\Lambda (n)$ 为滤波器的失配度，右半边 $\frac{1}{N}\sum_i|x(n-i)|^2$ 显然是远端信号能量的估计；因此，这一整项似乎与 AEC 滤波后的远端信号残留能量相关的某个量；下面我们证明，它的确就是远端信号的残留能量的估计，推导过程如下：

$\begin{aligned}\sigma_r^2 &= \mathbb E\{ |\sum_i \hat w_ix_i- \sum_i w_ix_i) |^2\} \\&=\mathbb E\{\sum_i | (\hat{w}_i - w_i)x_i |^2 \}\\&= \sum_i\mathbb E\{|\hat{w}_i-w_i|^2 \cdot|x_i|^2\} +\sum_{i\ne j} \mathbb E \{(\hat{w}_i-w_i)^*(\hat{w}_j-w_j)x_i^*x_j\} \\&=(\sum_i|\delta_i|^2)\mathbb E\{|x_i|^2\}\\&= \Lambda(n)\sigma_x^2\\&=\Lambda(n)\frac{1}{N}\sum_i |x(n-i)|^2 \quad as \ N\to \infty\end{aligned}$

其中：由语音是白噪声这一假设可知 “交叉项” 为 0。

进而有： $\mu_{opt}(n) = \frac{1}{1+\frac{\sigma_v^2}{\sigma_r^2}}=\frac{\sigma_r^2}{\sigma_r^2+\sigma_v^2} = \frac{\sigma_r^2}{\sigma_e^2}$

这里又用到了残差信号与远端残留信号的不相关性，得到总误差的能量 $\sigma_e^2$ 就等于组成误差的两个分量的能量和 $\sigma_v^2+\sigma_r^2$ 。

结论：在以上的假设下，理论上的最优学习率 近似等于 残留回声 与 总误差信号 的能量比。

在实际操作中，迭代第 n 步有：（其中 $\hat{\sigma}_r(n),\hat{\sigma}_e(n)$ 分别是 $\sigma_r,\sigma_e$ 在第 n 步得估计）

$\hat{\mu}_{opt}(n) = \min(\frac{\hat{\sigma}_r^2(n)}{\hat{\sigma}_e^2(n)}, 1) \quad (10)$

4 量的估计

4.1 残留回声估计的上限

当条件 $\mathbb E\{\Lambda(n+1)\} = \Lambda(n)$ 满足时，我就达到了统计意义上的收敛（此时，滤波器在整体上达到了平稳状态，但系数之间或有细微变化），在公式 $(*)$ 中，令 $\mathbb E\{\Lambda(n+1)\} = \Lambda(n)$ ，并用远端信号得方差 $\sigma_x^2$ 替代其估计 $\frac{1}{N}\sum_i|x(n-i)|^2$ 后，可解得：

$\Lambda(n) \approx \frac{\sigma_v^2}{\sigma_x^2(\frac{2}{\mu}-1)} \quad (11)$

因此，如果当前滤波器失配度满足上式（11），则在统计意义下，使用最优学习率，将不再能改善滤波器失配度。最优学习率为 $\hat{\mu}_{opt} = \frac{\hat{\sigma}_r^2}{\hat{\sigma}_e^2}$ 代入（11）式中，考虑到误差信号得方差 $\sigma_e^2$ 可以估计得很准确，因此我们直接使用它得期望值代入，而残留回声方差 $\sigma_r^2$ 则用估计值 $\hat{\sigma}_r^2$ 代入，即将 $\frac{\hat{\sigma}_r^2}{\sigma_e^2}$ 作为学习率代入（11）式得：

$\begin{aligned}(11)式 &\Leftrightarrow \Lambda(n)\sigma_x^2 \approx\frac{\sigma_v^2}{\frac{2}{\mu} - 1} \\&\Leftrightarrow \sigma_r^2 \approx \frac{\sigma_v^2}{\frac{2\sigma_e^2}{\hat{\sigma}_r^2}-1} \\&\Leftrightarrow \sigma_r^2 \approx \frac{\sigma_v^2}{\frac{2(\sigma_r^2+\sigma_v^2)}{\hat{\sigma}_r^2}-1}\end{aligned}$

上式继续整理，得到如下方程：

$\begin{aligned}(11)式&\Leftrightarrow2(\sigma_r^2)^2+2\sigma_r^2\sigma_v^2-(\sigma_r^2+\sigma_v^2)\hat{\sigma}_r^2 = 0 \\&\Leftrightarrow (2\sigma_r^2-\hat{\sigma}_r^2)(\sigma_r^2-\sigma_v^2) = 0\end{aligned}$

所以，残差得估计值和期望值应满足： $\sigma_r^2 \approx \min(\frac{1}{2}\hat{\sigma}_r^2, \sigma_v^2) \quad (12)$

因此，在达到平稳后，残留回声的估计方差，不应该大于 2 倍的真实方差，即不能高估这个方差超过 2 倍（3 dB）。

显然，前面对近端语音 $v(n)$ 和远端语音 $x(n)$ 的白噪声假设是不合理的；但是，我们可以到频域进行估计，理由是：

a. 经观察，同一个频点沿着时间轴上的相关性比时域数据的相关性要低，这使得我们的假设在频域更加可靠。

b. 在频域，我们可以对每个频点单独计算最优学习率： $\hat{\mu}_{opt}(k, l) \approx \frac{\sigma_r^2(k,l)}{\sigma_e^2(k,l)} \quad (14)$

4.2 残留回声的估计

a. 假设滤波器有一个和频率无关的泄露系数 $\eta(l)$ ，则可以通过下式来对残留回声的方差进行估计：

$\hat{\sigma}_r^2(k,l) = \hat{\eta}(l)\hat{\sigma}_{\hat{Y}}^2(k,l) \quad (15)$

下面我们来证明（15）式的合理性，显然根据 $x(n)$ 的白噪假设，有：

$\begin{aligned}\mathbb E \{ \hat{y}^2 \} &= \mathbb E \{(\sum_i \hat{w}_i x(n-i))^2\}\\&= \sum_i \hat{w}_i^2 \mathbb E \{ x(n-i)^2\}\\&= (\sum_i \hat{w}_i^2) \sigma_x^2\end{aligned}$

由前面关于残留回声功率的结论，并代入上式有：

$\begin{aligned}\sigma_r^2 &= \Lambda \sigma_x^2 \\&= \frac{\Lambda}{\sum_i \hat{w}_i^2} \mathbb E\{\hat{y}^2\} \\&= \frac{\sum_i \delta_i^2}{\sum_i \hat{w}_i^2} \sigma_{\hat y}^2\end{aligned}$

在上式中，我们记 $\hat{\eta} = \frac{\sum_i \delta_i^2}{\sum_i \hat{w}_i^2}$ ；我们注意到 $\hat{\eta}$ 的分子和分母都是和滤波器权重的整体相关的一个值，是关于 echo path 以及 echo path change 的一个量；而在实际中回声路径被认为是缓慢变化的，所以我们可以认为 $\hat{\eta}$ 是一个缓慢变化的量。而另一项 $\sigma_{\hat y}^2$ 是语音的能量，因而是变化迅速的。

该过程同样也能说明频域的学习率也能分成这两个部分。

虽然在上式中对残留回声方差的估计，包含了不可计算的滤波器失配度 $\Lambda$ ，但是从式子来看， $\hat{\eta}$ 的确可是视为滤波器的 ERLE。并且作者给出了下图，说明 $\hat{\eta}$ 的倒数与 ERLE近似程度。

ERLE与

\eta

的一致性

使用上式来估计残留回声有两个好处：

1）将残留回声的估计转为对两个量的估计，一个是变化缓慢但不易估计的泄露系数 $\hat{\eta}$ ，另一个是变化快速但容易估计的回声功率 $\sigma_{\hat y}^2$ 。

2）在 double-talk 到来时，可以迅速做出反应：如下式的最优学习率，我们使用回声的瞬时功率 $|\hat{Y}(k,l)|^2$ 作为它的方差，同样使用误差信号的瞬时功率 $|E(k,l)|^2$ 作为误差信号的方差：

$\hat{\mu}_{opt}(k,l) = \min(\hat{\eta}(l)\frac{|\hat{Y}(k,l)|^2}{|E(k,l)|^2}, \mu_{max}) \quad (16)$

如此，就将学习率的自适应解耦成泄露系数 $\hat{\eta}$ 和回声功率与误差功率之比这两部分，正是后者提供了对 double-talk 的快速反应能力。（而泄露系数需要较长时间才能迭代到目标值）

另外，系统初始化时，滤波器系数为 0，所以 $\hat{Y}$ 也是 0，因此最优学习率总是 0；为了解决这个问题，考虑在初始化后使用固定的学习率（代码中是0.25）进行自适应滤波，等到滤波器系数收敛到一定程度时，再使用这里提出的自适应学习率算法。

4.3 泄露系数的估计

综上，我们就将残留回声方差估计的问题，转为对泄露系数 $\eta$ 的估计问题：

1）利用信号的非平稳性

2）利用估计的回声信号和误差信号之间的线性回归值来估计

3）基于的事实是：残留的回声与估计的回声是高度相关的，且估计的回声与误差信号中的噪声是无关的。

SpeexDSP 代码中的计算过程：

1）先使用一阶去 DC 滤波器去除估计的回声功率和误差功率谱的 DC分量，得到零均值的 $P_y,P_e$ ，代码中使用递归平均来估计均值，然后分别减去均值）

2）再使用下式计算平滑的相关：

$\begin{aligned}R_{EY} (k,l) &= (1-\beta(l))R_{EY}(k,l) + \beta(l) P_Y(k)P_E(k)\\R_{YY}(k,l) &= (1-\beta(l))R_{YY}(k,l) + \beta(l)P_Y(k)^2\end{aligned}$

其中， $\beta(l)$ 是一个动态的平滑系数，用于估计泄露系数，如下计算：

$\beta(l)= \beta_0\min(\frac{\hat{\sigma}_{\hat Y}^2}{\hat{\sigma}_e^2}, 1)$ ，其中 $\beta_0$ 是一个基础平滑系数

在 Speex 代码中，和前面的学习率一样，使用瞬时功率代替 $\beta(l)$ 估计式中的方差。

当远端没有信号时，有 beta = 0，所以 $R_{EY}$ 和 $R_{YY}$ 保持不变，即泄露系数保持不变，这样到下一次远端有语音（回声）过来时，能有一个合理的初始值。

3）最后将所有的频点上的相关性加起来，得到频点无关的泄露系数的估计：

$\hat{\eta}(l) = \frac{\sum_k R_{EY}(k,l)}{\sum_k R_{YY}(k,l)}$

5 双讲分析

分析前面的最优学习率公式，可以发现该系统能够处理 double-talk 场景：

a. 当近端语音突然出现时，公式（16）中的分母（瞬时功率）会陡然增大，从而使得学习率变得很小；这会持续到 double-talk 结束；

b. 假设近端没有语音，但是有平稳噪声，则（16）中的分母保持不变，此时学习率收到回声功率和泄露系数的影响：

i）当回声功率增大时，学习率变大，使得收敛迅速；

ii）当滤波器系数收敛时，泄露系数 $\eta$ 变小，从而学习率变小，使得收敛过程更平稳；

c. double-talk 时，残留回声与误差的相关性很小，从而保持学习率很小，不会造成滤波器发散；

d. Echo-path 发生变化时，滤波器系数完全失效，残留回声与误差有很大的相关性，从而使得学习率迅速达到 1，从而加快收敛。

6 代码解析

6.1 循环卷积到线性卷积

在代码实现时，使用FFT计算卷积，会产生循环卷机效应。为了看懂代码，非电子、信息、通信专业的工程师还是要搞清楚这一点，才能理解代码，下面先来分析一下为什么出现循环卷积：

$\begin{aligned}\mathcal F[{x(n)}] &= X(k) \\\mathcal F[w(n)] &= W(k)\\\end{aligned}$

其中： $\mathcal F$ 为傅立叶变换， $x(n),n=0,...,N-1$ 为某一帧数据， $w(n),n=0,...,N-1$ 为该时刻的滤波器权重。

期望使用卷积定理的等价形式 $\mathcal F^{-1}[\mathcal F[x(n)]\mathcal F[w(n)]]$ 来计算时域卷积运算（可参考任一本dsp方面的书籍），对于离散傅立叶变换的特殊情况，我们将其展开：

$\begin{aligned}y(m) &= \frac{1}{N}\sum_{k=0}^{n-1} X(k)W(k)e^{j2\pi km/N} \\&= \frac{1}{N}\sum_{k=0}^{n-1}\left (\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}\right) \left (\sum_{l=0}^{N-1}w(l)e^{-j2\pi kl/N}\right)e^{j2\pi km/N} \\&=\sum_{n=0}^{N-1}\sum_{l=0}^{N-1} \left(x(n)w(l) \frac{1}{N}\sum_{k=0}^{N-1}e^{-j2\pi k(n+l-m)/N} \right)\end{aligned} \tag{6.1.1}$

其中等比数列计算为： $\frac{1}{N} \sum_{k=0}^{N-1}e^{-j2\pi k(n+l-m)/N} = \begin{cases} 1 ,\quad n+l-m=pN, p\in\mathcal Z \\ 0 ,\quad\text{otherwise}\end{cases}$

因此，只有当 $l = m-n+pN$ 时，才包含有效的计算，上 $(6.1.1)$ 式为：

$y(m)=\sum_{n=0}^{N-1}x(n)w(m-n+pN), \quad p\in \mathcal Z, \ \mathbf {s.t.}\ 0 \le m-n+pN \le N-1$ ，这是一个循环卷积！见下面简单示例：

例如：当 $N=8$ 时，我们计算 $y(3)$ 的值有 $y(3) = \sum_{n=0}^{7}x(n)w(3-n+p\cdot 8)$ ，如下图所示 $x(n)$ 与 $w(l)$ 的对齐情况：

$\begin{bmatrix}x_0 &x_1& x_2& x_3&x_4&x_5&x_6&x_7\\w_3&w_2&w_1&w_0&w_7&w_6&w_5&w_4\end{bmatrix}$

这是一个非因果运算，显然我们的滤波器不允许这种情况发生，为了避免循环卷积，可以拼接前后两帧数据，并将权重的后半段置为 0 即可，多出来的 0 将不会产生有效的计算，如下图所示，产生 $y(0),\dots,y(7)$ 的过程：

$\begin{bmatrix}*&*&\cdots&*&x_0 & x_1& x_2& x_3&x_4&x_5&x_6&x_7&x_8&x_9&x_{10}&x_{11}&x_{12}&x_{13}&x_{14}&x_{15} \\0&0&\cdots&0&w_7&w_6&x_5&w_4&w_3&w_2&w_1&w_0& \rightarrow y(0) \\&0&\cdots&0&0&w_7&w_6&x_5&w_4&w_3&w_2&w_1&w_0& \rightarrow y(1) \\&&&&&&&&&&&&&\cdots\\&&&&&&&&0&0&\cdots&0&w_7&w_6&x_5&w_4&w_3&w_2&w_1&w_0& \rightarrow y(7) \\\end{bmatrix}$

注意：利用 $\mathcal F^{-1}[\mathcal F[x(n)]\mathcal F[w(n)]]$ 计算卷积，实际上的计算过程永远都是循环卷积，只不过在这里，我们通过将权重的后半段强行置 0 后，只有后半段的 $x(n)$ 在进行有效的计算，而此时得到的结果恰好就是线性卷积的结果。

（TODO）

最后编辑于：2024.03.08 20:00:22

禁止转载，如需转载请通过简信或评论联系作者。

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,793评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 87,567评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 151,342评论 0赞 338
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,825评论 1赞 277
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,814评论 5赞 368
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,680评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,033评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,687评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 42,175评论 1赞 300
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,668评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,775评论 1赞 332
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,419评论 4赞 321
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,020评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,978评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,206评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,092评论 2赞 351
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,510评论 2赞 343