Audio matlab challenge: Reverberation, FDN, Filtering

Question

To give a bit of context, i am trying to implement in matlab from scratch the following signal diagram, a Feedback Delay Network (FDN). pic: FDN

With an appropriate matrix, indifferent to delay lengths, virtually white noise comes out when fed a dirac impulse.

I've managed to do this in code, but my goal is another and hence my question. I want to apply a filter h(z) after each delay line z^-m. pic: h(z)

More specifically, i want to apply a third-octave cascaded graphic equalizer after each delay line. The purpose is to create frequency dependent attenuation on the whole structure, and consequently delay dependent. I've successfully designed the filter in the form of SOS, but my problem is: how do I apply it within the structure? I assume to use sosfilt() somewhere with what I have, but I'm not sure.

I haven't reduced the order of the system for sake of purpose. The order is 16 (16x16 matrix, 16 delay lines, 31x16 biquad filters)

The first code refers to the lossless FDN, safely runnable which generates white noise. I have commented my failed attempt to introduce the filtering in the loop saying: % Filtering

Unfortunately, I can't post all GEQ entries, but I'll leave 8 in the end corresponding to the first 8 delays.

So, the question is how do I code to filter the white noise, implementing frequency dependent attenuation in the whole FDN structure. Also, although it may be computationally inefficient, I'd prefer to apply this without higher level functions and based on what I already have, i.e: applicable in GNU Octave

Edit: Assuming you have to apply the bandpass 2nd order filtering sample by sample using the difference equation, how would you recursively do it for 31 bands in series? One is shown in the second code section.

% Feedback Delay Network - Mono

fs = 44100;
in = [ 1; 0 ];            % Dirac Impulse
in = [in; zeros(3*fs,1)]; % Space for reverb
 
% Householder matrix N=16
A = eye(4) - ( 2/4 * ones(4,1) * ones(1,4) );
MTX = .5 * [A -A -A -A; -A A -A -A; -A -A A -A; -A -A -A A];

N = size(MTX,1);    % Matrix order

delays = [1500 1547 1602 1668 1745 1838 1947 2078 2232,...
   2415 2623 2890 3196 3559 3989 4500]; % N=16 delays in samples

load('GEQ.mat');    % Third octave graphic equalizer calculated based
                    % on an atenuation-per-sample and scaled by delay. 
                    % To be applied before or after each delayline
                   
% Initialize signals
delayBuffer = max(delays) + max(delays)/10;
bufferN = zeros(delayBuffer,N);   % Delay buffers
FB = zeros(1,N);               % Feedback matrix output
in_dl = zeros(1,N);            % Delay input
out_dl = zeros(1,N);           % Delay output
nSample = length(in);          % Number of samples
out = zeros(nSample,1);        % Output
 

% FDN Computation
for sample = 1:nSample     % each sample
    for n = 1:N            % each delayline
         
       in_dl(n) = in(sample,1) + FB(n); % Input and Feedback matrix sum
        
       [out_dl(n),bufferN(:,n)] = funcDelay( in_dl(n), bufferN(:,n),...
                           sample, delays(n) ); % Delaying
       % Filtering
       %out_dl(n) = sosfilt( GEQ(:,:,n), out_dl(n) );  
    end
     
    out(sample,1) = 1/sqrt(2) * sum(out_dl); % Delay output sum
     
    FB = out_dl * MTX; % Feedback matrix output recalculation
end

% Used function
function [out,buffer] = funcDelay(in,buffer,n,delay)
  
  % Circular buffer indices
  len = length(buffer);
  indexC = mod(n-1, len) + 1;        % Current
  indexD = mod(n-delay-1, len) + 1;  % Delay
  
  out = buffer(indexD,1);
  
  % Stores output on appropriate index
  buffer(indexC,1) = in; 
end
 %sound(out,fs,16)

Second code section: applying filter difference equation to signal. Suggestions for coding it for 31 filters recursively?

in = (rand(1,100)*2)-1; % Example noise 100 samples
samples = length(in);
out = zeros(1,samples);

b0 = GEQ(1,1,1);      % Coeffs extracted from actual GEQ
b1 = GEQ(1,2,1);    a1 = GEQ(1,5,1);
b2 = GEQ(1,3,1);    a2 = GEQ(1,6,1);

out(1) = b0 * in(1);
out(2) = b0 * in(2) + b1 * in(1) - a1 * out(1);
for n = 3:samples
    out(n) = b0*in(n) + b1*in(n-1) + b2*in(n-2) - a1*out(n-1) - a2*out(n-2);
