How do I implement a low pass filter difference equation to filter sample by sample, samplewise

Hi Guys, Im currently working on a reverberation alogrithm.

Heres the structure:

My supervisor told me that I have to implement a difference equation to filter the samples samplewise by just putting the samples in the equation and feed the result back to the system etc...

Now my question is how and where in the loop I can do that? The Loop Filters in the picture should be all second order low pass filters. How do I get the difference equation for such a filter?

Heres my code:

in = [ 1; 0 ]; % Dirac Impulse
Fs = 44100;
in = [in; zeros(3*Fs,1)]; % Space for reverb
% Define delay parameters
       %Hadamard Matrix as feedback matrix
H = 0.75*hadamard(16);



delayTime = [0.02, 0.03, 0.05, 0.07];% in seconds 0.01 = 10ms 0.02 = 20ms 0.03 = ....
%feedbackGain = [0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65]; % adjust to control feedback level 0.5,0.53,0.55,0.58,0.6,0.63,0.65,0.68,0.7,0.72,0.75,0.77,0.8,0.81,0.85,0.88
feedbackGain =[0.1,-0.15,-0.2,0.25,-0.3,0.35,0.4,-0.45,-0.5,0.55,0.6,-0.65,-0.7,-0.75,-0.8,0.85]; %0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85
wetGain = [0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7]; % adjust to control wet signal level
% Convert delay time to samples
delaySamples = [101,233, 439,613, 853,1163, 1453,1667 , 1801,2053, 2269,2647, 3001,3607, 4153,4591]; %0.02s, 0.03, 0.04, 0.05 881, 1009, 1321,1511, 1777, 1999, 2203, 2351
                                                               %101, 439, 853, 1453 , 1801, 2269, 3001, 4153

%Summe delays = 30644


% delaySamples = [1200, 1800, 2200,2500];
% Create an empty array to store delayed audio
delayedAudio = zeros(size(in));
delayedAudio2 = zeros(size(in));
delayedAudio3 = zeros(size(in));
delayedAudio4 = zeros(size(in));
delayedAudio5 = zeros(size(in));
delayedAudio6 = zeros(size(in));
delayedAudio7 = zeros(size(in));
delayedAudio8 = zeros(size(in));
delayedAudio9 = zeros(size(in));
delayedAudio10 = zeros(size(in));
delayedAudio11 = zeros(size(in));
delayedAudio12 = zeros(size(in));
delayedAudio13 = zeros(size(in));
delayedAudio14 = zeros(size(in));
delayedAudio15 = zeros(size(in));
delayedAudio16 = zeros(size(in));




delayedAudios =[delayedAudio,delayedAudio2,delayedAudio3,delayedAudio4,delayedAudio5,delayedAudio6,delayedAudio7,delayedAudio8,delayedAudio9,delayedAudio10,delayedAudio11,delayedAudio12,delayedAudio13,delayedAudio14,delayedAudio15,delayedAudio16];
delayedAudios(1,:) = 1;
d = designfilt('lowpassfir','FilterOrder',2 ,'CutoffFrequency',2000,'SampleRate', Fs);

% Apply delay effect

for i = 1:length(in)-2*max(delaySamples)
     for k = 1:16
        if delayedAudios(i+delaySamples(k),k) == 0
        delayedAudios(i+delaySamples(k),k) = filterfeedbackGain(k)*delayedAudios(i,k);
        end
        
     end      

     for k = 1:16
        for j = 1:16
           if delayedAudios(i+delaySamples(j)+delaySamples(k),k) == 0 && k~= j 
            delayedAudios(i+delaySamples(j)+delaySamples(k),k) = H(j,k)*delayedAudios(i+delaySamples(j),j);
           end
        end
     end     
end


%Original and Mix together

for k = 1:16
    delayedAudios(:,k) = wetGain(k)*delayedAudios(:,k);
end

