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];
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