end

Thanks!

GEQ(:,:,1) = [0.6444   -1.2882    0.6438    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9984    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9980    1.0000   -1.9978    0.9979;
    1.0000   -1.9975    0.9975    1.0000   -1.9973    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9958    0.9960    1.0000   -1.9957    0.9959;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9935    0.9938    1.0000   -1.9929    0.9932;
    1.0000   -1.9915    0.9920    1.0000   -1.9912    0.9917;
    1.0000   -1.9892    0.9900    1.0000   -1.9887    0.9895;
    1.0000   -1.9861    0.9873    1.0000   -1.9856    0.9868;
    1.0000   -1.9825    0.9845    1.0000   -1.9810    0.9830;
    1.0000   -1.9771    0.9802    1.0000   -1.9757    0.9789;
    1.0000   -1.9702    0.9752    1.0000   -1.9685    0.9735;
    1.0000   -1.9609    0.9688    1.0000   -1.9587    0.9667;
    1.0000   -1.9483    0.9608    1.0000   -1.9458    0.9584;
    1.0000   -1.9311    0.9508    1.0000   -1.9281    0.9478;
    1.0000   -1.9070    0.9381    1.0000   -1.9039    0.9350;
    1.0000   -1.8738    0.9228    1.0000   -1.8698    0.9187;
    1.0000   -1.8275    0.9043    1.0000   -1.8215    0.8980;
    1.0000   -1.7608    0.8807    1.0000   -1.7538    0.8732;
    1.0000   -1.6659    0.8520    1.0000   -1.6580    0.8432;
    1.0000   -1.5308    0.8178    1.0000   -1.5209    0.8059;
    1.0000   -1.3382    0.7756    1.0000   -1.3278    0.7616;
    1.0000   -1.0671    0.7226    1.0000   -1.0607    0.7118;
    1.0000   -0.7061    0.6745    1.0000   -0.6929    0.6388;
    1.0000   -0.2324    0.6083    1.0000   -0.2311    0.5703;
    1.0000    0.3354    0.5587    1.0000    0.3047    0.4869;
    1.0000    0.9408    0.5246    1.0000    0.8392    0.4163;
    1.0000    1.5310    0.6212    1.0000    1.2251    0.3584];
GEQ(:,:,2) = [0.6356   -1.2706    0.6350    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9984    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9980    1.0000   -1.9978    0.9979;
    1.0000   -1.9975    0.9975    1.0000   -1.9973    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9959    0.9960    1.0000   -1.9957    0.9958;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9935    0.9938    1.0000   -1.9929    0.9932;
    1.0000   -1.9915    0.9920    1.0000   -1.9912    0.9917;
    1.0000   -1.9892    0.9900    1.0000   -1.9887    0.9895;
    1.0000   -1.9861    0.9873    1.0000   -1.9856    0.9868;
    1.0000   -1.9825    0.9845    1.0000   -1.9810    0.9830;
    1.0000   -1.9771    0.9803    1.0000   -1.9757    0.9789;
    1.0000   -1.9702    0.9752    1.0000   -1.9684    0.9734;
    1.0000   -1.9609    0.9689    1.0000   -1.9587    0.9666;
    1.0000   -1.9483    0.9608    1.0000   -1.9458    0.9583;
    1.0000   -1.9311    0.9509    1.0000   -1.9280    0.9478;
    1.0000   -1.9070    0.9382    1.0000   -1.9039    0.9350;
    1.0000   -1.8739    0.9228    1.0000   -1.8697    0.9186;
    1.0000   -1.8276    0.9044    1.0000   -1.8214    0.8979;
    1.0000   -1.7609    0.8808    1.0000   -1.7537    0.8731;
    1.0000   -1.6660    0.8522    1.0000   -1.6579    0.8431;
    1.0000   -1.5310    0.8180    1.0000   -1.5208    0.8057;
    1.0000   -1.3384    0.7758    1.0000   -1.3276    0.7614;
    1.0000   -1.0672    0.7227    1.0000   -1.0606    0.7116;
    1.0000   -0.7063    0.6751    1.0000   -0.6927    0.6382;
    1.0000   -0.2324    0.6088    1.0000   -0.2311    0.5697;
    1.0000    0.3359    0.5598    1.0000    0.3042    0.4858;
    1.0000    0.9423    0.5263    1.0000    0.8375    0.4146;
    1.0000    1.5349    0.6247    1.0000    1.2195    0.3537];
