# SFR

**Repository Path**: dirry/SFR

## Basic Information

- **Project Name**: SFR
- **Description**: 基于频域建模与自适应Tikhonov正则化的多通道声场重构
- **Primary Language**: Unknown
- **License**: MulanPSL-2.0
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2025-05-25
- **Last Updated**: 2025-06-09

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# 基于频域建模与自适应Tikhonov正则化的多通道声场重构

### A Multi-Channel Sound Field Reconstruction Algorithm Based on Frequency-Domain Modeling and Adaptive Regularization

> 作者：WSY 
>
> 日期：2024年5月   

---

## 项目简介

本项目实现了一套基于 **频域建模（Frequency-Domain Modeling）** 的多通道声场重构系统。系统首先通过已知的激励信号与响应数据计算扬声器到麦克风之间的频域传递函数，并利用正则化最小二乘法求解其逆传递函数。之后，将目标声场信号与逆传递函数进行卷积，生成用于重构声场的扬声器输出信号。整个过程在频域完成建模，在时域完成信号合成。系统还支持对输出结果进行仿真重建，通过与目标声场在1/3倍频程声压级上的误差评估重构精度。算法引入 Tikhonov 正则化以提高数值稳定性，并结合幅度控制机制实现对参数 β 的自适应调节，从而在控制扬声器最大输出幅值的同时保持最小声场重构误差。该方法适用于复杂声场的精准再现，具有良好的稳定性与鲁棒性。

应用背景包括虚拟声场构建、沉浸式音频再现、空间音响仿真、主动声学控制等场景。

---

## 文件结构

- `SFR_Base.py`：算法核心类，包含声场建模、逆滤波器计算、声场仿真验证等方法。
- `SFR_main.py`：主程序入口，支持命令行调用，提供帮助文档、主流程、调试功能等。
- `temp/`：用于存储中间结果、日志等。
- `wavPinkNoise_24bit_48k_20s.wav`：随机白噪音文件(见链接)
- `Channel_64.wav`：随机白噪音响因文件(见链接)
- `Airport_Handset.wav`：需要重构的目标音频文件(见链接)
- `Speaker_out.wav`：重构出来的结果音频文件(见链接)
- `wavLink.txt`：以上涉及到的所有wav文件的链接
- `README.md`:项目说明文档

---

## 核心算法流程

