频率规整/频率弯折（Frequency Warping）原理 WFIR/WIIR滤波器设计matlab代码

最新推荐文章于 2024-05-04 21:55:24 发布

learnrocks100

最新推荐文章于 2024-05-04 21:55:24 发布

阅读量647

点赞数 5

分类专栏： DSP HRTF滤波器设计文章标签：频率弯折 warping frequency Warped FIR Warped IIR HRTF 滤波器设计

本文链接：https://blog.csdn.net/cyz_2014/article/details/136744524

版权

DSP 同时被 2 个专栏收录

13 篇文章 4 订阅

订阅专栏

HRTF滤波器设计

2 篇文章 0 订阅

订阅专栏

本文介绍了频率弯折（FrequencyWarping）在音频领域的应用，如HRTF滤波器设计，重点探讨了如何从FIR/IIR模型转换到WFIR和WIIR滤波器，包括权重函数和非线性频率刻度的使用，以及MATLAB中的实现代码示例。

摘要由CSDN通过智能技术生成

频率弯折（Frequency Warping）广泛应用在HRTF滤波器设计、扬声器均衡、以及Warped线性预测编码等音频领域中。

在HRTF滤波器设计过程中，我们最直接的方式是使用FIR模型。也可以将HRTF 的FIR模型转换为IIR模型，这里面有Prony方法（对应matlab的prony函数）、Yule-walker方法（对应matlab的yulewalker函数），还有一些方法如BMT（Balanced Model Truncation）。这些方法都是针对线性频率的。

人耳处理听觉信息是根据非均匀临界带频率分辨率，这意味着对HRTF建模也要遵循非线性的频率刻度。有两种方法可用于接近非线性频率刻度：

使用权重函数，使得高频的误差更大，低频拟合得更好；
使用非线性频率分辨率设计滤波器，即频率弯折（frequency warping）

1 频率弯折的原理

1.1 计算WFIR的系数

$H(z) = \sum_{n=0}^{\infty }h(n)z^{-n}$ ，

$h(n)$ 是脉冲响应， $H(z)$ 是传递函数

用 $z^{-1} = D(\bar{\textup{z}}) = \frac{\bar{\textup{z}}^{-1} + \lambda}{1 + \lambda \bar{\textup{z}}^{-1}}$ 替换 $H(z)$ 得到：

$H(\bar{z}) = \sum_{n=0}^{\infty }h(n)(\frac{\bar{\textup{z}}^{-1} + \lambda}{1 + \lambda \bar{\textup{z}}^{-1}})^{n}$

将 $H(\bar{z})$ 按照 $\bar{z}$ 展开得：

$H_{W}(\bar{z}) = \sum_{n=0}^{\infty }h_{W}(n)\bar{z}^{-n} = H(\bar{z})$

$h_{W}(n)$ 即为WFIR的系数

1.2 构造WFIR滤波器结构进行滤波

上图中的 $D_{1}(z) = \bar{z} = \frac{\textup{z}^{-1} - \lambda}{1 - \lambda \textup{z}^{-1}}$ ， $\beta$ 即是WFIR的系数 $h_{W}(n)$

1.3 获取 $h_{W}(n)$ 的方法

计算 $h_{W}(n)$ 有多种方法，这里介绍一种简单的方法：

将上图中的 $D_{1}(z)$ 替换为 $D(\bar{\textup{z}}) = \frac{\bar{\textup{z}}^{-1} + \lambda}{1 + \lambda \bar{\textup{z}}^{-1}}$ ， $\beta$ 替换为 $h(n)$ ，in是一个脉冲信号，out是所得到的WFIR系数。

2 频率弯折frequency warping的matlab实现代码

% 读取HRIR
[h1, fs1] = audioread('full\elev0\L0e270a.wav');
[~, h] = rceps(h1); % 获取最小相位序列

lambda = 0.7364;

% --------- 计算hw(n)-------------------
b = [lambda, 1];
a = [1, lambda];

imp = zeros(300,1);
imp(1) = 1;

Nw = 512;

ytemp = zeros(length(imp), Nw);
yout = zeros(size(imp));
ytemp(:, 1) = imp;

for i=2:Nw
    ytemp(:, i) = filter(b, a, ytemp(:, i-1));
end

for i=1:Nw
    yout = yout + ytemp(:, i)*h(i);
end
hw2 = yout;
%------------------计算hw(n)-------------

%------------------转换为WIIR------------
[BB, AA] = prony(hw2, 30, 30);
[zzb, ppb, kkb] = tf2zpk(BB, AA);
% unwarping
zz = (zzb + lambda) ./ (1 + lambda*zzb);
pp = (ppb + lambda) ./ (1 + lambda*ppb);

scale = db2mag(-2.7);
SOS = zp2sos(zz, pp, kkb*scale);
%------------------转换为WIIR------------


[BBn, AAn] = prony(h, 30, 30);
% fvtool(SOS, h,1);
% fvtool(BBn, AAn, h,1);
%------------绘图-------------------
[HH1, ww1] = freqz(h, 1, 512);
[HHprony, ww2] = freqz(BBn, AAn, 512);
[HHwarp, ww] = freqz(SOS, 512);
fhz1 = ww1/pi*fs1/2;
fhz2 = ww2/pi*fs1/2;
figure;
semilogx(fhz1,20*log10(abs(HH1)), fhz2,20*log10(abs(HHprony)),fhz2,20*log10(abs(HHwarp)), 'linewidth', 1)
xlabel('Frequency/Hz');
ylabel('Magnitude/dB');
legend('origin FIR', 'Prony', 'WIIR and Prony');
grid on;

实验结果：