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

Question

我想平滑脉冲响应音频文件。该文件的 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), 'smoothSpectrum:invalidX', 'X must be real.');
    assert(all(f>=0), 'smoothSpectrum:invalidF', 'F must contain positive values.');
    assert(Noct>=0, 'smoothSpectrum:invalidNoct', '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,F,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,f,Noct);
% plot
semilogx(f,Y,'LineWidth',0.7,f,Z,'LineWidth',2.2);
xlim([20,20000])
grid on

PS。我有 Octave GNU，所以我没有 Matlab 工具箱可用的功能。

Here is the test IR audio file.

Answer 1

我想我找到了。由于音频文件（实数）的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,f,Noct);

% plot
semilogx(f,Y,'LineWidth',0.7,f,Z,'LineWidth',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);

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

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

GNU Octave: 1/N Octave Smoothing of actual FFT Data (not the representation of it)

octave

fft

gnu

smoothing