GEQ(:,:,3) = [0.6255   -1.2504    0.6249    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9984    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9980    1.0000   -1.9978    0.9979;
    1.0000   -1.9975    0.9975    1.0000   -1.9973    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9959    0.9960    1.0000   -1.9957    0.9958;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9935    0.9938    1.0000   -1.9929    0.9932;
    1.0000   -1.9915    0.9920    1.0000   -1.9912    0.9917;
    1.0000   -1.9892    0.9900    1.0000   -1.9887    0.9895;
    1.0000   -1.9861    0.9873    1.0000   -1.9856    0.9868;
    1.0000   -1.9825    0.9845    1.0000   -1.9810    0.9830;
    1.0000   -1.9771    0.9803    1.0000   -1.9757    0.9788;
    1.0000   -1.9702    0.9752    1.0000   -1.9684    0.9734;
    1.0000   -1.9610    0.9689    1.0000   -1.9587    0.9666;
    1.0000   -1.9483    0.9609    1.0000   -1.9458    0.9583;
    1.0000   -1.9312    0.9509    1.0000   -1.9280    0.9477;
    1.0000   -1.9071    0.9382    1.0000   -1.9038    0.9349;
    1.0000   -1.8739    0.9229    1.0000   -1.8697    0.9185;
    1.0000   -1.8277    0.9045    1.0000   -1.8213    0.8978;
    1.0000   -1.7610    0.8810    1.0000   -1.7536    0.8730;
    1.0000   -1.6662    0.8523    1.0000   -1.6577    0.8429;
    1.0000   -1.5312    0.8182    1.0000   -1.5206    0.8055;
    1.0000   -1.3385    0.7761    1.0000   -1.3274    0.7612;
    1.0000   -1.0674    0.7229    1.0000   -1.0605    0.7115;
    1.0000   -0.7066    0.6757    1.0000   -0.6925    0.6376;
    1.0000   -0.2324    0.6095    1.0000   -0.2311    0.5690;
    1.0000    0.3363    0.5610    1.0000    0.3037    0.4845;
    1.0000    0.9440    0.5282    1.0000    0.8355    0.4126;
    1.0000    1.5395    0.6288    1.0000    1.2128    0.3482];
GEQ(:,:,4) = [0.6136   -1.2265    0.6130    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9984    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9980    1.0000   -1.9978    0.9979;
    1.0000   -1.9975    0.9975    1.0000   -1.9973    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9959    0.9960    1.0000   -1.9957    0.9958;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9935    0.9939    1.0000   -1.9929    0.9932;
    1.0000   -1.9915    0.9920    1.0000   -1.9912    0.9917;
    1.0000   -1.9892    0.9900    1.0000   -1.9887    0.9895;
    1.0000   -1.9861    0.9874    1.0000   -1.9855    0.9868;
    1.0000   -1.9825    0.9845    1.0000   -1.9809    0.9829;
    1.0000   -1.9771    0.9803    1.0000   -1.9757    0.9788;
    1.0000   -1.9702    0.9753    1.0000   -1.9684    0.9734;
    1.0000   -1.9610    0.9689    1.0000   -1.9586    0.9665;
    1.0000   -1.9484    0.9609    1.0000   -1.9457    0.9582;
    1.0000   -1.9312    0.9510    1.0000   -1.9279    0.9476;
    1.0000   -1.9072    0.9383    1.0000   -1.9037    0.9348;
    1.0000   -1.8740    0.9230    1.0000   -1.8696    0.9184;
    1.0000   -1.8278    0.9046    1.0000   -1.8211    0.8977;
    1.0000   -1.7612    0.8811    1.0000   -1.7534    0.8728;
    1.0000   -1.6663    0.8525    1.0000   -1.6575    0.8427;
    1.0000   -1.5314    0.8184    1.0000   -1.5204    0.8053;
    1.0000   -1.3388    0.7764    1.0000   -1.3272    0.7609;
    1.0000   -1.0675    0.7232    1.0000   -1.0604    0.7112;
    1.0000   -0.7069    0.6764    1.0000   -0.6922    0.6368;
    1.0000   -0.2325    0.6103    1.0000   -0.2310    0.5681;
    1.0000    0.3369    0.5624    1.0000    0.3030    0.4830;
    1.0000    0.9460    0.5304    1.0000    0.8332    0.4102;
    1.0000    1.5449    0.6336    1.0000    1.2047    0.3416];
