Why do I still have these resonances in the lower frequences in the spectrograma and how do get rid of them in my fdn reverb?

Hello, as in the title I have problems with resonances in lower frequences and it makes the sound terrible. What exactly is the problem in my code? My code + spectrogram is attached below. I would really appreciate some help.

What I already tried: I tried to change the matrix to a hadamard and also tried to change the delay lines: smaller, bigger whatever i try its the same. Only the resonances change in frequences. The low pass filter I implemented already cut 95% of the resonances.

Heres my code and the spectrogram:

Code for plotting:

[x, fs] = audioread('Hier ist es.wav');
x = x(1:fs*1);
Nx = length(x);
nsc = 100;
nov = floor(nsc/2);
nff = max(256,2^nextpow2(nsc));

spectrogram(x,hanning(nsc),nov,nff);

% maxerr = max(abs(abs(t(:))-abs(x(:))));


ylabel("Frequency (Hz)")
xlabel("Time (s)")
title("Spectrogram delta impulse")

As you can see, in the frequences of 6850,4900,3100,2200 and 2400 for example there are horizontal lines, which I guess are resonances?

classdef FDNreverb2 < audioPlugin
    properties
        preDelayT = 0;              % pre delay [ms]
        reverbTime = 3.0;           % reverb time [s]
        roomSize = 15;
        mix = 65;                  % mix of wet and dry signal [Percent], 100 % -> only wet
        order = enumFDNorder.order8;                      % order of FDN, should be power of 2
        in_coeff = 1/2;         % coeff for in and output lines --- just for now one value for all
        out_coeff = 0.6;
    end

    properties (Access = private)
        N = 8;                  % order of FDN 
        a = 1.1;                    % multiplying factor for delays
        cSound = 343.2;             % speed of sound
        
        primeDelays = [2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 ...
            67 71 73 79 83 89 97 101 103 107 109 113 127 131];

        % following values have to be set manually for Interface
        maxPreDelay = 0.2;          % max pre delay [s] (200ms)
        maxRoomSize = 20.0;         % max room size [m] 

        % following variables initialized in reset method
        
        fs;         % sample rate
        A;                  % feedback matrix
        b;                  % input gain coefficients
        c;                  % output gain coefficients
        f;                  % feedback lines after matrix
        v;                  % signal before delay - v_i(n) = b_i * x(n) + f_i(n)
        s;                  % output lines
        d;                  % signal to be sent to matrix A       
        buffDelayLines;     % buffer delay lines, initialized in reset method
        m;                  % delay in samples of each delay line
        g;                  % gain coefficients of delay lines
        maxDelay;           % maximum delay in delay lines
        buffInput;          % input buffer for pre delay
        pBuffDelayLines;    % pointer for delay lines buffer
        pBuffInput;         % pointer for input buffer 
        preDelayS;          % preDelay in samples
        alpha;
        v_prev2;
        v_prev1;
        v_filt_prev1;
        v_filt_prev2;
        
    end

    properties(Constant)
        PluginInterface = audioPluginInterface( ...
            audioPluginParameter('preDelayT', ...
            'DisplayName', 'Pre Delay [ms]', ...
            'Mapping', {'int', 0, 200}, ...
            'Style', 'rotary', ...
            'Layout', [1,1]), ...
            audioPluginParameter('roomSize', ...
            'DisplayName', 'Room Size [m]', ...
            'Mapping', {'lin', 1.0, 20.0}, ...
            'Style', 'rotary', ...
            'Layout', [1,2]), ...
            audioPluginParameter('reverbTime', ...
            'DisplayName', 'Reverb Time [s]', ...
            'Mapping', {'lin', 0.1, 5.0}, ...
            'Style', 'rotary', ...
            'Layout', [1,3]), ...
            audioPluginParameter('order', ...
            'DisplayName', 'Order of FDN', ...
            'Mapping', {'enum', '8', '16', '32'}, ...
            'Layout', [3,2]), ...   
            audioPluginParameter('in_coeff', ...
            'DisplayName', 'Input Gain', ...
            'Mapping', {'lin', 0.0, 1.0}, ...
            'Style', 'rotary', ...
            'Layout', [3,3]), ...
            audioPluginParameter('out_coeff', ...
            'DisplayName', 'Output Gain', ...
            'Mapping', {'lin', 0.0, 1.0}, ...
            'Style', 'rotary', ...
            'Layout', [3,4]), ...
            audioPluginParameter('mix', ...
            'DisplayName', 'Mix', ...
            'Mapping', {'int', 0, 100}, ...
            'Style', 'rotary', ...
            'Layout', [3,1]), ...
            audioPluginGridLayout( ...
            'RowHeight', [100 100 100 100 100 100], ...
            'ColumnWidth',[100, 100 100 100 100 100], ...
            'RowSpacing', 10, ...
            'ColumnSpacing', 10, ...
            'Padding', [10 10 10 10]) ...
            );
    end

    methods

        function plugin = FDNreverb
            init(plugin)
        end

        function out = process(plugin, in) 
            out = zeros(size(in));
            %Butterworth coefficients
            [bx, ax] = butter(2, 3000/(plugin.fs/2)); % 2. Ordnung, Grenzfrequenz 1000 Hz
            b0 = bx(1);
            b1 = bx(2);
            b2 = bx(3);
            a1 = ax(2);
            a2 = ax(3);
             rawVn(1) = b0*plugin.v(1);
             rawVn(2) = b0*plugin.v(2) + b1*plugin.v(1) - a1* rawVn(1);
            
            
            
            % Initialize the previous output for the filter
            prevY = zeros(plugin.N, 1);
            
            for i = 1:size(in,1)
                % Summieren der L/R - Kan�le
                inL = in(i,1);
                inR = in(i,2);
                inSum = (inL + inR)/2;
                plugin.buffInput(plugin.pBuffInput + 1) = inSum;
                
                % loop over delay lines
                for n=1:plugin.N
                    % d_n = gain * delayed v_n
                    for k=1:plugin.N
                        plugin.d(k) = plugin.g(k) * plugin.buffDelayLines(k, mod(plugin.pBuffDelayLines + plugin.m(k), plugin.maxDelay +1) + 1);
                    end
                    % f_n = A(n,:) * d'
                    plugin.f(n) = plugin.A(n,:) * plugin.d(:);
                    % v_n with pre delay
                    rawVn = plugin.b(n) * plugin.buffInput(mod(plugin.pBuffInput + plugin.preDelayS, (plugin.maxPreDelay * plugin.fs + 1)) + 1) + plugin.f(n);
    
                    % Filter anwenden
                    plugin.v(n) = b0 * rawVn + b1 * plugin.v_prev1(n) + b2 * plugin.v_prev2(n) ...
                                - a1 * plugin.v_filt_prev1(n) - a2 * plugin.v_filt_prev2(n);
                    
                    % Update der vorherigen Filterwerte
                    plugin.v_prev2(n) = plugin.v_prev1(n);
                    plugin.v_prev1(n) = rawVn;
                    plugin.v_filt_prev2(n) = plugin.v_filt_prev1(n);
                    plugin.v_filt_prev1(n) = plugin.v(n);
                    
                    plugin.buffDelayLines(n, plugin.pBuffDelayLines + 1) = plugin.v(n);

                    % output lines
                    plugin.s(n) = plugin.c(n) * plugin.d(n);
                
                    out(i,:) = out(i,:) + real(plugin.s(n));
                    
                end

                % Assign to output
                out(i,1) = plugin.mix/100 * out(i,1) + (1.0 - plugin.mix/100) * in(i,1);
                out(i,2) = plugin.mix/100 * out(i,2) + (1.0 - plugin.mix/100) * in(i,2);

                calculatePointer(plugin);
            end 
        end

        function calculatePointer(plugin)
            if (plugin.pBuffDelayLines==0)
                plugin.pBuffDelayLines = plugin.maxDelay;
            else
                plugin.pBuffDelayLines = plugin.pBuffDelayLines - 1;
            end

            if (plugin.pBuffInput==0)
                plugin.pBuffInput = plugin.maxPreDelay * plugin.fs;
            else
                plugin.pBuffInput = plugin.pBuffInput - 1;
            end
        end

        function set.in_coeff(plugin, val)
            plugin.in_coeff = val;
            plugin.b = ones(plugin.N,1) * val;
        end

        function set.out_coeff(plugin, val)
            plugin.out_coeff = val;
            plugin.c = ones(plugin.N,1) * val;
        end

        function set.order(plugin, val)
            plugin.order = val;
            switch (plugin.order)
                case enumFDNorder.order8
                    plugin.N = 8;
                case enumFDNorder.order16
                    plugin.N = 16;
                case enumFDNorder.order32
                    plugin.N = 32;    
            end
            reset(plugin)
        end

        function set.reverbTime(plugin, val)
            plugin.reverbTime = val;
            calculateGainCoeffs(plugin, plugin.reverbTime)
            calculateDelays(plugin)
            %sprintf('Set Reverb Time: %f', plugin.reverbTime)
            %disp(['Set Reverb Time: ', num2str(plugin.reverbTime), ' s.']);
        end

        function set.preDelayT(plugin, val)
            plugin.preDelayT = val;
            calculatePreDelayS(plugin, plugin.preDelayT)
            %disp(['Set Predelay: ', int2str(plugin.preDelayT), ' ms.']);
        end

        function set.mix(plugin, val)
            plugin.mix = val;
            %disp(['Set Mix value: ', int2str(plugin.mix), ' %.']);
        end

        function calculateMatrix(plugin)
            plugin.A = 1*generateFDNmatrix(plugin.N);
        end
 
        function calculateDelays(plugin)
            % % calculate sample delays dependent on room size (val) and
            % % sample rate (fs)
            % % m_1 = trace of sound / cSound * fs
            % 
            M = ceil(0.15 * plugin.reverbTime * plugin.fs);
            % %disp(['M = ', int2str(M)]);
             plugin.m = zeros(plugin.N,1);
             interval = M/(plugin.N*4);
             %test = 0;
             for i=1:plugin.N
                 tmp = interval/2 + (i-1) * interval;
                 e = floor(0.5 + log(tmp)/log(plugin.primeDelays(i)));
                 plugin.m(i) = plugin.primeDelays(i)^(e);
             end
           
        end

        function calculateGainCoeffs(plugin, val)
            for i=1:plugin.N
                plugin.g(i) = 10^((-3) * (plugin.m(i)/plugin.fs) / val);
            end
            % Frequenzabhängige Dämpfung für tiefe Frequenzen
            cutoff_freq = 3000;  % z.B. 3000 Hz als Grenzfrequenz
            for i = 1:plugin.N
                freq = plugin.fs / plugin.m(i);  % Frequenz der Verzögerungslinie
                if freq < cutoff_freq
                plugin.g(i) = plugin.g(i) * 0.5;  % Dämpfe Gain für tiefe Frequenzen
                end
            end
        end

        function calculatePreDelayS(plugin, val)
            plugin.preDelayS = val * plugin.fs / 1000;
        end

        function init(plugin)
            plugin.fs = getSampleRate(plugin);

            % initialize buffDelayLines & pointer
            %plugin.maxDelay = floor(plugin.maxRoomSize/plugin.cSound * plugin.fs * plugin.a^(plugin.N)); %probably higher than actual max delay as actual delays get rounded down
            %plugin.maxDelay = ceil(0.15 * 5 * plugin.fs);   % maximum possible delay for max reverb time = 5s
            plugin.maxDelay = 131^2;
            plugin.buffDelayLines = zeros(plugin.N, plugin.maxDelay + 1);
            plugin.pBuffDelayLines = plugin.maxDelay;
            
            %initialize filter elements
            plugin.v_prev1 = zeros(plugin.N,1);
            plugin.v_prev2 = zeros(plugin.N,1);
            plugin.v_filt_prev2 = zeros(plugin.N,1);
            plugin.v_filt_prev1 = zeros(plugin.N,1);
            
            % initialize 
            plugin.A = zeros(plugin.N);
            plugin.b = ones(plugin.N,1) * plugin.in_coeff; % input gain coefficients
            plugin.c = ones(plugin.N,1) * plugin.out_coeff; % output gain coefficients
            plugin.f = zeros(plugin.N,1);             % feedback lines after matrix
            plugin.v = zeros(plugin.N,1);             % signal before delay - v_i(n) = b_i * x(n) + f_i(n)
            plugin.s = zeros(plugin.N,1);             % output lines
            plugin.d = zeros(plugin.N,1);             % signal to be sent fsto matrix A
            plugin.m = zeros(plugin.N,1);             % delay in samples of each delay line
            plugin.g = zeros(plugin.N,1);             % gain coefficients of delay lines
            plugin.buffInput = zeros((plugin.maxPreDelay * plugin.fs) + 1, 1);   % input buffer for pre delay, 0.2 = max pre delay of 200ms
            plugin.preDelayS = 0;                     % pre delay in samples
            plugin.pBuffInput = plugin.maxPreDelay * plugin.fs;    % pointer for input buffer

            % calculate sample delays of delay lines
            calculateDelays(plugin)

            % calculate gain coeffs of delay lines
            calculateGainCoeffs(plugin, plugin.reverbTime)
            
            % calculate FDN matrix based on set order
            calculateMatrix(plugin) 
        end

        function reset(plugin)
            init(plugin)
        end

        function generateFDNmatrix(order)
            % generates FDN matrix with evenly distributed eigenvalues at
            % the unit circle
            eig_nr = order/2;
            FDNmatrix = zeros(order, order);
            M = orth((rand(order,order)));
            BC_n = zeros(order, order, eig_nr);
            DtD_n = zeros(order, order, eig_nr);
            k1 = 1;
            for k=2:2:eig_nr*2
                x = M(:,k-1);
                y = M(:,k);
                x = x / sqrt(2);
                y = y / sqrt(2);
                B = x * x.';
                C = y * y.';
                D = x * y.';
                BC_n(1:order,1:order,k1) = B + C;
                DtD_n(1:order,1:order,k1) = D.' - D;
                w = pi * k / (eig_nr * 2 + 2);
                temp = 2 * BC_n(:,:,k1) * cos(w) + 2 * DtD_n(:,:,k1) * sin(w);
                FDNmatrix = FDNmatrix + temp;
                k1 = k1 + 1;
            end
        end
    end
