GNU Octave：实际 FFT 数据的 1/N Octave 平滑不是它的表示

如何解决GNU Octave：实际 FFT 数据的 1/N Octave 平滑不是它的表示

我想要平滑脉冲响应音频文件。该文件的 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);

似乎有效！ :) 如果有更多知识渊博的人可以确认思想和代码是好的，那就太好了:)