GEQ(:,:,5) = [0.5999   -1.1993    0.5993    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9985    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9980    1.0000   -1.9978    0.9979;
    1.0000   -1.9975    0.9975    1.0000   -1.9973    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9959    0.9960    1.0000   -1.9957    0.9958;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9935    0.9939    1.0000   -1.9928    0.9932;
    1.0000   -1.9915    0.9920    1.0000   -1.9912    0.9917;
    1.0000   -1.9892    0.9900    1.0000   -1.9887    0.9894;
    1.0000   -1.9861    0.9874    1.0000   -1.9855    0.9868;
    1.0000   -1.9826    0.9846    1.0000   -1.9809    0.9829;
    1.0000   -1.9772    0.9803    1.0000   -1.9756    0.9788;
    1.0000   -1.9703    0.9753    1.0000   -1.9683    0.9733;
    1.0000   -1.9611    0.9690    1.0000   -1.9586    0.9665;
    1.0000   -1.9485    0.9610    1.0000   -1.9456    0.9582;
    1.0000   -1.9313    0.9511    1.0000   -1.9278    0.9475;
    1.0000   -1.9072    0.9384    1.0000   -1.9037    0.9348;
    1.0000   -1.8741    0.9231    1.0000   -1.8695    0.9183;
    1.0000   -1.8280    0.9048    1.0000   -1.8210    0.8975;
    1.0000   -1.7614    0.8813    1.0000   -1.7532    0.8726;
    1.0000   -1.6665    0.8527    1.0000   -1.6573    0.8425;
    1.0000   -1.5316    0.8187    1.0000   -1.5201    0.8050;
    1.0000   -1.3390    0.7767    1.0000   -1.3270    0.7605;
    1.0000   -1.0677    0.7234    1.0000   -1.0602    0.7109;
    1.0000   -0.7072    0.6773    1.0000   -0.6918    0.6359;
    1.0000   -0.2325    0.6112    1.0000   -0.2310    0.5672;
    1.0000    0.3376    0.5640    1.0000    0.3022    0.4813;
    1.0000    0.9484    0.5331    1.0000    0.8304    0.4073;
    1.0000    1.5511    0.6393    1.0000    1.1953    0.3338];
GEQ(:,:,6) = [0.5839   -1.1672    0.5833    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9985    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9981    1.0000   -1.9978    0.9979;
    1.0000   -1.9975    0.9975    1.0000   -1.9972    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9959    0.9960    1.0000   -1.9957    0.9958;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9936    0.9939    1.0000   -1.9928    0.9931;
    1.0000   -1.9915    0.9920    1.0000   -1.9912    0.9917;
    1.0000   -1.9892    0.9900    1.0000   -1.9886    0.9894;
    1.0000   -1.9861    0.9874    1.0000   -1.9855    0.9868;
    1.0000   -1.9826    0.9846    1.0000   -1.9808    0.9828;
    1.0000   -1.9772    0.9804    1.0000   -1.9756    0.9787;
    1.0000   -1.9703    0.9753    1.0000   -1.9683    0.9733;
    1.0000   -1.9611    0.9691    1.0000   -1.9585    0.9664;
    1.0000   -1.9485    0.9611    1.0000   -1.9456    0.9581;
    1.0000   -1.9314    0.9512    1.0000   -1.9277    0.9475;
    1.0000   -1.9073    0.9385    1.0000   -1.9036    0.9347;
    1.0000   -1.8742    0.9232    1.0000   -1.8693    0.9182;
    1.0000   -1.8281    0.9050    1.0000   -1.8208    0.8973;
    1.0000   -1.7616    0.8815    1.0000   -1.7530    0.8724;
    1.0000   -1.6668    0.8530    1.0000   -1.6571    0.8422;
    1.0000   -1.5319    0.8191    1.0000   -1.5198    0.8046;
    1.0000   -1.3393    0.7771    1.0000   -1.3266    0.7601;
    1.0000   -1.0678    0.7237    1.0000   -1.0600    0.7106;
    1.0000   -0.7076    0.6783    1.0000   -0.6914    0.6348;
    1.0000   -0.2325    0.6123    1.0000   -0.2310    0.5660;
    1.0000    0.3384    0.5659    1.0000    0.3013    0.4792;
    1.0000    0.9513    0.5363    1.0000    0.8271    0.4039;
    1.0000    1.5586    0.6460    1.0000    1.1838    0.3244];