![SFR流程图](https://github.com/user-attachments/assets/7089b814-23ea-4094-835c-db8523a8978b)



本系统采用时域卷积形式实现声场重构，主要包括以下步骤：

### 1️.获取Random_in(t)和Response(t)

`Random_in(t)`信号为**白噪声/随机输入信号**,`shape`为`[8,960000]`,即$[l,fs\cdot Dur]$,

其中:$l$表示**声道数目(喇叭数目)**,$fs$表示**采样频率(程序设定为`48khz`)**,$Dur$表示**播音/采样时间(程序设定为`20s`)**

`Random_in(t)`信号的后`fs`个数据全为0,也就是将白噪声音频的最后一秒**置零**,目的在于**防止扬声器或系统末端音频冲击**，可以作为一个**“自然分界”或“哨兵标志”**有利于后续程序自动判断噪声信号是否已经播放完毕.

`Response(t)`信号为播放白噪声信号后采集的**响应信号**,`shape`为`[8,8,960000]`,即$[l,M,fs\cdot Dur]$,`M`表示用于采集的麦克风数目

`Random_in(t)`和`Response(t)`信号可从`White_IO.mat`文件读取,也可自己定义`Rando_in(t)`,再利用程序进行采集`Response(t)`

****

### 2.计算传递函数 H(f)和h(t)

显然我们认为实际实验环境为一个**线性时不变系(LTI)**系统,故可利用**PSD自功率谱密度** 和 **CPSD互功率谱密度** 计算 频域传递函数,用于估计线性时不变系统的**频率响应（transfer function）** 

对于 L 个扬声器和 M 个麦克风，系统的频域传递函数表示为：

$H_{ml}(f) = S_{xy}^{ml}(f) / S_{xx}(f)$

其中：

- $S_{xy}^{ml}(f)$：第 l 个扬声器与第 m 个麦克风之间的互功率谱（使用 CSD 计算）
- $S_{xx}(f)$：激励信号的自功率谱（使用 Welch 方法计算）

得到`H(f)`后利用**逆傅里叶变换**求解时域传递函数`h(t)`,但是需要注意,为了使频谱满足**赫尔曼对称**,即$X[−f]=X^∗[f]$ ,进行逆傅里叶变换的时候需要对频域传递函数`H(t)`进行**水平翻转**并进行**共轭处理**,得到完整的**双边频谱** 也就是包括正频率和负频率部分,才能保证时域信号为正确的实信号,程序中此类**对称共轭再逆傅里叶变换**的方法为自定义方法,同一集成于函数`symmetry_IFFT(self, Var, n)`

****

### 3.求解逆传递函数 C(f)和c(t)

通过正则化最小二乘法求解频域逆滤波器,为避免所得逆传递函数矩阵**病态**,引入**Tikhonov**正则化.：

$C(f) = (Hᴴ(f)·H(f) + β·I)^{-1} · Hᴴ(f)$

其中：

- `Hᴴ(f)`：H 的共轭转置
- `I`：单位矩阵
- `β`：正则化因子，控制逆解稳定性

团队经实验发现,此处**Tikhonov正则化算子**最优为**单位矩阵**,也就是此处的**TikHonov正则化特化为L2正则化**.

`β`经测试最优为$8e^{-5}$,故设置**正则化因子**的**初始值**为该值,后续根据结果自动调整

之后使用对称 IFFT 变换得到时域逆卷积核：

$c_{ml}(t) = IFFT_{symmetric}[C_{ml}(f)]$

这里的$IFFT_{symmetric}$指的就是此前计算`h(t)`用到的**对称共轭再逆傅里叶变换**的方法，即函数`symmetry_IFFT(self, Var, n)`

---

### 4.扬声器信号合成

目标声场为多通道时间序列$p_1(t)$，通过逆滤波器与目标信号进行卷积，得到扬声器输入：

$x_l(t) = ∑_{m=1}^{M} p_{1_m}(t) * c_{ml}(t)$

使用有效长度的 `valid` 卷积，避免边缘效应。

---

### 5.声场仿真与误差验证

将计算出的扬声器信号 $x_l(t)$ 与传递函数 $h_ml(t)$ 卷积，预测麦克风端重构声场：

$p̂_m(t) = ∑_{l=1}^{L} x_l(t) * h_{ml}(t)$

误差以 SPL（声压级）为指标，采用 1/3 倍频程积分方式比较：

$SPL_{error_{dB}} = mean(| SPL_{target} - SPL_{predicted} |)$

倍频程 SPL 计算使用 CPB 模块：

$SPL_i = 10·log10( ∑_{f ∈ Bi} |X(f)|² · Δf / (4 · 10^{-10}) )$

---

### 6.正则化参数 β 的自动调节

程序会自适应调整 **β**，使扬声器最大输出幅度趋于 1，同时优化误差：

- 若输出幅度 > 1：加大 β,对信号进行抑制
- 若输出幅度 < 1：减小 β,对信号进行增益

当出现连续的$β_1$和$β_2$使得$β_1\le 1$和$β_2\le 1$并且$\sqrt{(β_{1}-1)^{2}+(β_{2}-1)^{2}}\approx0$

则确定$[\min(\beta_1,\beta_2),\max(\beta_1,\beta_2)]$为**正则化因子**$\beta$的**粗调范围**

后续**细调**操作主要使用**二分搜索**

****

## 示例调用

```bash
# 全流程计算：包含目标声场，输出 Speaker_out.wav
python SFR_main.py wavPinkNoise.wav Channel_64.wav Airport_Handset.wav sfr.pkl Speaker_out.wav 14400 1e-4 0.5

# 调试模式：加载已有模型，重新计算逆滤波
python SFR_main.py -D 1e-3 0.5 sfr.pkl Airport_Handset.wav Speaker_out.wav

# 单独计算 1/3 倍频程数据
python SFR_main.py -F_PCB random_in.wav 48000 random_in_CPB.data
```

---

## 输出结果说明

- `sfr.pkl`：保存了传递函数、逆滤波器、采样率、滤波器长度等中间变量
- `Speaker_out.wav`：可直接播放的扬声器信号
- `log_data.txt`：记录每次迭代的 β、最大输出幅度、误差

---

## 依赖库

以下是项目中使用的主要第三方库及其版本：

| 库名称       | 版本         | 用途说明                              |
| ------------ | ------------ | ------------------------------------- |
| numpy        | 1.26.4       | 数组运算、信号处理基础                |
| scipy        | 1.15.2       | 信号处理：Welch功率谱、卷积、伪逆等   |
| soundfile    | 0.12.1       | 音频读写（支持 WAV 文件的加载与保存） |
| librosa      | 0.10.2.post1 | 提取 RMS 能量用于信号起止检测         |
| scikit-learn | 1.5.0        | KMeans 聚类用于起始点自动检测         |

安装布署依赖库：

```bash
pip install numpy scipy librosa soundfile scikit-learn
```

---

## 注意事项

- 所有 wav 文件应为 48kHz 单精度浮点型
- 默认输出路径在用代码目录下的 `./Temp`
- 软件有效期截至 2025 年 10 月 26 日，届时请联系作者续期，由于代码开源，可在`SFR_BASE`文件中的`check_EXPIRE_DAYS`函数中修改`EXPIRE_DAYS`变量

---

## 性能评估

系统性能通过对预测声场与目标声场的差异进行评价，核心指标为 SPL（声压级）误差，流程清晰、评估直观，适用于多通道声场控制系统的泛化评测。

---

### 1️.SPL 平均误差（SPL Error）

使用 1/3 倍频程分析方法对目标声场与重构声场在各麦克风位置进行声压级对比，定义误差如下：

$SPL_{error_{dB}} = mean( | SPL_{target} - SPL_{predicted} | )$

其中，SPL 使用 CPB 方法估算，频带覆盖 16 Hz 至 20 kHz，符合国际标准与听觉感知模型。


![result](https://github.com/user-attachments/assets/ee2c1663-eec8-47bb-9b12-6b5172050187)

---

### 2️.仿真评估流程

重构声场精度通过以下步骤完成验证：

1. **预测声场计算**：将扬声器信号与系统传递函数卷积：

   $p̂_m(t) = ∑_{l=1}^{L} x_l(t) * h_ml(t)$

2. **倍频程 SPL 转换**：使用 `CPB()` 提取声压级能量谱

3. **误差分析**：调用 `Validate()` 函数输出所有通道 SPL 差值的平均值

---

### 3️.指标意义

- SPL 平均误差$\lt  1.5 dB$：声场重构具有高度保真性
- 通道间误差一致性好：说明逆滤波器稳定、鲁棒
- β 自适应调节：有效平衡输出幅度与仿真误差

---

### 4️.日志输出样例（`log_data.txt`）

```
iter,beta,max_amp,Error_in_dB
1,1.00e-04,1.10,2.31
2,2.00e-04,0.95,1.72
3,1.50e-04,1.01,1.36
...
```

每轮迭代记录当前正则化因子、输出最大幅度、误差值，便于调参追踪与结果分析。



## 联系方式

- 联系邮箱：`202207020624@sust.edu.cn`
- 院校背景：陕西科技大学 计算机科学与技术
- 研究方向：信号处理、声场重构、机器学习与信号恢复