问题描述
我想要平滑脉冲响应音频文件。该文件的 FFT 显示它非常刺耳。我想平滑音频文件,而不仅仅是它的情节,以便我有一个更平滑的 IR 文件。 我发现 a function 显示平滑的 FFT 图。这种平滑如何应用于实际的 FFT 数据,而不仅仅是它的绘图?
[y,Fs] = audioread('test\test IR.wav');
function x_oct = smoothSpectrum(X,f,Noct)
%sMOOTHSPECTRUM Apply 1/N-octave smoothing to a frequency spectrum
%% Input checking
assert(isvector(X),'smoothSpectrum:invalidX','X must be a vector.');
assert(isvector(f),'smoothSpectrum:invalidF','F must be a vector.');
assert(isscalar(Noct),'smoothSpectrum:invalidNoct','NOCT must be a scalar.');
assert(isreal(X),'X must be real.');
assert(all(f>=0),'F must contain positive values.');
assert(Noct>=0,'NOCT must be greater than or equal to 0.');
assert(isequal(size(X),size(f)),'smoothSpectrum:invalidInput','X and F must be the same size.');
%% Smoothing
% calculates a Gaussian function for each frequency,deriving a
% bandwidth for that frequency
x_oct = X; % initial spectrum
if Noct > 0 % don't bother if no smoothing
for i = find(f>0,1,'first'):length(f)
g = gauss_f(f,f(i),Noct);
x_oct(i) = sum(g.*X); % calculate smoothed spectral coefficient
end
% remove undershoot when X is positive
if all(X>=0)
x_oct(x_oct<0) = 0;
end
end
endfunction
function g = gauss_f(f_x,F,Noct)
% GAUSS_F calculate frequency-domain Gaussian with unity gain
%
% G = GAUSS_F(F_X,NOCT) calculates a frequency-domain Gaussian function
% for frequencies F_X,with centre frequency F and bandwidth F/NOCT.
sigma = (F/Noct)/pi; % standard deviation
g = exp(-(((f_x-F).^2)./(2.*(sigma^2)))); % Gaussian
g = g./sum(g); % normalise magnitude
endfunction
% take fft
Y = fft(y);
% keep only meaningful frequencies
nfft = length(y);
if mod(nfft,2)==0
Nout = (nfft/2)+1;
else
Nout = (nfft+1)/2;
end
Y = Y(1:Nout);
f = ((0:Nout-1)'./nfft).*Fs;
% put into dB
Y = 20*log10(abs(Y)./nfft);
% smooth
Noct = 12;
Z = smoothSpectrum(Y,Noct);
% plot
semilogx(f,Y,'linewidth',0.7,Z,2.2);
xlim([20,20000])
grid on
附注。我有 Octave GNU,所以我没有 Matlab 工具箱提供的功能。
Here is the test IR audio file.
解决方法
我想我找到了。由于音频文件的FFT(实数)是对称的,两边实部相同但虚部相反,我想到了这样做:
- 取 FFT,保留其一半,然后应用平滑函数而不将幅度转换为 dB
- 然后复制平滑的 FFT,并只反转虚部
- 将这两个部分结合起来,这样我就有了与开始时相同的对称 FFT,但现在它变得平滑了
- 对此应用逆 FFT 并取实部并将其写入文件。
代码如下:
[y,Fs] = audioread('test IR.wav');
function x_oct = smoothSpectrum(X,f,Noct)
x_oct = X; % initial spectrum
if Noct > 0 % don't bother if no smoothing
for i = find(f>0,1,'first'):length(f)
g = gauss_f(f,f(i),Noct);
x_oct(i) = sum(g.*X); % calculate smoothed spectral coefficient
end
% remove undershoot when X is positive
if all(X>=0)
x_oct(x_oct<0) = 0;
end
end
endfunction
function g = gauss_f(f_x,F,Noct)
sigma = (F/Noct)/pi; % standard deviation
g = exp(-(((f_x-F).^2)./(2.*(sigma^2)))); % Gaussian
g = g./sum(g); % normalise magnitude
endfunction
% take fft
Y = fft(y);
% keep only meaningful frequencies
NFFT = length(y);
if mod(NFFT,2)==0
Nout = (NFFT/2)+1;
else
Nout = (NFFT+1)/2;
end
Y = Y(1:Nout);
f = ((0:Nout-1)'./NFFT).*Fs;
% smooth
Noct = 12;
Z = smoothSpectrum(Y,Noct);
% plot
semilogx(f,Y,'LineWidth',0.7,Z,2.2);
xlim([20,20000])
grid on
#Apply the smoothing to the actual data
Zreal = real(Z); # real part
Zimag_neg = Zreal - Z; # opposite of imaginary part
Zneg = Zreal + Zimag_neg; # will be used for the symmetric Z
# Z + its symmetry with same real part but opposite imaginary part
reconstructed = [Z ; Zneg(end-1:-1:2)];
# Take the real part of the inverse FFT
reconstructed = real(ifft(reconstructed));
#Write to file
audiowrite ('smoothIR.wav',reconstructed,Fs,'BitsPerSample',24);
似乎有效! :) 如果有更多知识渊博的人可以确认思想和代码是好的,那就太好了:)