GEQ(:,:,7) = [0.5656   -1.1306    0.5651    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9985    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9981    1.0000   -1.9978    0.9978;
    1.0000   -1.9975    0.9976    1.0000   -1.9972    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9959    0.9960    1.0000   -1.9957    0.9958;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9936    0.9939    1.0000   -1.9928    0.9931;
    1.0000   -1.9915    0.9920    1.0000   -1.9912    0.9917;
    1.0000   -1.9893    0.9900    1.0000   -1.9886    0.9894;
    1.0000   -1.9861    0.9874    1.0000   -1.9855    0.9868;
    1.0000   -1.9827    0.9847    1.0000   -1.9808    0.9828;
    1.0000   -1.9773    0.9804    1.0000   -1.9755    0.9787;
    1.0000   -1.9704    0.9754    1.0000   -1.9682    0.9732;
    1.0000   -1.9612    0.9691    1.0000   -1.9584    0.9663;
    1.0000   -1.9486    0.9611    1.0000   -1.9455    0.9580;
    1.0000   -1.9315    0.9513    1.0000   -1.9276    0.9473;
    1.0000   -1.9074    0.9386    1.0000   -1.9035    0.9345;
    1.0000   -1.8744    0.9234    1.0000   -1.8692    0.9180;
    1.0000   -1.8283    0.9052    1.0000   -1.8206    0.8971;
    1.0000   -1.7618    0.8818    1.0000   -1.7527    0.8721;
    1.0000   -1.6670    0.8533    1.0000   -1.6568    0.8419;
    1.0000   -1.5323    0.8195    1.0000   -1.5194    0.8041;
    1.0000   -1.3397    0.7776    1.0000   -1.3262    0.7595;
    1.0000   -1.0681    0.7241    1.0000   -1.0598    0.7102;
    1.0000   -0.7080    0.6795    1.0000   -0.6910    0.6335;
    1.0000   -0.2326    0.6136    1.0000   -0.2309    0.5646;
    1.0000    0.3393    0.5682    1.0000    0.3002    0.4768;
    1.0000    0.9546    0.5400    1.0000    0.8231    0.3999;
    1.0000    1.5672    0.6537    1.0000    1.1702    0.3133];
GEQ(:,:,8) = [0.5444   -1.0883    0.5439    1.0000   -1.9989    0.9989;
    1.0000   -1.9987    0.9987    1.0000   -1.9987    0.9987;
    1.0000   -1.9984    0.9985    1.0000   -1.9983    0.9983;
    1.0000   -1.9980    0.9981    1.0000   -1.9978    0.9978;
    1.0000   -1.9975    0.9976    1.0000   -1.9972    0.9973;
    1.0000   -1.9967    0.9968    1.0000   -1.9966    0.9967;
    1.0000   -1.9959    0.9960    1.0000   -1.9957    0.9958;
    1.0000   -1.9949    0.9951    1.0000   -1.9944    0.9946;
    1.0000   -1.9936    0.9939    1.0000   -1.9928    0.9931;
    1.0000   -1.9915    0.9920    1.0000   -1.9911    0.9916;
    1.0000   -1.9893    0.9901    1.0000   -1.9886    0.9894;
    1.0000   -1.9862    0.9874    1.0000   -1.9855    0.9867;
    1.0000   -1.9827    0.9847    1.0000   -1.9807    0.9827;
    1.0000   -1.9773    0.9805    1.0000   -1.9755    0.9786;
    1.0000   -1.9705    0.9755    1.0000   -1.9681    0.9731;
    1.0000   -1.9613    0.9692    1.0000   -1.9583    0.9662;
    1.0000   -1.9487    0.9612    1.0000   -1.9454    0.9579;
    1.0000   -1.9316    0.9514    1.0000   -1.9275    0.9472;
    1.0000   -1.9076    0.9387    1.0000   -1.9033    0.9344;
    1.0000   -1.8745    0.9235    1.0000   -1.8690    0.9179;
    1.0000   -1.8286    0.9054    1.0000   -1.8203    0.8968;
    1.0000   -1.7621    0.8821    1.0000   -1.7524    0.8718;
    1.0000   -1.6674    0.8537    1.0000   -1.6564    0.8415;
    1.0000   -1.5327    0.8200    1.0000   -1.5190    0.8036;
    1.0000   -1.3401    0.7782    1.0000   -1.3258    0.7589;
    1.0000   -1.0683    0.7246    1.0000   -1.0595    0.7098;
    1.0000   -0.7085    0.6810    1.0000   -0.6904    0.6319;
    1.0000   -0.2327    0.6152    1.0000   -0.2309    0.5630;
    1.0000    0.3403    0.5708    1.0000    0.2990    0.4740;
    1.0000    0.9586    0.5445    1.0000    0.8184    0.3951;
    1.0000    1.5774    0.6630    1.0000    1.1537    0.2999];

