为了提供一些上下文，我试图在matlab中从头开始实现以下信号图，即反馈延迟网络(FDN(。图片：FDN
通过适当的矩阵，与延迟长度无关，当输入dirac脉冲时，实际上会产生白噪声。
我已经在代码中做到了这一点，但我的目标是另一个，因此也是我的问题。我想在每个延迟线z^-m之后应用滤波器h(z(。图片：h(z(
更具体地说，我想在每个延迟线之后应用第三倍频程级联图形均衡器。其目的是在整个结构上产生与频率相关的衰减，从而产生与延迟相关的衰减。我已经成功地以SOS的形式设计了过滤器，但我的问题是：如何在结构中应用它？我想在某个地方使用sosfilt((，但我不确定。
我没有为了目的而降低系统的顺序。顺序为16(16x16矩阵、16条延迟线、31x16双四元滤波器(
第一个代码是指无损FDN，它可以安全运行，产生白噪声。我评论了我在循环中引入过滤的失败尝试，说：%过滤
不幸的是，我不能发布所有GEQ条目，但我会在最后留下8个对应于前8个延迟的条目。
因此，问题是如何编码来过滤白噪声，在整个FDN结构中实现频率相关衰减。此外，尽管它可能在计算上效率低下，但我更喜欢在没有更高级别函数的情况下应用它，并基于我已经拥有的功能，即：适用于GNU Octave
编辑：假设你必须使用差分方程逐个样本地应用带通二阶滤波，你将如何递归地对31个串联频带进行滤波？其中一个显示在第二个代码部分。
% 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)
第二个代码部分：将滤波器差分方程应用于信号。建议递归地为31个过滤器编码？
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
谢谢！
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];
因此，从FDN的dirac脉冲中过滤噪声的一种逐样本、效率相当低的方法是再添加2个缓冲器和计算31个级联双二阶滤波器差分方程的方法(下面有任何提高计算速度的建议评论(
% 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