GNU倍频程:实际FFT数据的1/N倍频程平滑(而不是它的表示)

Question

我想平滑一个脉冲响应音频文件。该文件的FFT显示，它是非常尖尖的。我想要平滑的音频文件，而不仅仅是它的情节，使我有一个更平滑的IR文件。我已经找到了一个函数，它显示了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工具箱提供的功能。

这是测试IR音频文件。

​

​

Answer 1

我想我找到了。由于音频文件(实数)的FFT是对称的，两边都有相同的实部，但虚部相反，我想这样做：

• 采取快速傅立叶变换，保留其中的一半，并应用平滑函数，而不将幅度转换为dB。
• 然后复制一份平滑的快速傅立叶变换，然后倒置想象的部分。
• 将这两个部分结合起来，使我具有与我在开始时相同的对称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);

)如果更有知识的人能够确认思想和代码是好的，那就太好了:)