对音频文件应用频率采样过滤器

Question

我正在尝试使用来自给定数据文件的频率采样来应用双带阻滤波器。 我使用的方法如下

1. 给定频率采样文件的反 fft (ifft)

2. 第 1 步给出的实数值循环移位

3. 对第 2 步的结果执行 fft

4. 使用卷积对音频文件应用过滤器。 （频域）

问题是带阻频率（925Hz 和 2090Hz）仍然存在。我的代码有什么问题还是我遗漏了什么？

[wave,fs]=audioread('audio.WAV'); 
data=importdata ('freqSampling.txt')

y=(ifft(data,401))
x=real (y)

r=circshift (x,200)

f=fft (r,4096);

new_sound=conv (wave, f)
sound(new_sound,fs,16);

谁能帮我解决这个问题？

Answer 1

前两步

1. 给定频率采样文件的反 fft (ifft)
2. 第 1 步给出的实数值循环移位

如果freqSampling.txt 文件正确指定了所需的完整双面频谱，则应该正确地为您提供根据您的规范构建的滤波器的时域系数（请参阅下面的“验证规范”）。如果频率规范包含陡峭的瞬变，ifft 点的数量也可能需要调整/增加。但是，如您在第 4 步中指出的那样在频域中执行卷积并不对应于典型的滤波操作，而是等效于两个信号的时域相乘。

时域过滤

根据第 2 步的结果，您可以使用 conv 直接在时域中过滤您的 wave 数据：

new_sound = conv(r, wave);

或filter:

new_sound = filter(r, 1, wave);

conv 将为您提供完整的length(wave)+length(r)-1 卷积，而filter 是一个更面向信号处理的函数，返回第一个length(wave) 卷积样本（并且还可以处理conv 所做的递归滤波器不直接支持）。

频域过滤

或者在频域中执行过滤

1. 使用至少为length(r)+length(wave)-1 的大小对步骤 2 的结果执行 FFT
2. 使用与步骤 3 中相同的大小对 wave 数据执行 FFT
3. 使用乘法（频域）对音频文件应用过滤器。
4. 计算第 5 步结果的逆 FFT (ifft)
5. 取第 6 步结果的实部（技术上不需要，但由于较小的数值舍入误差而经常需要）

这可以通过以下方式实现：

N = length(wave)+length(r)-1;
wave_fd = fft(wave, N);              % step 3
filter_fd = fft(r, N);               % step 4

filtered_fd = wave_fd .* filter_fd;  % step 5
new_sound = real(ifft(filtered_fd)); % step 6 & 7

请注意，您还可以使用the overlap-add method 以更小的块执行此频域过滤操作。

验证规范

根据您的 cmets，从 freqSampling.txt 文件导入的 data 可以通过以下方式重建：

N = 401;
data = ones(N,1);
data(19:23) = [2 1 0 1 2]/3;
data(51:56) = [2 1 0 1 2]/3;
data(N-[2:(N+1)/2]+2) = data([2:(N+1)/2]);

为了验证这是否会过滤所需的频率，我们可以将此规范绘制为频率的函数。为此，我们需要使用的采样率 (fs)，根据您的图表，这似乎是 22050。然后你应该能够用以下方式绘制这些：

hold off; plot([0:N-1]*fs/N, data);
hold on;  plot([925 925;2090 2090]', [0 1.2;0 1.2]', 'k:');
axis([0 3000 0 1.2]);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
legend('Specs', 'Tones');

这应该给出一个看起来像这样的情节：

基于此，规范似乎没有提供任何音调频率的衰减。可以通过以下方式构建更好的拟合：

N = 401;
data = ones(N,1);
data(round(925*N/fs)+1+[-2:2]) = [2 1 0 1 2]/3;  % data([16:20])
data(round(2090*N/fs)+1+[-2:2]) = [2 1 0 1 2]/3; % data([37:41])
data(N-[2:(N+1)/2]+2) = data([2:(N+1)/2]);

生成如下所示的验证图：

P.S.：根据您信号的频域图，第二个音似乎更接近 2600Hz，而不是指示的 2090Hz。