Answer 1

And so a sample by sample, rather inefficient way of filtering the noise out of a dirac impulse from a FDN would be to add 2 more buffers and means of calculating difference equations of 31 cascaded biquad filters (Any suggestions for improving calculation speed comment below)

% Feedback Delay Network - Mono
% Creates impulse response based on reverberation times

fs = 44100;
in = [ 1; 0 ];            % Dirac Impulse
in = [in; zeros(3*fs,1)]; % Space for reverb
 
% Householder matrix N=16
A = eye(4) - ( 2/4 * ones(4,1) * ones(1,4) );
MTX = .5 * [A -A -A -A; -A A -A -A; -A -A A -A; -A -A -A A];

N = size(MTX,1);    % Matrix order

delays = randi([dmin dmax],N); % N delays

load('GEQ.mat');    % Third octave graphic equalizer calculated based
                    % on an atenuation-per-sample and scaled by delays.
                    % SOS Form Size 31x6xN 
                   
% Initialize signals
delayBuffer = max(delays) + max(delays)/10;
bufferN = zeros(delayBuffer,N);   % Delay buffers
   buffFin = zeros(31,3,N);         % New buffers for filter outputs
   buffFout = zeros(31,3,N);
FB = zeros(1,N);               % Feedback matrix output
in_dl = zeros(1,N);            % Delay input
out_dl = zeros(1,N);           % Delay output
   out_fdl = zeros(1,N);            % Filtered delay output
nSample = length(in);          % Number of samples
out = zeros(nSample,1);        % Output
 
% FDN Computation
for sample = 1:nSample     % each sample
    for n = 1:N            % each delayline
         
       in_dl(n) = in(sample,1) + FB(n); % Input and Feedback matrix sum
       
       % Delaying
       [out_dl(n),bufferN(:,n)] = funcDelay( in_dl(n), bufferN(:,n),...
                           sample, delays(n) );
       % Filtering
       [out_fdl(n),buffFin(:,:,n),buffFout(:,:,n)] = funcFilterGEQ(...
            GEQ(:,:,n), out_dl(n), buffFin(:,:,n), buffFout(:,:,n));  
    end
     
    out(sample,1) = sum(out_fdl); % Filtered delay output sum
     
    FB = out_fdl * MTX; % Feedback matrix output recalculation
end


% Used functions
function [out,buffer] = funcDelay( in, buffer, n, delay)
  
  % Circular buffer indices
  len = length(buffer);
  indexC = mod(n-1, len) + 1;        % Current
  indexD = mod(n-delay-1, len) + 1;  % Delay
  
  out = buffer(indexD,1);
  
  % Stores output on appropriate index
  buffer(indexC,1) = in; 
end

function [outfilt,buffin,buffout] = funcFilterGEQ( GEQ, in, buffin, buffout)
                          % Sample based filter cascading
nBands = size(GEQ,1);

out = zeros(1,nBands+1);
out(1) = in;              % Stores value pre-filtered
                          
for b = 1:nBands          % Performs series bandpass filtering
    [out(b+1),buffin(b,:),buffout(b,:)] = funcBandpassFilt( out(b),...
        GEQ(b,1:3), GEQ(b,5:6), buffin(b,:), buffout(b,:) );
end

outfilt = out(end);       % Outputs final value post-filtered
end

function [out,buffin,buffout] = funcBandpassFilt( in, bs, as, buffin, buffout)
                          % Sample based filtering
buffin(3) = buffin(2);    % Valid for biquad filters
buffin(2) = buffin(1);    
buffin(1) = in;           % Sequential indexing

buffout(1) = bs(1)*buffin(1) + bs(2)*buffin(2) + bs(3)*buffin(3)...
                                - as(1)*buffout(2) - as(2)*buffout(3);
out = buffout(1);         
                          % Outputs calculation based on 3 latest values
buffout(3) = buffout(2);
buffout(2) = buffout(1);  
buffout(1) = 0;
end