%outputAudioM = delayedAudios(:,16) + delayedAudios(:,15) + delayedAudios(:,14) + delayedAudios(:,13) + delayedAudios(:,12) +delayedAudios(:,11) +delayedAudios(:,10) +delayedAudios(:,9)+delayedAudios(:,8) +delayedAudios(:,7) +delayedAudios(:,6) +delayedAudios(:,5) +delayedAudios(:,4) +delayedAudios(:,3) +delayedAudios(:,2 )+delayedAudios(:,1);
outputAudioM = sum(delayedAudios,2);
outputAudioM(1) = outputAudioM(1)*0.02;




% %Write the output audio to a new file and plot
audiowrite('output_audio_with_delay.wav', outputAudioM, Fs);
plot(outputAudioM);

Answer 1

[Mathieu NOE]

hello

see this second order recursion - here used with a step input (works for any input signal)

I used the b,a coefficients of a second order low pass Butterworth filter (up to you to change filter type and cut off frequency)

all this can by simply implemented using gain and delay blocks in simulink

[b,a] = butter(2,0.1);
b0 = b(1);
b1 = b(2);
b2 = b(3);
a1 = a(2);
a2 = a(3);

% step response (step input at sample = 11) 
% manual for loop coding IIR filter
x = [zeros(10,1); ones(60,1)];
samples = length(x);
y(1) = b0*x(1) + 0 + 0  + 0 + 0;
y(2) = b0*x(2) + b1*x(1) + 0  - a1*y(1) - 0;
for k = 3:samples 
    y(k) = b0*x(k) + b1*x(k-1) + b2*x(k-2)  - a1*y(k-1) - a2*y(k-2);
end

figure(1)
plot(x,'k')
hold on
plot(y,'r')

[Muhsin Zerey]

Hi Mathieu, thanks for your answer. Ive seen you manually assigned y(1) and y(2). Do you mind explain to me why you did that and does it has a relevance to my code?

[Mathieu NOE]

hello

simply because matlab cannot work with zero or negative index , so if you have a recursion like above

y(k) = b0*x(k) + b1*x(k-1) + b2*x(k-2) - a1*y(k-1) - a2*y(k-2);

you must take care of how to generate the first 2 outputs (as for k = 1 and 2 , k - 2 is <= 0 )

[Muhsin Zerey]

Thanks, that makes sense to me now.

In my code all my samples are delayed and each with a different delayTime. Can I just use your code to filter sample after sample starting at 1 or do I have to consider the delaytime too? What is your opinion?

[Mathieu NOE]

the recursion is valid for any time delay , but if you keep the same a and b coefficients , this would correspond to a different cut off frequency of your digital filter

if you need to have the same exact filtering for different sampling rates , then you have to compute a and b coefficients for each sampling rate and use them accordingly on the corresponding data

[Muhsin Zerey]

Thank you very much. That is really helpful.

But now I struggle a bit with implementing that transfer function to my algorithmus since I dont have input x and output y. I only have a 16x16 matrix called delayedAudios where everything is stored. Do you know how I can implement that butterworth filter to my code?

[Mathieu NOE]

the same logic, with minor modifications, can be used for a multichannel x input

here I generated 16 steps , with some time delta between the channels

clc
clearvars

[b,a] = butter(2,0.1);
b0 = b(1);
b1 = b(2);
b2 = b(3);
a1 = a(2);
a2 = a(3);
% step response (step input at sample = 11) 
% manual for loop coding IIR filter

% multi channel input (70 samples / 16 channels
channels = 16;
for ci = 1:channels
    x(:,ci) = [zeros(10+ci,1); ones(60-ci,1)];
end

% filtering the multi channel input
samples = size(x,1);
y(1,:) = b0*x(1,:) + 0 + 0  + 0 + 0;
y(2,:) = b0*x(2,:) + b1*x(1,:) + 0  - a1*y(1,:) - 0;
for k = 3:samples 
    y(k,:) = b0*x(k,:) + b1*x(k-1,:) + b2*x(k-2,:)  - a1*y(k-1,:) - a2*y(k-2,:);
end

figure(1)
plot(x,'k')
hold on
plot(y,'r')