外耳/中耳传输函数的模拟

一、实验要求

参照ANSI S3.4-2007标准中的数据，利用Matlab或者Python，分别实现外耳和中耳的滤波功能（滤波器阶数自定，一般建议不超过3阶）。要求给出两个子函数，并用实现对一段实际音频的滤波和效果图示。

二、实验原理

1、外耳计算模型

建模实现从外耳耳廓到中耳鼓膜的频率响应（Transfer function of ear canal）

精细版本：头相关传递函数HRTF

简化版本：用低阶FIR滤波器模拟相似的频响效果

数据源自ANSI S3.4-2007标准。

表2.1.1 自由场到鼓膜的传递函数（正面发生率）。第一列表示声音频率（Hz），第二列表示鼓膜处电平减去自由场（无听者时）测到的电平（dB）。

图2.1.2 自由场到鼓膜的频率传输函数。纵坐标为鼓膜水平减去自由场的水平（dB），横坐标以kHz为单位。

2、中耳计算模型 

从中耳鼓膜到内耳耳蜗的频率响应（Transfer function of middle ear）

源自ANSI S3.4-2007标准

中耳建模通常有几种方法： 1，基于功能的模拟电路方法 2，基于切片的有限元建模分析 3，基于生理物理测量的数字滤波器

表2.2.1 中耳传递函数。第一列表示声音频率（Hz），第二列表示耳蜗相对于鼓膜的有效电平，斜体值表示未经验证的范围。

 

 图2.2.2 中耳频率传输函数。纵坐标为中耳的有效增益（dB），横坐标为频率（kHz），虚线为未验证范围。

3、外/中耳建模的数字滤波方法

已知频响，设计滤波器

基于窗函数的FIR设计方法： 1，设计频响 2，IDFT到时域 3，截取时域信号 4，加窗函数，获得实际的冲激响应 5，计算FIR滤波器系数 6，对冲激响应DFT得到实际频响

Matlab提供两个函数fir1（标准频响）和 fir2（任意频响）。

Python的SciPy中，signal库里也提供相应的firwin1和firwin2函数。

三、实验程序

%外耳传输频响

clear;

clc;

input=audioread('a.wav');

n=60;  %阶数

f1=[20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 750 800 1000 1250 1500 1600 2000 2500 3000 3150 4000 5000 6000 6300 8000 9000 10000 11200 12500 14000 15000 16000 20000];

[f,ps]=mapminmax(f1,0,1); %归一化

m1=[0 0 0 0 0 0 0 0 0.1 0.3 0.5 0.9 1.4 1.6 1.7 2.5 2.7 2.6 2.6 3.2 5.2 6.6 12.0 16.8 15.3 15.2 14.2 10.7 7.1 6.4 1.8 -0.9 -1.6 1.9 4.9 2.0 -2.0 2.5 2.5];

m=10.^(m1./20);

b=fir2(n,f,m); %已知频响的fir滤波器

output=filter(b,1,input);

[h,w] = freqz(b,1,128);   %求解数字滤波器频率响应

[h1,w1] = freqz(input,1,128);

[h2,w2] = freqz(output,1,128);

subplot(211)

plot(f,m,w/pi,abs(h))

legend('Ideal','fir2 Designed')

title('Transfer function of ear canal')

subplot(212)

plot(w1/pi,abs(h1),w2/pi,abs(h2));

legend('Input','Output')

ylabel('Amplititude');xlabel('Frenquency');title('Filter effect of ear canal');



%中耳传输频响

clear;

clc;

input=audioread('a.wav');

n=60;  %阶数

f1=[20.0 25.0 31.5 40.0 50.0 63.0 80.0 100.0 125.0 160.0 200.0 250.0 315.0 400.0 500.0 630.0 750.0 800.0 1000.0 1250.0 1500.0 1600.0 2000.0 2500.0 3000.0 3150.0 4000.0 5000.0 6000.0 6300.0 8000.0 9000.0 10000.0 11200.0 12500.0 14000.0 15000.0 16000.0 18000.0 20000.0];

[f,ps]=mapminmax(f1,0,1); %归一化

m1=[-39.6 -32.0 -25.85 -21.4 -18.5 -15.9 -14.1 -12.4 -11.0 -9.6 -8.3 -7.4 -6.2 -4.8 -3.8 -3.3 -2.9 -2.6 -2.6 -4.5 -5.4 -6.1 -8.5 -10.4 -7.3 -7.0 -6.6 -7.0 -9.2 -10.2 -12.2 -10.8 -10.1 -12.7 -15.0 -18.2 -23.8 -32.3 -45.5 -50.0];

m=10.^(m1./20);

b=fir2(n,f,m); %已知频响的fir滤波器

output=filter(b,1,input);

[h,w] = freqz(b,1,128);   %求解数字滤波器频率响应

[h1,w1] = freqz(input,1,128);

[h2,w2] = freqz(output,1,128);

subplot(211)

plot(f,m,w/pi,abs(h))

legend('Ideal','fir2 Designed')

title('Transfer function of middle ear')

subplot(212)

plot(w1/pi,abs(h1),w2/pi,abs(h2));

legend('Input','Output')

ylabel('Amplititude');xlabel('Frenquency');title('Filter effect of middle ear');

四、实验效果

图4.1 外耳传输函数和音频滤波效果图

图4.1 中耳传输函数和音频滤波效果图

不足：报告中不需要列出代码。原始数据要插值补充两头的频率。频响的展示要用dB幅度。滤波效果缺乏分析，以及听感描述。 优点是考虑了分贝数据到线性数据的转换。