end

Answer 1

[Ayush]

Hi @Muhsin Zerey,

I understand that you are implementing a Feedback Delay Network (FDN) reverb, and you’re experiencing issues with resonances at specific lower frequencies. This is a common challenge with FDN reverbs due to the nature of the delay lines and feedback matrix.

Here are some suggestions to address the resonances:

1. Feedback matrix: you can try using normalized Hadamard matrix or another orthogonal matrix like “householder” or “givens” rotation matrix.
2. The choice of prime numbers for delay line lengths is good for decorrelation, but the specific choice can still lead to resonances. You can try experimenting with different sets of prime numbers or slightly perturb the delays to break up resonant patterns.
3. You can try adding a spectral shaping stage after the FDN processing to further smoothen out resonances. This could be a dynamic EQ.
4. You can use a multi-band filter or a more complex damping curve to dampen high frequencies.

Here’s one addon, you can try implementing a simple orthogonal matrix for the feedback matrix using QR decomposition.

Here’s the code for your reference:

function generateFDNmatrix(plugin) 
    % Create an orthogonal matrix using QR decomposition 
    [Q, ~] = qr(randn(plugin.N)); 
    plugin.A = Q; 
end 

You can also try adjusting the window size (nsc), overlap (nov) and FFT points (nff) to get a clearer picture of the frequency content.

To read more about ”qr” decomposition, you can refer here: https://www.mathworks.com/help/matlab/ref/qr.html

Hope it helps!